OpenEyeSim: A biomechanical model for simulation of closed-loop visual perception

Abstract

We introduce OpenEyeSim, a detailed three-dimensional biomechanical model of the human extraocular eye muscles including a visualization of a virtual environment. The main purpose of OpenEyeSim is to serve as a platform for developing models of the joint learning of visual representations and eye-movement control in the perception–action cycle. The architecture and dynamic muscle properties are based on measurements of the human oculomotor system. We show that our model can reproduce different types of eye movements. Additionally, our model is able to calculate metabolic costs of eye movements. It is also able to simulate different eye disorders, such as different forms of strabismus. We propose OpenEyeSim as a platform for studying many of the complexities of oculomotor control and learning during normal and abnormal visual development.

Introduction

Vision during natural tasks is an active process involving different kinds of eye movements. These eye movements are produced through the activation of six extraocular eye muscles (EOMs) for each eye. Despite the importance of eye movements for natural vision, there is a surprising shortage of closed-loop computational models that describe how the activation of EOMs moves the eyes, how this changes visual input, how this in turn drives the next eye movement, and so forth. Existing models typically focus exclusively on either the motor-control part or the visual-processing part. Examples of the former are models trying to explain how the characteristic velocity profile of saccades comes about. Examples of the latter are models that describe how the brain may estimate image motion or binocular disparity, how it may recognize objects in images, and how neural responses are affected by attention mechanisms. The lack of closed-loop models is problematic for several reasons. First, neural responses in visual cortex are thought to be modulated by feedback from the eye muscles (Trotter, Celebrini, Stricanne, Thorpe, & Imbert, 1992). Second, the normal development of visual perception during infancy critically depends on the proper self-calibration of visual sensorimotor loops for stereoscopic vision and motion vision. Third, developmental disorders of vision such as strabismus and amblyopia, which affect several percent of the general population but whose origins are poorly understood (Quoc & Milleret, 2014), may be best viewed as failures of this self-calibration of sensorimotor loops. 

What has held back the development of closed-loop models addressing these issues is the lack of a proper modeling platform and software environment to facilitate their development. Here we offer a solution to this problem: We present OpenEyeSim, an open software environment for developing closed-loop models of visual perception. OpenEyeSim combines, for the first time, a sophisticated model of the biomechanics of the oculomotor plant with the ability to render a virtual environment and calculate how movements of the eyes change the visual input. In the following, we briefly review previous models of the biomechanics of the oculomotor plant that have led to the development of OpenEyeSim. 

Related work on biomechanical modeling

Early one-dimensional models of oculomotor mechanics (Clark & Stark, 1974; Robinson, 1981) were based on experimental measurements of EOMs, such as their tension–length relationship (Robinson, 1964). These models focused only on horizontal eye movements, but they laid the groundwork for subsequent studies. The first model that included mechanics in three dimensions was developed by Robinson (1975). This model, based on measurements of contractile muscle force, elastic muscle force, muscle lengths, and orbital geometry, was limited to static mechanics. It attempted to represent realistic muscle paths and empirical EOM innervation–tension–length relationships. It was suitable for studies on binocular alignment in normal and abnormal cases such as strabismus, as well as surgery planning. Later models—Simonsz's (Simonsz & Spekreijse, 1996) and the SQUINT model (Miller & Robinson, 1984)—included additional improvements. After further improvements, the last model was converted to the software Orbit 1.8 (Miller, Pavlovski, & Shaemeva, 1999), which is suitable for studying different abnormalities and planning surgery. 

Most of the 3-D models already mentioned (Miller & Robinson, 1984; Robinson, 1975; Simonsz & Spekreijse, 1996) can simulate only static eye positions. To study the neural control of eye movements, models using simplified properties of EOMs have been developed (Quaia & Optican, 1998; Raphan, 1998). These models do not take into account the anatomical variations of different EOMs, such as muscle length, cross-sectional areas, and muscle forces. In some previously developed models, muscle forces were assumed to be proportional to the innervation, whereas real muscle forces are complex functions that depend on muscle length, velocity, and innervation. To our knowledge, only the model presented by Wei, Sueda, and Pai (2010) incorporates proper muscle dynamics and produces realistic eye movements. However, all these models only cover the biomechanical component, neglecting any visual processing. 

We present a novel 3-D biomechanical model of EOMs with embedded visualization of a virtual environment. Our model uses the biomechanical simulator OpenSim (Delp et al., 2007) and extends its abilities in several ways. It incorporates several important features: 

•  
  Realistic EOM force dynamics, which to our knowledge have only been presented in one previous 3-D model of eye muscles (Wei et al., 2010)
•  
  Muscle pulleys—extraocular connective tissues which stabilize muscle paths
•  
  Realistic muscle paths, based on measurements in humans
•  
  Visualization of a virtual environment, which is crucial for simulations of closed-loop visuomotor control and learning

In the following sections we describe the architecture of our model, as well as the simulator. Then we describe different components of our model and show their correspondence to measurements on humans. Finally, we show that the model can generate different kinds of eye movements and demonstrate how it can be used to simulate different forms of strabismus. 

Methods

OpenSim

Our model is based on the biomechanical simulator OpenSim (Delp et al., 2007). This is an open-source software that already implements a variety of muscle models and includes many freely available musculoskeletal models. It is used for studies on musculoskeletal dynamics and for modeling different disorders. One of its components, which is the most important for our studies, is a model of the forward dynamics. This allows us to use muscle activations as input, from which the simulator calculates all arising forces in the model as well as their interaction and provides as output a simulated movement. 

We have developed an extension to OpenSim to allow the simulation of closed-loop visually guided behaviors. It enables us to get pictures from virtual cameras placed in the eyes. We have written additional classes, which allow for fast off-screen rendering. This enables us to perform our simulations in parallel. With the visualizer embedded in OpenSim, this is not possible, because of the architecture of internal classes. 

Architecture of the model

OpenEyeSim includes six EOMs per eye: the four rectus and the two oblique muscles. These six EOMs are controlled by cranial nerves and generate the forces which rotate the globe during different kinds of eye movements. The four rectus muscles are straight muscles: medial, lateral, inferior, and superior rectus (Figure 1). The medial and lateral rectus muscles form an agonist/antagonist pair and produce eye movements in the horizontal plane. The agonist/antagonist pair of superior and inferior rectus muscles are responsible for producing eye movements in the vertical direction, but due to the location of the muscle insertion points they also affect the torsional component—i.e., rotations of the eyeball around the line of sight. 

Figure 1

 

Architecture of the model.

The four straight muscles and the superior oblique muscle originate at the annulus of Zinn, a tendon ring lying at the back of the orbit. The superior oblique muscle passes through the cartilaginous trochlea, which alters its path by 54°. The inferior oblique originates near the nasal bone in the orbital wall anteroinferior and connects posterior to the eyeball. The specific locations of the oblique muscles enable them to perform torsional eye movements as their main action. It has been suggested that EOMs are highly coupled to the orbital wall through special connective tissues (Miller, 1989). Further studies have shown strong evidence for the presence of such structures, as well as their significant role in the kinematics and neural control of the orbital plane (Demer, Oh, Clark, & Poukens, 2003; Miller, 2007). These structures, called muscle pulleys, were fully described by Miller (2007). He also suggested that the well-known Listing's law describing torsional eye movements is implemented in the muscle-pulley structures and not via brain-level control. The idea of muscle pulleys has also been supported by a few experimental studies. Ghasia and Angelaki (2005) showed that cyclovertical motoneurons do not modulate their firing during an eccentric pursuit, as would be necessary if the brain stem implemented Listing's law. Klier, Meng, and Angelaki (2006) stimulated the abducens nerve and nucleus, downstream of all neural circuits that might contribute to the implementation of Listing's law, and found that eye movements nevertheless had Listing's kinematics, proving that ocular-plant mechanics are capable of implementing Listing's law without neural assistance. 

In our model we use a passive pulley model, which helps to avoid sliding of the muscle over the eyeball and constrains the transverse shifts of muscle paths. We have chosen the passive pulley model in order to keep the model simple and provide faster simulation speed with a possible reduction of model accuracy in tertiary gaze positions (Miller, 2007). Table 1 shows the positions of such muscle pulleys, as well as muscle origins and attachment points to the globe. These data are based on measurements on cadavers and have been revised by Miller and Robinson in their model (1984). A detailed description of the implemented geometry is given in the supplementary material (in the form of an OpenSim simulation model file .osim). 

Table 1

 

Location of point of origin and other parameters for the right eye (all in mm).

 

Coordinate system

The eye can rotate in three dimensions around fixed axes in the head-constrained coordinate system. Movement of the eye nasally is called adduction, and temporally, abduction. Elevation and depression correspond to eye movement up and down in the vertical direction. Torsional eye movements rotate the eye around the visual axis. A nasalward rotation is called intorsion and a temporalward rotation is called extorsion. The coordinate system for the right eye is shown in Figure 2. We describe eye rotations in the Fick coordinate system, which describes an eye position as the result of rotating the eyeball from its primary position via three successive rotations around different axes: an initial rotation around a head-fixed vertical axis (Z), a subsequent rotation around an eye-fixed horizontal axis (X), and a final rotation around the line of sight (Y). 

Figure 2

 

Rotational directions for the right eye in the Fick coordinate system. The vertical (Z) axis is head fixed, while the horizontal (X) axis is eye fixed.

Muscle model

In order to simulate musculoskeletal dynamics, we use the modern implementation of the Hill-type muscle model (Hill, 1938), which is widely applied in describing the force-generation mechanism. It contains three elements: an active contractile element, a passive elastic element, and a tendon element. 

We use the Millard model (based on the Hill-type model) implemented in OpenSim (Millard, Uchida, Seth, & Delp, 2013), which allows us to manually fit force-generation dynamics. EOMs have two separate layers—the orbital and the global. For simplicity, we do not model this separation into two layers. 

Muscle dynamics

Based on a traditional description of the Hill-type model, the active force dynamics rely on the active force–length relationship and the force–velocity relationship:  

 

where a ∈ [0, 1] is the muscle activation (0 = fully relaxed muscle, 1 = fully innervated muscle), f_FL is the force–length relationship, f_FV is the force–velocity relationship, l is the muscle length, and v is the contraction speed. The passive elastic force F_PE is calculated as a function of muscle length:  

 

where l is the muscle length and f_p is the passive force–length relationship (for a detailed explanation, see Supplementary Appendix A).  

These relationships for active and passive forces were measured experimentally (Robinson, O'meara, Scott, & Collins, 1969). Models of force–velocity have been proposed based on Hill's skeletal-muscle model, experimental data on rat EOMs, and maximum saccadic velocity of human eyes (Robinson, 1981). We adopt the parameters of our muscle model to fit the active force–length behavior of a fully innervated EOM, as well as passive force–length behavior of the SQUINT model (Miller & Robinson, 1984), which to our knowledge best approximates the experimental data. Fitting was done by using the fminsearch solver in MATLAB, which is based on the Nelder–Mead simplex method (Lagarias, Reeds, Wright, & Wright, 1998). Figure 3 shows the resulting fit. 

Figure 3

 

(A) Force–length relationship. The x-axis shows the muscle stretch L relative to the optimal fiber length L0. The y-axis shows the muscle force F relative to the force in the primary position F0. Active force corresponds to F_CE, passive force corresponds to F_PE. (B) Force–velocity relationship used in our model.

The tendon element in the three-element Hill model has a very high stiffness in our EOM model. Robinson (1981) has pointed out that the force mechanism of the EOM tendon element has not been quantified well due to the lack of physiological data. 

Dynamic properties of muscles

The behavior of the selected Hill-type muscle model depends mainly on a few parameters: the maximum isometric force, the maximum contraction velocity, the tendon slack length, and the optimal fiber length. In order to estimate the maximum isometric force, we took data of the medial rectus muscle from physiological experiments (Collins, Carlson, Scott, & Jampolsky, 1981) and estimated other forces by simply considering the ratio of cross-sectional areas of the muscles (Miller & Robinson, 1984; see Table 2). These cross-sectional area data were taken from measurements of dissected slices of eye muscles by Volkmann and revised by Robinson (1975). Tendon slack lengths and optimal fiber lengths are shown in Table 1. Due to different histological structures of EOMs compared to skeletal muscles, EOMs have a higher fraction of fast-twitch fibers (Wasicky, Ziya-Ghazvini, Blumer, Lukas, & Mayr, 2000) and a different force–velocity behavior. Also, these fibers have slightly different behavior compared to fibers in skeletal muscles. Our parameters, such as maximum contraction velocity and additional parameters explaining changes in force–velocity slopes, are based on the measurements of EOMs (Collins, 1971). This provides our model with the ability to generate eye rotations in a range up to 900°/s–1000°/s (Siegelbaum & Hudspeth, 2000), which is consistent with the peak velocity of saccadic eye movements. Also, because of a different structure of the neural control of eye movements, activation and deactivation delays are lower than in skeletal muscles, at about 5 ms (Robinson, 1968). 

Table 2

 

Maximum isometric force (g) and cross-sectional area relative to the lateral rectus muscle.

 

Other components of EOM dynamics

Additional components of our model are forces generated by passive orbital tissues. These tissues include all nonmuscular suspensory tissues, such as Tenon's capsule, the optic nerve, the fat pad, and the conjunctiva. The force–displacement curve of the net elasticity P has been initially measured by Robinson and can be represented as  

 

where e is the degree of rotation and k_p and k_c are coefficients. Robinson's measurements (Robinson et al., 1969) found k_p = 0.48 g/° and k_c = 1.56 g/°³. However, later studies by Collins (Robinson, 1981) showed that the proper value for k_p is 0.32 g/°. These forces try to keep an eyeball in its primary position and serve its stabilization. In our model we used the built-in functionality of OpenSim; we were therefore able to include only the first term in our simulations. However, neglecting the cubic component should not result in a big change in the model dynamics. For example, at 20° eye rotation the contribution of the cubic component is around 11% of total elasticity force.  

Metabolic costs

Using the OpenSim extension based on the work of Umberger et al. (2003), we are able to calculate the metabolic costs of eye movements. The rate of metabolic energy consumption is calculated as  

 

where A_i is the activation heat rate of muscle i, M_i is the maintenance heat rate, S_i is the shortening heat rate, and W_i is the mechanical heat rate (see Supplementary Appendix B for details). A calculation of these values requires additional parameters—muscle masses and ratios of slow-to fast-twitch fibers. Values for the ratio of slow- to fast-twitch fibers were roughly estimated based on histological studies on EOMs (Wasicky et al., 2000) and are 15.1%/84.9% for central muscle parts and 14.3%/85.7% for peripheral muscle parts. These studies showed that nearly 85% of muscle fibers found in EOMs have the same characteristics as fiber types found in skeletal muscles. Another 15% of fibers have different specific structures, which are not present in skeletal muscles. Because of the lack of information on these fiber types, we neglect their contribution (which has been shown to be present in some saccades; Porter, Baker, Ragusa, & Brueckner, 1995) and simply calculate metabolic costs based on the ratio of slow- to fast-twitch fibers of 85% of total volume. Figure 4 shows an example of such calculations of metabolic costs. Different lateral and medial rectus muscle activations that lead to the same static eye positions are indicated by white lines. The energy consumption for using such muscle activations for 1 s is presented by different colors. This figure reveals the redundancy in oculomotor control: Different muscle activations lead to the same eye position.  

Figure 4

 

Example calculations of metabolic costs. The x- and y-axis show muscle activations (0 = relaxed, 1 = fully innervated). The solutions of the Covariance Matrix Adaptation Evolution Strategy for −30, 0, and +30 horizontal rotations are marked with magenta.

Muscle paths

Wrapping objects are another important component of our model, which help us to simulate the proper dynamics. Using wrapping objects implemented in OpenSim allows us to model realistic force directions as the muscle runs over the eyeball. From the point of origin to the point of tangency, a force direction arising from simulation is aligned along one line, but from the point of tangency to the insertion point, the direction of forces is distributed as shown in Figure 5. 

Figure 5

 

Force direction over the wrapping object.

Calculation of muscle activations

In order to calculate activation levels of the eye muscles at different static positions we used a black-box optimization algorithm, the Covariance Matrix Adaptation Evolution Strategy (Hansen, 2006). This is a stochastic optimization algorithm which searches for local optima by sampling and evaluating points p͐_i in the search space. Such a point p͐_i is sampled from a multivariate Gaussian distribution 𝒩(p͐(t); Σ(t)), where p͐ is the mean, Σ is the covariance, and t is the iteration number. The sampling space is modified over iterations based on successful samples. The algorithm finds a solution whenever a sample reaches a predefined acceptance threshold. 

The search is initialized in the primary eye position, which corresponds to near-zero rotation values. From the search distribution, n samples are drawn:  

 

where each sample p͐_i is used as an input to the forward-dynamics module of the OpenSim simulator. The muscle activations are held constant until the eyeball has settled into the corresponding gaze position. The produced rotations, as well as muscle activations, are used to calculate a reward:  

 

where Display Formula , Display Formula , and Display Formula are rotation angles (in degrees) of the eyeball, which correspond to sampled solutions p͐_i; Display Formula is a corresponding activation of muscle j for sampled solution i; and s is a scaling factor. This reward function penalizes high activation levels while favoring eye rotations close to the desired. We have chosen a value of s = 0.5. As shown in Figure 4, solutions based on this reward function correspond to low metabolic costs. The samples p͐_i, their corresponding rewards q_i, and the search parameters p͐(t) and Σ(t) are used to generate new search parameters p͐(t + 1) and Σ(t + 1) to produce the next generation of samples.  

The search stops when the reward reaches the acceptance threshold. If the search algorithm converges to a local optimum above the acceptance threshold, the search distribution is reset and the search is restarted. All parameters of the search algorithm are provided in Supplementary Appendix C. 

Rendering component

Our model contains a graphical component, which allows us to perform closed-loop experiments. We are able to model different virtual environments and receive information from two cameras placed in the centers of the eyeballs. The interpupillary distance is set to 62 mm. The rendering speed using our off-screen rendering extension reaches 40 ms per iteration. For rendering we use an OpenGL engine. The main rendering class lies on the side of the OpenSim simulator, whereas in the standard version it is part of the library that describes component classes. The connection to OpenSim from this library is realized via a pipeline interaction between processes. From our experience, such connections become unstable for large numbers of rendering requests. In comparison, our rendering extension is very stable. Another important component which was added to the rendering classes is the Simple OpenGL Image Library, which enables us to apply textures to the objects. OpenSim developers provide such ability only in a Windows external application, which is only suitable for creating a precise biomechanical model, rather than closed-loop simulations. An example of a virtual environment with the rendered images from two eye cameras is shown in Figure 6. 

Figure 6

 

Example of a virtual scene and rendered images for the left and right eyes. Red crosses mark the fixation point.

Results

We begin by demonstrating how OpenEyeSim can be used to model different forms of strabismus. Then we show how OpenEyeSim can be used to model the time course of saccadic eye movements. 

Hess–Lancaster test

The Hess–Lancaster test is a widely used clinical test for assessing binocular alignment in the nine diagnostic positions of gaze. It may be used to quantify comitant and incomitant strabismus in children and adults. There are many similar testing methods (e.g., Harm-Screen, Lee-Screen) which use the same testing principle based on displaying different images to the two eyes. 

The main procedure of the Hess–Lancaster test is divided into three steps: 

1.  
   The test subject wears glasses with green–red filters with the red filter in front of the right eye, which is called the fixing eye. Then the patient is given a green-light pointer while the examiner uses a red-light pointer.
2.  
   The examiner places the red-light dot into different positions on the Hess screen covering the main working area of gaze positions. The tested subject is asked to bring the green-light dot (seen by the following eye) over the red-light dot. In the normal case, these dots should overlap.
3.  
   The same procedure is repeated with changed filters.

If a healthy subject fixates, the fixation can be reached with normal innervation. If, for example, a subject has a palsy of the right medial rectus muscle, the green pointer will point to a position that is different from the red-dot position. 

Figure 7A shows the Hess–Lancaster diagram. The blue dots are the gaze positions that the subject is asked to fixate with the fixing eye. The red points are the gaze positions that the subject fixates with the following eye. The difference between the positions of the blue and red points represents the deviation of binocular alignment. In this example we model an isolated medial rectus palsy. In this disorder, only the medial rectus muscle is affected. 

Figure 7

 

(A) Simulation of isolated medial rectus palsy. (B) Simulation of abducens nerve palsy. (C) Simulation of superior rectus palsy. (D) Simulation of Duane's syndrome. The data represent single simulations for a given change of parameters.

In order to simulate such disorders, we follow this procedure: 

1.  
   Calculate activation levels for eye muscles under normal conditions for different points on the Hess–Lancaster diagram. Such calculations are performed using the black-box optimization algorithm presented in Methods. The obtained muscle activation values for the nine points on the Hess–Lancaster diagram are given in Supplementary Appendix C.
2.  
   Scale the maximum activation levels of affected muscles to reduce calculated forces at different gaze positions.
3.  
   Produce forward-dynamics simulations with reduced activation signals to calculate the altered eye positions.

In all cases, the pattern of altered gaze positions resembles those seen in patients (Kaufmann & Steffen, 2012) and in the models of Haslwanter et al. (2005) and Wei et al. (2010). 

Isolated medial rectus palsy

An isolated medial rectus palsy is a case of an ocular motility disorder—internuclear ophthalmoplegia. The disorder happens because of injury or dysfunction in the medial longitudinal fasciculus, a tract which connects the paramedian pontine reticular formation of the contralateral side to the oculomotor nucleus of the ipsilateral side. 

Patients with this disorder have a problem with adduction, which leads to horizontal diplopia. Figure 7A shows an example of a simulated isolated medial rectus palsy. In order to simulate this, we reduced the maximum activation of the medial rectus muscle by 40%. 

Abducens nerve palsy

Abducens nerve palsy is an ocular motility disorder associated with the sixth cranial nerve. It causes improper contraction of the lateral rectus muscle, which abducts the eye. This palsy is the most common among the isolated motor-nerve palsies. 

Usually, individuals with this kind of palsy have horizontal diplopia and an esotropia in primary gaze. The deviation of binocular alignment is greater in secondary and tertiary positions. Typically, a check for sixth nerve palsy involves excluding palsies of other nerves. In sixth nerve palsy only deviations in the horizontal plane are observed; deviations in vertical and torsional directions are not observed. 

Abducens palsy is frequently caused by stroke, trauma, or postviral syndrome. Figure 7B shows an example of simulated abducens nerve palsy. In order to simulate this disorder we reduced the maximum activation of the lateral rectus muscle by 50%. 

Isolated palsy of the superior rectus muscle

Isolated palsy of the superior rectus is a rare ocular motility disorder associated with the third cranial nerve. The superior rectus muscle produces an elevation in a primary position. Because of the palsy, an overacting of the antagonist inferior rectus muscle moves the eye downwards. The elevation is better in adduction than in abduction. Due to the overacting of the inferior rectus and inferior oblique muscles, there is noticeable extorsion. 

Figure 7C shows an example of a simulated superior rectus muscle palsy. Simulations have been performed by reducing the maximum activation level of the superior rectus muscle by 25%. 

Duane's syndrome

Duane's syndrome is an ocular motility disorder characterized by inability of the eye to move outward (abduction) and limited ability to move inward (adduction). Duane's syndrome happens because of miswiring of the brain stem due to failures during development. The sixth cranial nerve is affected and does not develop properly. There is also unusual innervation of the third cranial nerve, which normally controls the medial rectus muscle. Additionally, upward or downward deviation may occur with attempted adduction due to a leash effect. Duane's syndrome is classified into three types: 

•  
  Limited abduction with or without esotropia
•  
  Limited adduction with or without exotropia
•  
  Limitation of both abduction and adduction and any form of horizontal strabismus

Figure 7D shows an example of the third type of Duane's syndrome. In order to simulate this disorder, we reduced the maximum activation levels of lateral rectus and medial rectus muscles by 40%. 

Response to ramp- and step-function inputs and noise

In order to illustrate the dynamics of our model, we modeled two different inputs using step and ramp functions. Inputs contain control signals for the following sequence of positions: [0; 0]°, [12; 6]°, [−18; −6]°, [24; 0]°, [0; 0]°. The resulting rotations, as well as muscle activations, are shown in Figure 8. We also modeled a slow ramp-like activation pattern of such positions: [0; 0]°, [12; −12]°, [−18, −6]°, [24, −6]°, [0; 0]°. These patterns could be used in a model of smooth-pursuit eye movements. The resulting rotations, as well as muscle activations, are shown in Figure 9A. We also simulated an ocular microtremor by injecting noise into muscle activations. The resulting rotations and activations are shown in Figure 9B. Corresponding movies are contained in the supplementary materials. 

Figure 8

 

(A) Simulated muscle activations for a step input signal (top) and corresponding rotation of an eyeball (bottom). (B) Simulated muscle activations for a ramp input signal (top) and corresponding rotation of an eyeball (bottom).

Figure 9

 

(A) Simulated muscle activations for a slow ramp input signal (top) and corresponding rotation of an eyeball (bottom). (B) Simulated muscle activations with injected noise (top) and corresponding rotation of an eyeball (bottom).

Saccades

We also tested our model's ability to generate realistic saccadic eye movements in response to biologically plausible muscle-activation signals. Since measurements of neural saccade commands in humans are not available, the neural signals of saccadic movements were taken from the experimental results obtained by Cullen on rhesus monkeys (Sylvestre & Cullen, 1999). 

Two different saccades were taken from these data: +13° and +22°, where the initial positions were near zero and the directions of the saccades were temporal. The corresponding average firing rates of oculomotor neurons and activations of the lateral rectus muscle are shown in Figure 10A and C (top). Muscle activations were obtained by normalizing the firing rates to [0 1], where zero activation means a totally relaxed muscle and one means fully innervated. Such normalization was done by applying an offset of 200 Hz, which corresponds to the firing rate at zero position, and scaling by a factor of 1/400, which corresponds to the difference between maximum firing rates from the data and firing rates corresponding to the zero position:  

 

 

Figure 10

 

(A) Firing rates for a saccade of 13° (top) and estimated activation of a lateral muscle (bottom). (B) Rotational components for the saccade of 13°. (C) Firing rates for a saccade of 22° (top) and estimated activation of a lateral muscle (bottom). (D) Rotational components for the saccade of 22°.

The results of this normalization are presented in Figure 10A and C (bottom). These control signals contain two main phases: the saccadic one, where this control signal increases sharply, and the fixation one, where the signal stays roughly constant with some tiny oscillations. Another component of force generation, the active-state tension, is also larger in monkeys; however, our normalization compensates for the difference between the properties of human and rhesus-monkey EOMs. Due to the lack of data, and for simplicity, the activations of other muscles have been fixed at a very small value of 0.05. These muscle activations do not produce strong forces, and with the ocular-tissue forces keep the eyeball at its zero position. Because we model horizontal saccades, an activation of superior and inferior eye muscles as well as oblique ones does not change the observed horizontal rotation. Because our saccades start at the zero position, activation of the medial rectus muscle does not exert any strong influence on the observed rotation. The resulting changes in the eye positions are presented in Figure 10B and D. The simulated trajectory is very similar to the trajectory measured in monkeys. The rotational components around the X and Y axes stay around zero. This confirms that the biomechanical part is modeled properly and the muscle paths configured correctly. 

A more thorough evaluation of the model would require simultaneous recordings of all six EOMs together with their neural control signals in humans. However, such data are currently unavailable. 

Discussion

We have introduced OpenEyeSim, a sophisticated biomechanical model of the six EOMs that incorporates a virtual environment for simulating closed-loop visual processing. To the best of our knowledge, OpenEyeSim is the first model of its kind. The goal of this article was to introduce OpenEyeSim as an open tool for vision scientists who are interested in modeling closed-loop visual information processing. The main features of OpenEyeSim are realistic EOM dynamics, muscle paths, and muscle pulleys, as well as the ability to render complex virtual environments providing visual stimulation. 

As we have demonstrated, OpenEyeSim can be used to investigate muscle activation patterns during static fixations and the control and dynamics of various eye movements. Our main reason for creating OpenEyeSim was to study the development of closed-loop visual perception. In particular, we are interested in active binocular vision and active motion vision and how the brain manages to calibrate such sensorimotor loops—or fails to do so during developmental disorders such as strabismus and amblyopia. This will be the topic of future work. While we made an effort to model the biomechanics of the oculomotor plant at a high level of realism, the model has a number of limitations. First, the maximum isometric forces of eye muscles were based on actual data for only the medial rectus muscle; values for other muscles were estimated based on their relative cross-sectional areas compared to the medial rectus muscle. Second, we used a passive pulley model and ignored the separation of EOMs into global and orbital layers. These and other improvements should be incorporated once more detailed measurements on the biomechanics of eye muscles become available. Third, for the forces generated by passive orbital tissues we considered only a simple linear model. For gaze locations far from the primary gaze position this is likely to become inaccurate. For studies of visual development, it would also be interesting to allow for online changes of kinematic and dynamic properties, such as the size of the eyeball, muscle properties, and interocular distance. While we have focused on the human oculomotor plant, it would also be very interesting to develop extensions for other organisms, in particular those which are commonly used in neuroscientific studies (mice, rats, cats, and primates). 

The visual-rendering component is based on a simple pinhole-camera model where the optical center is also the center of rotation. Much more realistic models of the imaging process in the human eye (and the eyes of other species) have been proposed and could be incorporated in the future (Artal, 2014; Emsley, 1948; Thibos, Ye, Zhang, & Bradley, 1992). From a developmental perspective it would also be interesting to be able to change the imaging properties of the model online. This could allow description of, e.g., the changing contrast sensitivity of infants during the period when sensorimotor loops for stereoscopic vision and motion vision are established and calibrated (Atkinson, Braddick, & Moar, 1977). 

In conclusion, with OpenEyeSim we hope we have created a versatile tool for vision scientists interested in modeling the manifold aspects of closed-loop visual perception. We invite our colleagues to use this tool and to help us improve and extend it in order to advance our understanding of the active nature of visual perception (https://simtk.org/projects/openeyesim). 

Acknowledgments

This work was supported in part by the German Federal Ministry of Education and Research under Grants 01GQ1414 and 01EW1603A, the Hong Kong Research Grants Council under Grant 618512, and the Quandt Foundation. 

Commercial relationships: none. 

Corresponding author: Alexander Priamikov.