*Article* **Simulating Extraocular Muscle Dynamics. A Comparison between Dynamic Implicit and Explicit Finite Element Methods**

**Jorge Grasa 1,2,\* and Begoña Calvo 1,2**


**Abstract:** The finite element method has been widely used to investigate the mechanical behavior of biological tissues. When analyzing these particular materials subjected to dynamic requests, time integration algorithms should be considered to incorporate the inertial effects. These algorithms can be classified as implicit or explicit. Although both algorithms have been used in different scenarios, a comparative study of the outcomes of both methods is important to determine the performance of a model used to simulate the active contraction of the skeletal muscle tissue. In this work, dynamic implicit and dynamic explicit solutions are presented for the movement of the eye ball induced by the extraocular muscles. Aspects such as stability, computational time and the influence of mass-scaling for the explicit formulation were assessed using ABAQUS software. Both strategies produced similar results regarding range of movement of the eye ball, total deformation and kinetic energy. Using the implicit dynamic formulation, an important amount of computational time reduction is achieved. Although mass-scaling can reduce the simulation time, the dynamic contraction of the muscle is drastically altered.

**Keywords:** finite element method; implicit FEM; explicit FEM; skeletal muscle

#### **1. Introduction**

The extraocular muscles (EOM) are responsible for the eye movements of the upper eyelid and the eyeball. The group that controls the eyeball contains six muscles: four muscles that run almost a straight course from origin to insertion and hence are called recti and two muscles that run a diagonal course, the oblique muscles [1]. The group that controls eye movement in the cardinal directions are the superior (responsible for elevation, incyclotorsion and adduction), inferior (responsible for depression, extorsion (outward, rotational movement) and adduction), lateral (responsible for abduction) and medial (responsible for adduction) rectus muscles. The movements of the extraocular muscles take place under the influence of a system of extraocular soft tissue pulleys in the orbit. The extraocular muscle pulley system is fundamental to the movement of the eye muscles, in particular to ensure conformity to Listing's law. Certain diseases of the pulleys (heterotopy, instability and hindrance of the pulleys) cause particular patterns of incomitant strabismus [2]. Simulating and analyzing eye movements is useful for assessing the role of these tissues and for exploring the equilibrium of the applied forces that can be impaired and lead to different pathologies [3,4].

The finite element method (FEM) has been widely used to simulate the behavior of skeletal muscle both passively and actively [5–10]. The vast majority of studies have focused on determining the essential parameters that best fit the experimental evidence [4,10,11]. Although different approximations and scenarios have been evaluated with the help of this numerical technique, the contraction of the muscle has been analyzed assuming quasi-static conditions with no inertia effects [4].

**Citation:** Grasa, J.; Calvo, B. Simulating Extraocular Muscle Dynamics. A Comparison between Dynamic Implicit and Explicit Finite Element Methods. *Mathematics* **2021**, *9*, 1024. https://doi.org/10.3390/ math9091024

Academic Editor: Vince Grolmusz

Received: 22 March 2021 Accepted: 29 April 2021 Published: 1 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

An important aspect to consider when reproducing the eyeball movements using FEM is the fast response of the muscles when activated and how the result of this contraction is translated to induce the system motion. This movement is achieved in a few tenths of a second [11]. When tackling these small time periods, realistic simulations should account for the inertia effect of the mass of the different elements. In such scenarios, the use of a dynamic formulation of the FEM is essential. Time integration algorithms for dynamic problems in FEM analysis can be classified as implicit or explicit. Basically, the implicit method computes the state of the model at each time increment based on the information of that same time increment and the previous time increment, while the explicit method uses the data of the previous time increment to solve the motion equations during the new time increment [12]. Implicit time integration schemes are more expensive per time step, but can obtain the solution for larger time steps and provide a control on the dynamic residual force vector, since they are usually used in conjunction with an iterative procedure within each time step [13]. The explicit algorithm can be solved directly without requiring iteration. This method is conditionally stable, and the critical time step for the operator (without damping) is a function of the material specification and the size of the smallest element of the mesh [12]. Increments larger than this critical time can cause the solution to be unstable and oscillations to occur in the model's response, which can lead to excessively distorted elements. To decrease the total computational time, a mass-scaling technique is commonly used, whereby the solver can use larger time increments by artificially increasing the density of the system. However, it is important to ensure that the added mass does not change the behavior of the model, which in simulating active muscle contraction is decisive. The choice between implicit and explicit methods with or without mass scaling has been the subject of many studies [12,14–16]

The aim of this study was therefore to compare dynamic implicit and explicit solutions using ABAQUS software [17] in the analysis of the contraction of the EOM for eyeball movements. More specifically, we compared the prediction of dynamic effects, potential convergence problems, the accuracy and stability of the calculations, the computational time between the two methods and the influence of mass-scaling in the explicit formulation.

The rest of the paper is organized as follows. Section 2 describes the formulation of the skeletal muscle behavior, the implicit and explicit algorithms, the implementation of the user subroutine and the description of the finite element model. In Section 3, selected results of the model comparing both algorithms are presented and then discussed in Section 4. Finally, in Section 5, the main conclusions are summarized.

#### **2. Materials and Methods**

#### *2.1. Muscle Contraction Model*

Let Ω<sup>0</sup> be a three-dimensional portion of the space representing the solid initial geometry of the muscle tissue. This region defines a set of points at a reference configuration which are identified by the position vector **X**. The motion of the solid defines the current configuration at time *t*, Ω*<sup>t</sup>* and can be described by the map **x** = *χ*(**X**) = **x**(**X**).

Let **F** = *<sup>∂</sup>***<sup>x</sup>** *<sup>∂</sup>***<sup>X</sup>** be the deformation gradient associated with the motion, where *J* ≡ *det***F** > 0 is the Jacobian of the transformation. A multiplicative decomposition of this deformation gradient into volume-changing and volume-preserving parts is established to handle the quasi-incompressibility constraint presented by soft tissues (*J* ∼= 1) [8]

$$\mathbf{F} = J^{\dagger} \mathbf{F}, \qquad \mathbf{F} = J^{-\dagger} \mathbf{F} \tag{1}$$

$$\mathbf{C} = \mathbf{F}^T \mathbf{F}, \qquad \mathbf{\tilde{C}} = J^{-\frac{2}{3}} \mathbf{C} = \mathbf{\tilde{F}}^T \mathbf{\tilde{F}} \tag{2}$$

$$\mathbf{b} = \mathbf{F} \mathbf{F}^T, \qquad \mathbf{\tilde{b}} = \mathbf{J}^{-\frac{2}{3}} \mathbf{b} = \mathbf{F} \mathbf{F}^T \tag{3}$$

where *J* 1 <sup>3</sup> **I** and **F¯** represent the volumetric and deviatoric deformation gradients, respectively. **C** and **b** are the right and left Cauchy–Green strain tensors and **C¯** and **b¯** their modified counterparts.

It is assumed in this work that the contraction process can be modelled as two fictitious steps [8,18] (see Figure 1). The first step is associated with the relative motion of the protein filaments myosin and actin during the power stroke of the cross bridges, and the second step relates to the elastic deformation of cross bridges. This contraction process can be expressed as a multiplicative decomposition of the deformation gradient **F¯**:

$$
\mathbf{F} = \mathbf{F}\_c \mathbf{F}\_a \tag{4}
$$

where **F¯** *<sup>a</sup>* is the deformation gradient associated with the contractile response induced by the actin and myosin translation, whereas **F¯***<sup>e</sup>* defines a deformation due to the cross bridges elasticity. The gradient **F¯** *<sup>a</sup>* represents the active contraction so it does not need to be integrable. Thus, infinitesimal parts of the tangent space Ω<sup>0</sup> are deformed independently, and the configuration they form after the motion may not be compatible. The gradient **F¯***<sup>e</sup>* guarantees the compatibility in the deformed configuration Ω*t*. Accordingly, let **C¯** *<sup>e</sup>* = **F¯** *<sup>T</sup> <sup>e</sup>* **F¯***<sup>e</sup>* = **F¯** <sup>−</sup>*<sup>T</sup> <sup>a</sup>* **C¯ F¯** <sup>−</sup><sup>1</sup> *<sup>a</sup>* be a deformation measure due to the titin and cross bridges motion which is not a state variable since it depends on **C¯** and **F¯** *<sup>a</sup>*.

**Figure 1.** Illustration of the contraction process modelled as two fictitious steps. **G***<sup>i</sup>* vectors located at a point **X** in a muscle fiber in the initial configuration Ω<sup>0</sup> transform into new vectors *λ*¯ *<sup>i</sup>***G***<sup>i</sup>* by the active contraction **F***a*. The intermediate step is associated with the relative motion of the protein filaments due to cross bridges power stroke. In the final step, the cross-bridges are deformed by **F***<sup>e</sup>* to restore the compatibility of the deformed configuration Ω*t*.

The active contraction occurs along the direction that is defined by the muscle fibers, so let us introduce this direction as **n**0, and let *λ*¯ *<sup>a</sup>* be the active stretch. Thus, the active contractile tensor, **F¯** *<sup>a</sup>*, can be written in the local coordinate system, **G***i*, as:

$$\bar{\mathbf{F}}'\_{a} = \bar{\lambda}\_{a} \mathbf{G}\_{1} \otimes \mathbf{G}\_{1} + \bar{\lambda}\_{a}^{-1/2} \mathbf{G}\_{2} \otimes \mathbf{G}\_{2} + \bar{\lambda}\_{a}^{-1/2} \mathbf{G}\_{3} \otimes \mathbf{G}\_{3} \tag{5}$$

where we assume the active contractile tensor **F¯** *<sup>a</sup>* to be isochoric. This local coordinate system varies along the fiber length and is represented at a particular point of the tissue in Figure 1.

The components of the contractile tensor expressed in the global system of coordinates **E***i*, **F¯** *<sup>a</sup>*, can be obtained as:

$$\mathbf{F}\_a = \mathbf{R}^T \mathbf{F}\_a^\prime \mathbf{R} \tag{6}$$

where **R** is the rotation tensor.

A strain energy density formulation decoupled into volume-changing and volumepreserving parts is commonly taken to formulate the elastic constitutive law for transversely isotropic materials such as skeletal muscle [8,18,19]. This energy is formulated in this work as:

$$\Psi = \Psi\_{\rm vol}(f) + \Psi\_p(\mathbf{\tilde{C}}, \mathbf{N}) + \Psi\_a(\mathbf{\tilde{C}}\_a, \mathbf{\tilde{A}}\_a, \mathbf{N}) \tag{7}$$

where **N** = **n**<sup>0</sup> ⊗ **n**0. Equivalently, Ψ can be expressed as a function of the invariants of the strain tensors:

$$\Psi = \Psi\_{\text{vol}}(I) + \Psi\_p(I\_1, I\_2, I\_4) + \Psi\_a(I\_{4\prime}X\_a) \tag{8}$$

where ¯*I*<sup>1</sup> = *tr***C¯** and ¯*I*<sup>2</sup> = <sup>1</sup> <sup>2</sup> ((*tr***C¯** )<sup>2</sup> <sup>−</sup> *tr***C¯** <sup>2</sup>) are the first and second modified strain invariants of the symmetric modified Cauchy–Green tensor **C¯** , and ¯*I*<sup>4</sup> <sup>=</sup> **<sup>n</sup>**<sup>0</sup> · **Cn¯** <sup>0</sup> <sup>=</sup> *<sup>λ</sup>*¯ <sup>2</sup> is the pseudo-invariant related to the anisotropy of the passive response (collagen fibers). Similarly, the active contribution of the strain energy function, Ψ¯ *<sup>a</sup>*, is expressed in terms of the pseudo-invariant associated to **C¯** *<sup>e</sup>* and the preferred direction **<sup>n</sup>**0, ¯*J*<sup>4</sup> <sup>=</sup> **<sup>n</sup>**<sup>0</sup> · **C¯** *<sup>e</sup>***n**<sup>0</sup> <sup>=</sup> *<sup>λ</sup>*<sup>2</sup> *e*. As shown, the anisotropy in the formulation is induced by a single orientation for both the passive and the active behavior. Although in fusiform muscles such as the EOM this is commonly accepted, in other muscle architectures two families of fibers should be considered to adopt a more suitable formulation [9,20].

The third term in Equation (8) represents the strain energy associated with the active response and, consequently, with the actin-myosin interaction. This term is written here as a function Ψ¯ *<sup>a</sup>* that relates to the energy stored in the cross-bridges, while *f*1(*λ*¯ *<sup>a</sup>*) is a function that accounts for the filament overlap and *f*2(*t*) for the muscle activation level:

$$\Psi = \Psi\_{vol}(f) + \bar{\Psi}\_p(\bar{I}\_1, \bar{I}\_2, \bar{I}\_4) + f\_1(\bar{\lambda}\_4) f\_2(t) \bar{\Psi}\_a'(\bar{f}\_4) \tag{9}$$

The function 0 < *f*1(*λ*¯ *<sup>a</sup>*) < 1 has been experimentally characterized in previous studies for different muscles [19,21] and fitted by a smooth exponential relationship:

$$f\_1(\bar{\lambda}\_4) = \exp^{\frac{-(\bar{\lambda}\_4 - \lambda\_{opt})^2}{2\xi^2}}\tag{10}$$

where *λopt* is the optimum length of the muscle at which isometric maximum stress is developed and *ξ* determines the curvature of the function. To formulate *f*2(*t*), we assume in this work that all muscle fibers are completely recruited in each contraction, and this function can be expressed as:

$$f\_2(t) = a \tanh^2(s\_1 t) \tag{11}$$

where 0 < *α* < 1 governs the activation level, *s*<sup>1</sup> regulates the initial slope of the function and *t* is the time variable. Finally, the energy stored in the cross-bridges is expressed in terms of the invariant associated to **C¯** *<sup>e</sup>* in the direction of the muscle fibers **n**<sup>0</sup> and a parameter *P*<sup>0</sup> related to the maximum active stress:

$$\Psi\_a' = \frac{1}{2} P\_0 (\bar{f}\_4 - 1)^2 \tag{12}$$

During the muscle contraction process, the second law of thermodynamics can be formulated in the shape of the Clausius–Planck inequality neglecting the thermal dissipation

rate. This inequality allows us to consider that some of the power produced internally is stored while another portion is dissipated [8]:

$$\mathcal{D}\_{\text{int}} = -\Psi + \frac{1}{2}\mathbf{S} : \dot{\mathbf{C}} + \frac{1}{2}\mathbf{S}\_a : \dot{\mathbf{C}}\_a \ge 0 \tag{13}$$

In Equation (13), **S***<sup>a</sup>* represents active stress and <sup>1</sup> <sup>2</sup>**S***<sup>a</sup>* : **<sup>C</sup>**˙ *<sup>a</sup>* the muscle power stroke [18]. Following the work of Hernández-Gascón et al. [8], the following constitutive relations are obtained:

$$\mathbf{S} = 2\frac{\partial \Psi}{\partial \mathbf{C}} + \mathbf{F}\_a^{-1} \left( 2\frac{\partial \Psi}{\partial \mathbf{C}\_\varepsilon} \right) \mathbf{F}\_a^{-T} \tag{14}$$

$$\left(\mathbf{P}\_d - 2\mathbf{F}\_d \frac{\partial \Psi}{\partial \mathbf{C}\_d} + 2\mathbf{C}\_\varepsilon \frac{\partial \Psi}{\partial \mathbf{C}\_\varepsilon} \mathbf{F}\_d^{-T}\right) : \mathbb{F}\_d \ge 0 \tag{15}$$

where **P***a* is the first Piola–Kirchoff active stress. Since contraction occurs along the muscle fiber only, Equation (15) reduces to:

$$\left[P\_a - \frac{\partial \Psi}{\partial \bar{\lambda}\_a} + \left(2\mathbf{C}\_c \frac{\partial \Psi}{\partial \mathbf{C}\_c} \mathbf{F}\_a^{-T}\right) : \frac{\partial \mathbf{F}\_a}{\partial \bar{\lambda}\_a}\right] \dot{\lambda}\_a \ge 0 \tag{16}$$

This expression leads to the following constitutive relation for the active contraction velocity ¯˙ *λa*:

$$P\_a - \frac{\partial \Psi}{\partial \vec{\lambda}\_a} + \left(2\vec{\mathbf{C}}\_{\epsilon} \frac{\partial \Psi}{\partial \vec{\mathbf{C}}\_{\epsilon}} \vec{\mathbf{F}}\_a^{-T}\right) : \frac{\partial \mathbf{F}\_a}{\partial \vec{\lambda}\_a} = \mathbb{C}\vec{\lambda}\_a \tag{17}$$

assuming that:

$$C = \frac{1}{v\_0} P\_0 f\_1(\lambda\_a) f\_2(t) \tag{18}$$

where *v*<sup>0</sup> is associated with the initial contraction velocity. The active stress *Pa* from Equation (17) is defined as a function of *P*0, *f*1(*λa*) and *f*2(*t*), and *ν* is a friction parameter that takes into account the relative sliding speed between actin and myosin:

$$P\_a = -\nu P\_0 f\_1(\vec{\lambda}\_a) f\_2(t) \tag{19}$$

Substituting Equations (18) and (19) and the last term of Equation (9) into Equation (17) leads to the expression for the contraction velocity:

$$\lambda\_a = v\_0 \left[ -\nu - \frac{1}{f\_1(\lambda\_a)} \frac{\partial f\_1(\lambda\_a)}{\partial \lambda\_a} \Psi\_a'(f\_4) + \left( 2\mathbf{C}\_\varepsilon \frac{\partial \Psi\_a'(f\_4)}{\partial \mathbf{C}\_\varepsilon} \mathbf{F}\_a^{-T} \right) : \frac{\partial \mathbf{F}\_a}{\partial \lambda\_a} \right] \tag{20}$$

Since Ψ¯ *<sup>a</sup>* depends on ¯*J*4, Equation (20) reduces to:

$$\vec{\lambda}\_a = \upsilon\_0 \left[ -\nu - \frac{1}{f\_1(\vec{\lambda}\_a)} \frac{\partial f\_1(\vec{\lambda}\_a)}{\partial \vec{\lambda}\_a} \Psi\_a'(f\_4) + 2 \frac{\lambda\_a^2}{\vec{\lambda}\_a} \frac{\partial \Psi\_a'(f\_4)}{\partial f\_4} \right] \tag{21}$$

Taking the first constitutive relation (Equation (14)) and the particular form of the strain energy density function, the expressions for the Cauchy stress tensor and the elasticity tensor can be derived. Both tensors must be provided to define the mechanical constitutive model in the implicit user material subroutine in ABAQUS/Standard [17], whereas only the definition of the Cauchy stress is needed in ABAQUS/Explicit [17].

From Equations (7) and (14), the second Piola–Kirchhoff stress tensor is found to be:

$$\mathbf{S} = \mathbf{S}\_{\text{vol}} + \mathbf{S}\_p + \mathbf{S}\_a \tag{22}$$

The Cauchy stress tensor is obtained by means of a weighted push-forward operation of **<sup>S</sup>**, *<sup>σ</sup>* <sup>=</sup> *<sup>J</sup>*−1*χ*∗(**S**) = *<sup>J</sup>*−1**FSF***T*:

$$\begin{split} \sigma &= \sigma\_{\text{vol}} + \bar{\sigma}\_{\text{p}} + \bar{\sigma}\_{\text{a}} = p\mathbf{1} + \frac{1}{\bar{\mathcal{I}}} dev\left[\mathbf{F}\tilde{\mathbf{S}}\_{\text{p}}\mathbf{F}^{\text{T}}\right] + \frac{1}{\bar{\mathcal{I}}} dev\left[\mathbf{F}\_{\text{e}}\tilde{\mathbf{S}}\_{\text{e}}\mathbf{F}\_{\text{e}}^{\text{T}}\right] \\ &= p\mathbf{1} + \frac{1}{\bar{\mathcal{I}}} dev[\tilde{\sigma}\_{\text{p}}] + \frac{1}{\bar{\mathcal{I}}} dev[\tilde{\sigma}\_{\text{e}}] \end{split} \tag{23}$$

with:

$$p = \frac{d\Psi\_{vol}(f)}{d\mathbf{J}}, \qquad d\boldsymbol{\varepsilon}\boldsymbol{v}[\cdot] = (\cdot) - \frac{1}{\mathfrak{J}}\boldsymbol{tr}[\cdot]\mathbf{1} \tag{24}$$

Differentiating Equation (22) with respect to **C** leads to the material elasticity tensor C, which can be divided into volumetric and deviatoric parts associated with the passive and active responses as follows:

$$\mathbf{C} = \mathbf{C}\_{\text{vol}} + \mathbf{\tilde{C}}\_{p} + \mathbf{\tilde{C}}\_{a} = 2\frac{\partial \mathbf{S}\_{\text{vol}}}{\partial \mathbf{C}} + 2\frac{\partial \mathbf{\tilde{S}}\_{p}}{\partial \mathbf{C}} + 2\frac{\partial \mathbf{\tilde{S}}\_{a}}{\partial \mathbf{C}} \tag{25}$$

The elasticity tensor in the spatial configuration, C, is obtained by a weighted pushforward operation of <sup>C</sup>, which can be expressed as <sup>C</sup> <sup>=</sup> *<sup>J</sup>*−1*χ*∗(C) and results in:

$$\mathbf{c} = \mathbf{c}\_{\text{vol}} + \mathbf{c}\_p + \mathbf{c}\_a \tag{26}$$

For a detailed explanation about these expressions and further information, the reader is referred to the works of Weiss et al. [22] and Hernández-Gascón et al. [8].

#### *2.2. Principle of Virtual Work and Finite Element Discretization*

The principle of virtual work that allows to establish the finite element formulation is derived from the balance of momentum of a body *V* with boundary *S* that can be written as:

$$
\int\_{S} \mathbf{t}dS + \int\_{V} \rho \mathbf{g}dV = \int\_{V} \rho \ddot{\mathbf{a}} \tag{27}
$$

where *t* is the stress vector, *<sup>ρ</sup>* is the material density, *g* is the gravity acceleration and **<sup>a</sup>**¨ is the accelerations vector. Applying the relation between the stress vector and the stress tensor *t* <sup>=</sup> *n<sup>σ</sup>* with *n* the normal surface vector and the divergence theorem [13], the following relation must be fulfilled at each material point:

$$
\nabla \cdot \sigma + \rho \mathbf{g} = \rho \ddot{\mathbf{a}} \tag{28}
$$

Multiplying this equation by a virtual displacement field *<sup>δ</sup>a* and integrating over the domain *<sup>V</sup>*: -

$$\int\_{V} \delta \mathbf{a} (\nabla \cdot \boldsymbol{\sigma} + \rho \mathbf{g} - \rho \ddot{\mathbf{a}}) = 0 \tag{29}$$

After applying the divergence theorem and some manipulations, the weak form of the equation of motion that represents the principle of virtual work is obtained:

$$\int\_{V} \boldsymbol{\sigma} : \delta \boldsymbol{\sigma} dV = \int\_{S} \rho (\mathbf{g} - \ddot{\mathbf{a}}) \delta \mathbf{a} dV + \int\_{S} \mathbf{t} \delta \mathbf{a} dS \tag{30}$$

with *e* = <sup>1</sup> <sup>2</sup> (**<sup>1</sup>** <sup>−</sup> *<sup>F</sup>*−*TF*−1) the Euler–Almansi strain tensor.

Equation (30) is the basis for the finite element discretization where the displacements at the nodes of the mesh elements are considered as the fundamental unknowns. The continuous displacement field *a* can be approximated at each element as:

$$\mathfrak{a} = \sum\_{k=1}^{n} \phi\_k(\zeta, \eta, \zeta) \mathfrak{u}\_k \tag{31}$$

where *φ<sup>k</sup>* are the interpolation functions of an element supported by *n* nodes and (*ξ*, *η*, *ζ*) the coordinates of the reference element. The proper definition of the interpolation functions using polynomials (linear in this work) and assembling element matrix contributions to the integration over all the solid volume in Equation (30) conduct to the well-known finite element equation system.

#### *2.3. Implicit Solution Method*

The dynamic response in ABAQUS/Standard for nonlinear models is obtained by direct time integration of all of the degrees of freedom of the finite element model [17]. The time step for implicit integration can be chosen automatically on the basis of the "half-increment residual" by monitoring the values of equilibrium residuals at *t* + Δ*t*/2 once the solution at *t* + Δ*t* has been obtained. The accuracy of the solution can be assessed and the time step adjusted appropriately.

The equilibrium equations are written at the end of the time step (time *t* + Δ*t*) and are calculated from the time integration operator. The finite element approximation is:

$$\mathbf{M}\ddot{\mathbf{u}} + \mathbf{I} = \mathbf{F} \tag{32}$$

where **F** is the vector of externally applied forces, **I** is the vector of internal element forces, **M** is the lumped mass matrix and **u**¨ is the accelerations vector.

The algorithm defined by Hilber et al. [23] is:

$$\mathbf{M}\ddot{\mathbf{u}}^{(i+1)} + (1+a)\mathbf{K}\mathbf{u}^{(i+1)} - a\mathbf{K}\mathbf{u}^{(i)} = (1+a)\mathbf{F}^{(i+1)} - a\mathbf{F}^{(i)}\tag{33}$$

where **u** is the displacement vector and maintaining Newmark's assumption that the acceleration **u**¨ varies linearly over the time step [13]:

$$\mathbf{u}^{(i+1)} = \mathbf{u}^{(i)} + \Delta t \dot{\mathbf{u}}^{(i)} + \Delta t^2 \left[ \left(\frac{1}{2} - \beta\right) \ddot{\mathbf{u}}^{(i)} + \beta \ddot{\mathbf{u}}^{(i+1)} \right] \tag{34}$$

$$\dot{\mathbf{u}}^{(i+1)} = \dot{\mathbf{u}}^{(i)} + \Delta t \left[ (1 - \gamma) \ddot{\mathbf{u}}^{(i)} + \gamma \ddot{\mathbf{u}}^{(i+1)} \right] \tag{35}$$

being **u**˙ the velocities vector and

$$\beta = \frac{1}{4}(1-a)^2; \qquad \gamma = 1/2 - a; \qquad -\frac{1}{3} \le a \le 0 \tag{36}$$

*α* = −0.05 is chosen by default in ABAQUS/Standard as a small damping term to remove the high frequency noise without affecting the lower frequency response [15]. Taking *α* = 0, the Newmark method is obtained. For the algorithmic implementation, it is necessary to obtain **u**¨ (*i*+1) from Equation (34) :

$$
\ddot{\mathbf{u}}^{(i+1)} = \frac{1}{\beta \Delta t^2} \Delta \mathbf{u} - \frac{1}{\beta \Delta t} \dot{\mathbf{u}}^{(i)} - \frac{1 - 2\beta}{2\beta} \ddot{\mathbf{u}}^{(i)} \tag{37}
$$

where <sup>Δ</sup>**<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**(*i*+1) <sup>−</sup> **<sup>u</sup>**(*i*). Substituting this expression into Equation (33) and taking **I**(*i*) = **Ku**(*i*) yields:

$$
\Delta \mathbf{u} = \left(\mathbf{K}^\*\right)^{-1} \mathbf{F}^\* \tag{38}
$$

with the algorithmic tangential stiffness matrix

$$\mathbf{K}^\* = (1+\alpha)\mathbf{K} + \frac{1}{\beta \Delta t^2} \mathbf{M} \tag{39}$$

and the vector **F**∗ defined as:

$$\mathbf{F}^\* = (1+a)\mathbf{F}^{(i+1)} - a\mathbf{F}^{(i)} - \mathbf{I}^{(i)} + \mathbf{M} \left( \frac{1}{\beta \Delta t} \dot{\mathbf{u}}^{(i)} + \frac{1-2\beta}{2\beta} \ddot{\mathbf{u}}^{(i)} \right) \tag{40}$$

Three factors should be considered when selecting the maximum allowable time step size: the rate of variation of the applied loading, the complexity of the nonlinear damping and stiffness properties and the typical period of vibration of the structure [17]. A maximum increment versus period ratio Δ*t*/*T* < 1/10 is recommended for obtaining reliable results [17], where *T* is the period of typical modes of response. The Hilber et al. [23] *α*-method time integration scheme for solving the implicit problem can be summarized as:

	- (a) Initialize the displacement increment Δ**u**<sup>0</sup> and the internal force **I** (*i*+1) <sup>0</sup> <sup>=</sup> **<sup>I</sup>**(*i*)
	- (b) Iterations *j* = 0, ... for finding "dynamic equilibrium" within the time step increment:
		- Compute tangential stiffness matrix: **K***<sup>j</sup>*
		- Compute the algorithmic stiffness matrix: **K**∗ *<sup>j</sup>* = (<sup>1</sup> <sup>+</sup> *<sup>α</sup>*)**K***<sup>j</sup>* <sup>+</sup> <sup>1</sup> *<sup>β</sup>*Δ*t*<sup>2</sup> **M**
		- Compute **F**∗ *<sup>j</sup>* = (<sup>1</sup> <sup>+</sup> *<sup>α</sup>*)**F**(*i*+1) <sup>−</sup> *<sup>α</sup>***F**(*i*) <sup>−</sup> **<sup>I</sup>** (*i*) *<sup>j</sup>* + **M** 1 *<sup>β</sup>*Δ*t***u**˙ (*i*) <sup>+</sup> <sup>1</sup>−2*<sup>β</sup>* <sup>2</sup>*<sup>β</sup>* **<sup>u</sup>**¨ (*i*)

*j*

	- **–** Compute the strain increment: Δ**u***j*+<sup>1</sup> → Δ*εk*,*j*+<sup>1</sup>
	- **–** Compute the stress increment: Δ*εk*,*j*+<sup>1</sup> → Δ*σk*,*j*+<sup>1</sup>
	- **–** Compute the total stress: *<sup>σ</sup>k*,*j*+<sup>1</sup> <sup>=</sup> *<sup>σ</sup>*(*i*) *<sup>k</sup>* + Δ*σk*,*j*+<sup>1</sup>

$$\mathbf{\dot{\color{red}{Velocities:}}} \mathbf{\dot{u}}^{(i+1)} = \dot{\mathbf{u}}^{(i)} + \Delta t \Big[ (1 - \gamma) \ddot{\mathbf{u}}^{(i)} + \gamma \ddot{\mathbf{u}}^{(i+1)} \Big]$$


#### *2.4. Explicit Solution Method*

The explicit dynamics analysis procedure in ABAQUS/Explicit is established by an explicit integration rule together with the use of diagonal or "lumped" element mass matrices. The equations of motion for the body are integrated using the explicit central difference integration rule as follows:

$$\mathbf{u}^{(i+1)} = \mathbf{u}^{(i)} + \Delta t^{(i+1)} \dot{\mathbf{u}}^{(i+\frac{1}{2})} \tag{41}$$

$$\dot{\mathbf{u}}^{(i+\frac{1}{2})} = \dot{\mathbf{u}}^{(i-\frac{1}{2})} + \frac{\Delta t^{(i+1)} + \Delta t^{(i)}}{2}\ddot{\mathbf{u}}^{(i)} \tag{42}$$

The superscripts *<sup>i</sup>*, *<sup>i</sup>* <sup>−</sup> <sup>1</sup> <sup>2</sup> and *<sup>i</sup>* <sup>+</sup> <sup>1</sup> <sup>2</sup> refer to the increment number and mid-increment numbers. The state of the analysis is advanced by assuming constant values for the velocities, **u**˙ , and the accelerations, **u**¨, across half time intervals [16]. The accelerations are computed at the start of the increment by:

$$
\ddot{\mathbf{u}}^{(i)} = \mathbf{M}^{-1} (\mathbf{F}^{(i)} - \mathbf{I}^{(i)}) \tag{43}
$$

As the lumped mass matrix is diagonalized, it is a trivial process to invert it, unlike the global stiffness matrix in the implicit solution method. Therefore, each time increment is computationally inexpensive to solve. Several possibilities of this lumping process are available, such as nodal quadrature, row-sum lumping, or a "special lumping technique" where only the latter method produces positive lumped masses for any element type [13]. A stability limit determines the size of the time increment:

$$
\Delta t \le \frac{2}{\omega\_{\max}}\tag{44}
$$

where *ωmax* is the maximum element eigenvalue. A conservative and practical method of implementing the above inequality is:

$$
\Delta t = \min \frac{L^{\varepsilon}}{c^{d}} \tag{45}
$$

where *L<sup>e</sup>* is the characteristic element length and *c<sup>d</sup>* is the dilatational wave speed:

$$
\varepsilon^d = \sqrt{\frac{\lambda + 2\mu}{\rho}}\tag{46}
$$

*λ* and *μ* are the Lamé elastic constants and *ρ* is the material density.

To maintain efficiency of the analysis, it is important to ensure that the sizes of the elements are as regular as possible since one small element reduces the time increment for the whole model [16]. If the inertia effects in the model are negligible or can be considered as quasi-static, it is not useful to maintain the stable time increment as the simulation would take too long. One method that can be used to artificially reduce the runtime of the simulation involves scaling the density of the material in the model. According to Equations (45) and (46), when the density is scaled by a factor, *f* 2, the runtime is reduced by a factor *f* . ABAQUS introduces in the explicit solver a bulk viscosity parameter, which introduces damping associated with the volumetric straining to improve the high speed dynamics simulations [17]. Two types of bulk viscosity parameter are considered: the linear set to 0.03 in this study and the quadratic set to 1.2.

The central difference time integration scheme for non-linear problems employed in most explicit computer codes [13] can be resumed in the following steps:

	- (a) Solve for total displacements: **u**(*i*+1) = **u**(*i*) + Δ*t* (*i*+1)**u**˙ (*i*<sup>+</sup> <sup>1</sup> 2 )
	- (b) Compute the displacement increment: <sup>Δ</sup>**<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**(*i*+1) <sup>−</sup> **<sup>u</sup>**(*i*)
	- (c) For each integration point j:
		- Compute the strain increment: Δ**u** → Δ*ε<sup>j</sup>*
		- Compute the stress increment: Δ*ε<sup>j</sup>* → Δ*σ<sup>j</sup>*
		- Compute the total stress: *σi*+<sup>1</sup> *<sup>j</sup>* = *σ<sup>j</sup>* + Δ*σ<sup>j</sup>*
	- (d) Compute the internal force vector: **I**(*i*+1)
	- (e) Solve for the new accelerations: **u**¨ (*i*+1) = **M**−<sup>1</sup> **<sup>F</sup>**(*i*+1) <sup>−</sup> **<sup>I</sup>**(*i*+1)
	- (f) Compute the velocities at new mid-time: **u**˙ *<sup>i</sup>*<sup>+</sup> <sup>3</sup> <sup>2</sup> = **u**˙ *<sup>i</sup>*<sup>+</sup> <sup>1</sup> <sup>2</sup> + <sup>Δ</sup>*<sup>t</sup>* (*i*+1)+Δ*t* (*i*) <sup>2</sup> **<sup>u</sup>**¨ (*i*)

#### *2.5. Development of User Material Subroutine*

The active behavior of the muscle is not provided in the libraries of the commercial finite element codes. It is therefore necessary to implement the active behavior in the form

of a user-defined stress update algorithm. This was implemented in the finite element code ABAQUS/Standard by means of a UMAT. Additionally, for ABAQUS/Explicit, a VUMAT was also developed. Much of the coding involved in the two algorithms is the same but there are several key issues that must be addressed to maintain consistency of the results between the two solvers. These subroutines, written in Fortran, implemented the behavior of the material in the form of a stress update algorithm that is called at each integration point for every iteration during the finite element simulation. At these integration points, it is also necessary to define the anisotropy of the material to form the different tensors and to obtain the set of invariants. The subroutines were able to read a discretized fiber orientations from an external file during the first time increment.

The most important difference between the two programmed subroutines is that the explicit one does not update the tangential stiffness matrix. Nevertheless, when writing the implicit subroutine this matrix should be accurately represented to obtain correct solutions. The initial time increment used in ABAQUS/Standard is chosen by the user, and subsequent increments are controlled by an automatic incrementation control. To determine the size of the initial time increment in ABAQUS/Explicit a bogus set of tiny strain increments are passed to the VUMAT at the start of the analysis. From the stress response of the material, a conservative value for the stable time increment is calculated [16].

#### *2.6. Eyeball and EOM Finite Element Model*

The geometrical data of the model were obtained from the database of Mitsuhashi et al. [24], which was created from a whole-body set of 2 mm interval MRI images of a male volunteer. Some adjustments to the surfaces were made by hand to remove discontinuities and increase smoothness. The contours of the right eyeball and the four recti muscles were defined. As intorsion and extorsion movements were not considered in this study, the oblique muscles were not included in the model [10] To prevent the eye muscles from slipping away while the globe rotates, connective tissue surrounds the globe and stabilizes the muscles acting as pulleys and serving as functional EOM origins [25].

Figure 2 shows the finite element mesh used in the model. Solid hexahedral finite elements were used to mesh the eyeball and the EOMs, whereas one-dimensional truss elements simulated the action of the connective tissue pulleys. In Table 1, the number of nodes and elements of the solid parts of the model are shown together with their volume and mass according to the eyeball density [26] and muscle density [27].

**Table 1.** Mesh size, volume and mass of the model parts.


The set of mechanical properties used for the muscle material behavior is included in Table 2. Both passive and active parameters were adapted from a previous work [19]. The passive properties of connective tissues such as the tendon ends of the EOMS and the pulleys were taken from Calvo et al. [5]. For the latter, circular sections of 1 mm radius were applied to the truss elements and a rigid body constraint was applied to the eyeball.

**Figure 2.** Finite element mesh of the right eyeball, the EOMs and the zones of tendinous insertions obtained from Mitsuhashi et al. [24] incorporating connector elements to simulate the action of soft tissue pulleys.

**Table 2.** Parameters considered in the model for the material behavior of the skeletal-muscle tissue [19].


To account for the anisotropy present in the muscles due to the presence of fibers, a set of directions was generated inside the volume of the EOMs (Figure 3). These directions define both the passive behavior of the tissue and the direction in which the active contraction occurs.

Each of the ends of the EOMs in contact with the eyeball was tied, fixing the three degrees of freedom of both surfaces. The nodes of the other end of the muscles were clamped, restraining all movements. No contact was considered between the surfaces of the model since no interaction of the different parts of the model was detected in the movement. The truss elements connected nodes of the EOMs to a common node that was pinned. Finally, a rigid solid constraint was applied to the eyeball fixing an instant center of rotation at the center of the sphere.

**Figure 3.** Representation of the muscle and collagen fiber orientations.

#### **3. Results**

As the maximum time increment in a dynamic implicit algorithm depends on the typical period of vibration of the system, the natural frequencies and mode shapes were obtained for the model of the right eyeball, the pulleys and the EOMs. The initial four mode shapes are represented in Figure 4 which correspond to bending shapes of the four different muscles at frequencies from 55.33 to 84.16 Hz. The lower natural frequency corresponds to the inferior EOM muscle which is the muscle with the largest volume. Taking the largest characteristic frequency obtained, a maximum recommended increment for an implicit dynamic simulation should be less than 0.001 s.

The model was analyzed under four different movement scenarios characterized by the evolution of the activation function in Equation (11) along time. Figure 5a shows the activation function *f*2(*t*) considering four *s*<sup>1</sup> parameters that will induce eyeball movements from a very slow one *s*<sup>1</sup> = 10 to an intended nearly instantaneous *s*<sup>1</sup> = 100. Although all the muscles contribute in a certain way to the eyeball movements (elevation–depression and abduction–adduction), the results presented in this work are those obtained for the single activation of the medial EOM (adduction), as shown in Figure 5b. In this figure, the distribution of the maximum principal stress is presented. As can be observed, the maximum values are reached in the medial EOM at the final point of the simulation at *t* = 0.3 s.

The inertia of the system provides a response slower than the activation function, as shown in Figure 6a. The displacement in the *x* direction is represented for a point located at the center of the pupil area in the simulated eyeball. Implicit and explicit results are presented for the four activation function parameters. Comparing both algorithms, differences under 3% at the end of the simulations were found. Figure 6b shows the *x* displacement field at the end of the adduction movement considering *s*<sup>1</sup> = 10 for the implicit simulation.

To compare the performance of both methodologies, the evolution of the total strain energy (Figure 7a) and total kinetic energy (Figure 7b) was analyzed in the model. The total strain energy accumulated by the model for both algorithms and for all the activation signals is nearly the same at the beginning of the simulation. In contrast, at the end of the simulation, the models calculated with the explicit algorithm show an extra level of internal energy compared with the implicit simulations. As can be observed, the total energy is higher when the contraction velocity increases induced by the activation function. These higher levels of total energy are a consequence of a greater amount of kinetic energy when increasing the contraction velocity. The total kinetic energy is shown in Figure 7b and again the explicit results outperform the implicit ones.

**Figure 4.** Initial four mode shapes obtained for the right eyeball and EOMs system. These modes correspond with bending of the EOMs at natural frequencies of: (**a**) 55.33 Hz for the inferior EOM; (**b**) 72.57 Hz for the lateral EOM; (**c**) 72.84 Hz flexion mode for the superior EOM; and (**d**) 84.16 Hz for the medial EOM.

**Figure 5.** (**a**) Evolution of the activation function *f*2(*t*) with *α* = 1 for different *s*<sup>1</sup> parameter values simulating four contraction velocities; and (**b**) maximum principal stress distribution for the medial EOM activation with *s*<sup>1</sup> = 10.

**Figure 6.** (**a**) Evolution of the displacement of a node located at the center of the surface of the mesh where the pupil is located for all the activation signals considered using both algorithms; and (**b**) displacement field in the *x* direction for an implicit simulation with *s*<sup>1</sup> = 10 at *t* = 0.3 s.

**Figure 7.** (**a**) Evolution of the total strain energy for all the activation signals considered using both algorithms; and (**b**) evolution of the total kinetic energy for all the activation signals considered using the implicit algorithm.

Finally, a comparative analysis of the computational time is summarized in Table 3. As can be seen, the implicit algorithm takes only a 5% of the time spent with the explicit algorithm without mass scaling. A series of global mass scaling factors was applied to the model. Increasing this factor, higher values of fictitious mass are added, penalizing the range of motion of the system. Although a factor of 100 notably reduces the simulation time to levels near that of the implicit algorithm, the maximum displacement at the end of the simulation differs by 36.5% with respect to that obtained with no mass scaling.


**Table 3.** Computational time in percent relative to the dynamic explicit simulation without mass scaling (m.s.) for both algorithms and incorporating different m.s. factors.

#### **4. Discussion**

When studying extraocular mechanics, muscle activation and deformation are important parameters to characterize the movement of the eyeball [11]. In this paper, the activation of the medial EOM is analyzed under four activation signals which induce increasing contraction velocities. This function was simplified unlike more realistic previous models [8,9,19] to reduce the computation of unnecessary terms. This was assumed considering that all muscle fibers are recruited during activation (tetanic contraction). The range of motion associated with the region that represents the pupil in the model is in good agreement with those reported previously in the literature [4,10], although other authors simulated even larger angles of rotation [11]. The maximum horizontal or *x* displacement in Figure 6a can be translated according to the position of the instant center of rotation to a rotational angle of 7.4 degrees, which is far from the 20 degrees simulated in the work of Wei et al. [11]. Differences in the maximum isometric stress and in the force length relationship could be responsible for this disagreement since the properties of the muscle active behavior incorporated in this model were taken from a previous work [19]. The use of a single model with a particular geometry also limits the comparison with previous results but, on the other hand, it proves the potential of the methodology to develop a functional model based on medical image.

The results obtained for the four activation signal paths (Figure 6a) show that the predictions of the two algorithms differ at the end of the simulations. Larger time increments in the implicit method could lead to underestimating the contraction velocity (Equation (21)) implemented in the user subroutine, and consequently the muscle contracts more slowly.

Using the explicit analyses, the mass-scaling option is available to increase the stable time increment by artificially adding mass to the system. Although mass-scaling could decrease the mean computational time in the simulation of the eyeball movement (see Table 3), the range of motion is reduced drastically. It has been suggested in the literature that mass-scaling results are acceptable when the proportion of kinetic energy to strain energy is less than 5% [12]. In this case, the dynamic effect is negligible and problems can be solved with quasi-static solutions. As can be observed in Figure 7, this ratio is not satisfied in this model at the beginning of the simulations. Future analysis should consider increasing the maximum isometric stress developed by the activated muscle to explore whether it is possible to balance the addition of mass to obtain more realistic results.

As pointed out by other authors [12,16] and indicated by the results in Table 3, for simpler loading conditions, the implicit method takes a shorter solution time. In the case of loading conditions involving contact between the muscles and the eyeball or even incorporating the orbital fat, the explicit method will be the preferable choice [14,15]. Furthermore, the problem solved with this method can be easily parallelized in separated computer processors since the inverse of the lumped mass matrix can be split in decoupled set of equations.

#### **5. Conclusions**

In this paper, a comparison between implicit and explicit dynamic algorithms is presented and applied to model the 3D motion of the eyeball subjected to the action of EOMs. Our high-speed simulations showed that the dynamic implicit algorithm offers a substantial reduction in the required computational time in a model with no contact interactions between the surfaces. Although mass-scaling can provide a reduction in the computational time with the explicit algorithm, it is not recommended for high-speed movements taking into consideration the activation of the muscle tissue, due to the system increment of mass inertia.

**Author Contributions:** Conceptualization, methodology, and writing—review and editing, J.G. and B.C.; software, J.G.; and project administration, B.C. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Spanish Ministerio de Ciencia, Innovación y Universidades grant number DPI2017-84047-R and the Department of Industry and Innovation (Government of Aragon) through the research group Grant T24-20R (co-financed by Feder 2014-2020: Construyendo Europa desde Aragon). Part of the work was performed by the ICTS "NANBIOSIS" specifically by the High Performance Computing Unit (U27), of the CIBER in Bioengineering, Biomaterials & Nanomedicne (CIBER-BBN at the University of Zaragoza).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

