Modelling 3D saccade generation by feedforward optimal control

An interesting problem for the human saccadic eye-movement system is how to deal with the degrees-of-freedom problem: the six extra-ocular muscles provide three rotational degrees of freedom, while only two are needed to point gaze at any direction. Measurements show that 3D eye orientations during head-fixed saccades in far-viewing conditions lie in Listing’s plane (LP), in which the eye’s cyclotorsion is zero (Listing’s law, LL). Moreover, while saccades are executed as single-axis rotations around a stable eye-angular velocity axis, they follow straight trajectories in LP. Another distinctive saccade property is their nonlinear main-sequence dynamics: the affine relationship between saccade size and movement duration, and the saturation of peak velocity with amplitude. To explain all these properties, we developed a computational model, based on a simplified and upscaled robotic prototype of an eye with 3 degrees of freedom, driven by three independent motor commands, coupled to three antagonistic elastic muscle pairs. As the robotic prototype was not intended to faithfully mimic the detailed biomechanics of the human eye, we did not impose specific prior mechanical constraints on the ocular plant that could, by themselves, generate Listing’s law and the main-sequence. Instead, our goal was to study how these properties can emerge from the application of optimal control principles to simplified eye models. We performed a numerical linearization of the nonlinear system dynamics around the origin using system identification techniques, and developed open-loop controllers for 3D saccade generation. Applying optimal control to the simulated model, could reproduce both Listing’s law and and the main-sequence. We verified the contribution of different terms in the cost optimization functional to realistic 3D saccade behavior, and identified four essential terms: total energy expenditure by the motors, movement duration, gaze accuracy, and the total static force exerted by the muscles during fixation. Our findings suggest that Listing’s law, as well as the saccade dynamics and their trajectories, may all emerge from the same common mechanism that aims to optimize speed-accuracy trade-off for saccades, while minimizing the total muscle force during eccentric fixation.


Introduction
This paper presents a computational model for the control of gaze in a simplified and scaled 3D mechanical eye model with tendon-driven actuation that approximates the extra-ocular muscle geometry of the human eye (see Figs 1 and 2). We propose an optimal open-loop control scheme that, when appropriately configured with cost terms that consider the minimization of accuracy, energy, saccade duration, and force in the muscles, generates the main kinematic and dynamics properties of saccadic eye movements found in primates. Previous studies either have focused on simpler 1D models (vertical/horizontal saccades) or analysed the kinematics and dynamic of saccade behavior independently. Our proposal demonstrates the explanatory power of unifying optimal control principles in both the kinematics and dynamics of a mechanical eye model that executes 3D saccade trajectories.

Saccade control
In order to generate a goal-directed movement, the brain has to resort to sophisticated neural control strategies that take into account the excess number of degrees of freedom, the noncommutativity of rotational kinematics, and the complex dynamics of the end effectors. The human oculomotor system is a well-studied example. The human eye has six extra-ocular muscles that can each pull the eye in a different direction (Fig 1). Because the eye is fully encased within its bony orbit, it cannot translate [1], thereby reducing the six degrees of freedom for rigid-body motion to three rotational degrees of freedom. However, to direct the gaze-line at any direction, only two coordinates need to be specified, e.g., the azimuth (θ) and elevation (φ) angles, leaving the ocular rotation around the line of sight (cyclo-torsion, ψ) unspecified. If not controlled, however, this poses a potential problem for eye movements, as mathematics holds that sequences of rigid-body rotations in general do not commute. Thus, the order in which the rotations are executed determines the final body orientation. In principle, this could lead to an unwanted accumulation of ocular torsion and significant localization errors [2,3]. It was already recognized by Donders and Helmholtz that the eye restricts its torsional roll angle in a unique way, depending only on the gaze direction, regardless the path taken by the eye to spawned numerous studies of gaze behaviour during human-robot interaction (see [9], for a review; [10]), not only for letting robots understand human gaze signals, but also (and more important for the present paper) to allow robots to display legible gaze behavior for a human interlocutor. One line of research tries to generate human-like behaviors in robots by recording human motions in the execution of particular tasks [11][12][13]. Another line of research tries to replicate the observed properties of human gaze behavior by directly programming the theoretical control laws into the robots [14,15].
While learning the gaze kinematics from humans, or directly programming the theoretical neuroscientific findings in robots are both viable approaches for designing working systems, they are limited in two aspects. First, the underlying principles that rule their emergence are left uncovered. Second, existing models address only very simplified settings.
The present paper considers the generation of saccades with fully unconstrained kinematics and develops a methodology to generate control systems that both satisfy mechanical optimality criteria (energy and stress), and performance criteria (duration and accuracy), while mimicking the kinematics and transient properties observed in natural gaze behavior. To validate our approach we have built a scaled up robotic prototype driven by three motors attached to elastic tendons, schematically shown in Fig 2. The tendons are attached to an artificial eye ball in a way that, although not matching its biological counterpart, shares a similar geometry. The other end of the tendons are connected in pairs to symmetric points on rods that are fixed to each motor. Unlike the biological muscles, these elastic tendons can only apply force when stretched so they are actuated in pairs by the rotation of the motor rods. The tendons for vertical and torsional muscles pass through fixed routing points to adjust their pulling directions to provide enough ocular range. We built a simulation model of the system in order to study the role of optimal control methods in the emergence of relevant oculomotor behaviors.
The proposed approach may lead to more energy efficient and more durable robots, with more flexibility in replicating the complex repertoire of oculomotor behaviors exhibited by humans. Building humanoid-representative oculomotor robots may also be useful for groundtruth benchmarking of many eye tracking systems.

3D kinematics
As will be described in more detail below, a 3D eye orientation is most efficiently described by a single-axis rotation of the eye from the primary position. Thus, eye orientations can be parametrized by the direction of the rotation axis,n, knk ¼ 1, and the amount of rotation, ρ, around that axis needed to attain the orientation: f ðn; rÞ ¼ rn ¼ e ¼ ðe x ; e y ; e z Þ, the so-called Euler vector, with (x, y, z) a head-fixed, right-handed Cartesian reference frame (x = torsion, y = vertical, z = horizontal rotation). Due to a closer relationship to quaternions, it is convenient to adopt the equivalent half-radian rotation vector representation [16]: r ¼ ðr x ; r y ; r z Þ ¼ tan ðr=2Þn.
With the head upright and not moving, and gaze directed at infinity, Donders' law is further constrained by stating that the rotation axes describing all eye orientations are restricted to a plane (Listing's law (LL); [4][5][6]17]). Fig 3 shows an example of more than 3500 rotation axes, computed during fixations and saccades made by a monkey in the lab (head-restrained, spontaneous eye movements in the light). It can be appreciated from this figure that LL is obeyed with remarkable precision, as the standard deviation of the torsional component of monkey eye orientations is typically less than one deg (torsional range in [-2.0, +2.0] deg; Fig 3C), whereas the oculomotor range for the horizontal and vertical gaze directions runs between ±30-40 deg.
The literature has heavily debated the question whether these behavioural rules are due to a neural control strategy in the brain (e.g. [19][20][21][22][23]), or whether they are fully determined by the mechanical properties of the eye plant (e.g., [24][25][26]). For example, the pulling directions of the muscles, which eventually determine the eye's equilibrium orientation, could change in a particular way with the gaze direction, and automatically comply with Listing's law [24][25][26]. Certain filaments attached to the globe and surrounding the muscular tissue would function as pulleys [25], which not only prevent the muscles from excessive side-slip along the globe during eccentric viewing [27], but could also serve to mechanically implement Listing's law [24][25][26]. In support of this, [28] and [29] demonstrated that microstimulation of the abducens nerve, which innervates the LR muscle, generates eye movements that comply with Listing's law. These results indicate that, at least for the horizontal muscle system, the LR/MR pulleys aid in the mechanical implementation of Listing's law.
Clearly, mechanical and geometric properties of the extra-ocular muscles largely determine the direction of passive and active forces acting on the eye, and therefore on its instantaneous orientation. Note, however, that LL is only valid for eye movements with the head erect and gazing at far targets. For vestibular-evoked eye movements [30], head-unrestrained gaze shifts [31,32], near viewing [33,34], and microstimulation in several premotor brainstem areas [18], Listing's law is violated, and Donders' law applies. Although it has been hypothesized that the muscle pulleys could be actively controlled to shift the muscle's action dynamically, so as to accommodate programmed violations of Listing's law [35,36], neurophysiological evidence [37] as well as theoretical arguments [38] have rendered this hypothesis unlikely. Further, it should also be noted that the saccadic system is not perfect: often, the eye spontaneously jumps out of Listing's plane by a small saccade, leading to >2 deg of cyclotorsion (e.g., Fig 3C). Interestingly, the next spontaneous saccade is then endowed with an opposite torsional component that typically brings it closer to the plane, thus preventing the eye from walking away from LP during free gazing behaviour [23]. This indicates that saccades have three, rather than two degrees of freedom, and that the control of ocular cyclotorsion may possibly involve the cerebellum [18].
In our robotic prototype, the pulleys were kept head-fixed, so that a potentially active pulley mechanism is not considered. Instead, we here analyze the potential relationship between the emergence of Listing's law and the properties of the oculomotor plant on the basis of controlling a simple physical model.

Main sequence nonlinearity
Apart from the restrictions on the eye's cyclo-torsion, primate saccadic eye movements also display a number of characteristic dynamic behaviors. First, saccades obey stereotyped kinematic relations, known as the main-sequence [44]: saccade duration increases with amplitude, and the peak eye velocity saturates at large saccade amplitudes. Further, because the acceleration time of saccades is approximately constant (Fig 4A), the increase in saccade duration is primarily due to its deceleration phase, and may betray an internal dynamic feedback mechanism [45]. As a result, the skewness of saccade velocity profiles increases with saccade duration [43]. These properties are summarized as follows: where D is saccade duration; ΔR: saccade amplitude; V pk : peak eye velocity; T accel : time to reach V pk ; S: skewness of the velocity profile, and a, b, V 0 , α, c, d, k parameters. The two final relations hold for fast, as well as slow (drug-affected, fatigued, non-visual) saccades (Fig 4B  and 4C). Together, the relations of Eq (1) indicate the presence of a nonlinear dynamic (feedback) controller in the saccadic system [46]: for a linear controller, saccade duration and shape should be independent of amplitude, and the peak eye velocity should increase proportionally with amplitude. Because the oculomotor neurons are known to carry these properties in their bursting activity, the main-sequence behavior is not explained by a nonlinearity in the plant. Models of the saccadic system have therefore proposed nonlinear controllers at the level of the saccade-velocity activity carried by burst-neurons in the pons (for horizontal saccades [45]), and riMLF (for vertical/torsional saccades [41]).

Straight trajectories
Oblique saccades have approximately straight trajectories, which means that the horizontal, vertical (and torsional) velocity components are scaled versions of each other. Thus, if the signal to the plant for a saccade with amplitude ΔR, is encoded as a 3D angular velocity command, ω ΔR (t) = [ω x (t), ω y (t), ω z (t)] T (see below), the following holds: with α 2 + β 2 + γ 2 = 1. However, because of the saturating nonlinearity in the amplitude-velocity relation (Fig 4B), the cross-coupling in the component velocity profiles is not trivial if the horizontal, vertical and torsional velocity generators would be independently driven by their own dynamic displacement errors, Δz and Δy, and Δx, respectively: with f(�), g(�) and h(�) nonlinear.
For example, if a Δz = 10 deg horizontal component is part of an oblique saccade in LP, given by [Δx, Δy, Δz] = [0, 20, 10] deg, its duration (and velocity shape) matches the larger 20 deg upward vertical component, and therefore has a longer duration (and skewness) than a purely horizontal [0, 0, 10] deg saccade. This phenomenon is called component stretching, and constitutes yet another nonlinear property of the saccadic system. However, if the component velocity generators would all be driven by a common vectorial velocity command, ω ΔR (t), a simple trigonometric scaling of this signal ensures straight saccades and component stretching. This scheme has been introduced as the common-source model by [47].
Interestingly, neurophysiological evidence has suggested that the main-sequence nonlinearity ( Fig 4B) may reside upstream from the saccade-component velocity generators, at the level of the midbrain Superior Colliculus (SC). The SC could thus serve as the common nonlinear vector-velocity generator for saccades [48,49]. Furthermore, microstimulation at the SC motor map does not lead to violations of Listing's law [22], which suggests that the law is implemented downstream from the stimulation site, either involving brainstem-cerebellar neural circuitry [2,18], or at the level of the oculomotor plant (see above; [24][25][26]).
The primary objective of the present paper is to investigate whether a feedforward optimal control strategy for saccades, applied to a simplified 3D model of the eye, as illustrated in Fig 2, could account for all observed kinematic properties of saccades (Listing's Law, main-sequence nonlinearity, and straightness of the trajectories), and how the extra-ocular muscle geometry relates to Listing's law and the orientation of Listing's plane. Although several studies have investigated potential optimal control principles that could underlie the generation of saccades [50][51][52], these studies were all confined to the control of horizontal gaze shifts, and therefore did not address the 3D kinematics and dynamic problems that were central to the present study.

3D eye orientation
When an extra-ocular muscle contracts to exert a force, a torque is applied on the eyeball, which causes it to rotate around its center about an axis that is perpendicular to the plane defined by the force vector and the moment arm from the eye's center to the muscle insertion point. To describe changes in eye orientation in 3D, it is convenient to define the primary position as the origin of a head-fixed right-handed reference frame, from which any new eye orientation can be reached by a single-axis rotation. In this frame of reference, a positive torsional rotation (roll) is clockwise (rotation about the frontal +x axis), a positive vertical rotation (pitch) is downward (lateral +y axis), and a positive horizontal rotation (yaw) is leftward (vertical +z axis). There are many different ways to represent 3D rotations, but the most frequently employed ones in the oculomotor literature are unit quaternions and so-called Euler-Rodrigues rotation vectors (e.g., [6]). Here we will only summarize their relevant properties. For details, the reader is referred to [6].
Briefly, any single-axis rotation of a rigid body around the unit axisn ¼ ½n x ; n y ; n z � T by angle ρ can be parametrized by a unit quaternion (a 4-parameter complex algebraic object with length one): with q 0 = cos(ρ/2) the (real) quaternion's scalar part, I = [i, j, k] are the imaginary numbers for the (x, y, z) coordinates satisfying the relations ij = k, jk = i, ki = j, i 2 = j 2 = k 2 = ijk = −1, and q = sin(ρ/2) [n x , n y , n z ] T � [q x , q y , q z ] T its (imaginary) vector part. Consequently, multiplying two quaternions, p and q, yields: pq = p 0 q 0 − p � q + [p 0 q + q 0 p + p × q] � I. It is convenient to introduce the Euler-Rodrigues formulation [16,17], here denoted as the rotation vector, by: with r À 1 ¼ À r ðinverseÞ In this latter parametrization, Donders' law and Listing's law (see above) have the following formulations: where the latter defines Listing's plane, expressed in the primary frame of reference (see Fig 3B). It can be shown that the single-axis rotation vector from an initial eye position, r A , to any other final eye position, r B , is then described by the angular-velocity saccade vector (up to Oðr 3 Þ) [17]: with d � r B − r A (difference vector). Thus, if the initial and final eye positions are both in Listing's plane (r A,x = r B,x = 0): whenever angle (r A , d) 6 ¼ {0, π}. From Eq 9, the angular-velocity axis of the eye, ω k s AB , can be tilted out of LP by any η 2 [−ρ/2, +ρ/2], depending on the angle between the saccade difference vector, d and the initial eye orientation, r A . This geometric property is known as the half-angle rule [39].
Finally, from the quaternion relation between the angular velocity, ω � ω � I, the eye orientation, q, and the rate-of-change of eye orientation, _ q (the coordinate velocity; e.g. [19]): Expressed as rotation vectors, it follows that, up to Oðr 3 Þ: The rotation vector has some nice algebraic properties, which may also be useful for neural representations. E.g., the trajectory of a single-axis rotation for a saccade in LP from A to B is well approximated (up to Oðr 3 Þ) by a straight line in LP between the initial and final eye-position rotation vectors [17]: where σ(t) 2 [0, 1] (monotonic increase).

The eye-plant model
The methodology followed in this paper can be summarized as follows. First, we created a non-linear simulator model in Matlab's Simulink, which implemented Euler's equations for rotations of a rigid body, and the noncommutative 3D rotational kinematics (described above) to our simplified mechanical model of the eye (see Fig 2). Input to the system is a time sequence of signals from three independent motors, u(t), that each couple in a feed-forward way to the three antagonistic pairs of eye muscles, modelled as elastic tendons. The model's output is the instantaneous rotation vector of the model eye, r(t), that parametrizes its current 3D orientation. To identify a linear approximation for this system for the optimal control algorithms described below, the non-linear model equations were run through Matlab's systemidentification algorithms with independent pseudo-random binary signal sequences (PRBS). PRBS is a signal often used in systems identification. It is a deterministic signal composed of rectangular pulses of fixed amplitude but variable duration that has similar properties to white noise (flat spectrum), thus excites equally all the frequencies of interest in the system. The amplitude range of these signals was chosen such as to obtain good approximations of the model's output on the range of the desired saccade amplitudes. The identified model was then used to compute the optimal input sequences to control saccades, by applying different sets of combined performance costs. Finally, the optimized inputs were provided to the nonlinear simulator to generate 1500 randomly distributed simulated eye movements, of which we analyzed the kinematic and dynamic properties.
The simulator. To model the eye we did not attempt to fully replicate the anatomical and physiological features of the individual muscles. Instead, we only included their most essential features: their insertion points on the globe and approximate pulling directions, their elastic properties (note that the pull is caused by a lengthening of the string, rather than by a physiological contractile element; the mechanical effect on the eye's orientation, however, is identical), dynamic friction (a viscous force, caused by the drag of the optic nerve and orbital tissues through the fatty layers), and the way in which the muscles will change their direction of action when the eye orientation changes.
The input for the simulator (the angular positions of the drivers) was provided as a 3D command u(t) = [u 1 (t), u 2 (t), u 3 (t)] T , which acted as three independent signals for the three antagonistic muscle pairs (see Fig 5).
The model's output corresponded to the resulting rotation vector, r(t) = [r x (t), r y (t), r z (t)] T , which represents the 3D orientations of the eye that produced the resulting rotational movement trajectory.
The rotational dynamics for this rigid-body mechanical system, expressed in the inertial frame, are governed by Euler's equation: where τ T is the total torque on the eye, τ i,el are the elastic torques applied by the six muscles to the eye at their insertion points, Q i , and τ fric is the total frictional torque felt by the eye (viscosity, τ dyn and static friction, τ stat , see below; we ignored the force of gravity); I(q) is the eye's The common driving signals, u 1,2,3 , are symbolized by the grey bars, that rotate around the red axes, thus providing an antagonistic contraction/relaxation force on the associated muscles (cables with same colours) from their end-insertion points at P i (green dots). The six muscle insertion points, Q i , on the eye (red dots) form the (elliptical) blue surface; C is the center of the globe at (0,0,0). The default arrangement is nearly symmetrical with respect to the horizontal plane of the eye at z = 0.
https://doi.org/10.1371/journal.pcbi.1008975.g005 moment of inertia tensor (which may change as the eye rotates about an arbitrary axis around its center, C, and the mass distribution with respect to the primary reference frame changes too), α is the 3D instantaneous angular acceleration vector (in rad/s 2 ), and ω is the eye's angular velocity (in rad/s). The cross-product term vanishes when I(q) = aI 3 (where I 3 is the 3x3 unity matrix). For our simulator, we verified that its contribution was negligible when compared to the angular acceleration term, as the relative rms power of |ω(t) × (I ω(t))| vs. |Iα(t)| remained <2 � 10 −4 for saccade amplitudes up to 60 deg (S1(D) Fig). We therefore neglected the cross-product term in the simulator.
The torques were calculated from the simplified geometry, illustrated in Fig 5. Table 1 provides the coordinates of the insertion points in the (x, y, z) coordinate frame of the model's eye ball (re. its center), the coordinates of the via points, X i , for the pulleys of the SR, IR, SO and IO muscles (the LR and MR were directly connected to their muscle endings at P i ), and the resulting normalized torques exerted by the muscles when the eye is at its equilibrium orientation, r = (0,0,0). The insertion points at which the driving signals u exercise their action are represented by P i . The changes of the insertion coordinates, from (P i,0,x , P i,0,y , P i,0,z ) ( Table 1) to (P i,x , P i,y , P i,z ), of the vertical recti and obliques (i = [SR,IR], or [SO, IO]), which are controlled by u 1 and u 3 , respectively, in radians) were determined by the following relations: with L M = 20 cm. For the lateral recti (i = [LR,MR], with L R = 28 cm, controlled by u 2 ) this yielded: The locations of the muscle insertion points on the eye, specified at Q 0,i in the primary position (Table 1), will change as a result of the eye's rotation (given by the rotation matrix, R) Table 1. Coordinates of the insertion points (in cm) on the mechanical model of the right eye at rest (Q 0,i ), and the craniocentric via points ('pulleys') (X i ) with respect to the (head-fixed) center of the eye. Here, LR and MR have no pulleys; they directly connect to driver 2. The cranial insertion points, P i,0 at rest are given in columns 8-10. Columns 11-13: the unit torque-vector components, as computed from Eq (18) with the eye in the resting orientation. Note that the SO, IO, SR and IR tendons all pull with a prominent cyclo-torsional (x) and vertical (y) component.
Because of the changes in the eye-insertion points, the muscles will change their lengths and exert elastic forces F i through Hooke's law: with l 0,i the rest length of muscle i, and l i � l 0,i the current length of the muscle. κ i is the elasticity of the muscle, and was given the same value of 6.0 N/m for all tendons.t i is the unit vector of the tendon's pulling direction from Q i to X i for SO, IO, SR and IR, and to P i for ML and LR. The elastic torques (Fig 6 and Table 1) are then determined by Besides the elastic forces, the eye also feels frictional forces, because of its rotation through the fatty embedding within the orbit. In the model, we included both a small static friction (by Newton's 3rd law: equal and opposite to the total exerted torque) and an angular-velocity dependent dynamic friction. The dynamic frictional force was approximated by: where χ dyn = 0.02 Nms/rad is the dynamic friction coefficient. This parameter, which relates to the viscosity of the plant, was set to make the eye behave as an over-damped system. For the static friction we took: with τ comp = τ el + τ dyn , � a small dynamic cut-off velocity, |ω| is the amplitude of the angular velocity vector, and χ stat = 0.006 Nm is the static friction coefficient. Once all torques are known, the angular acceleration, α, can be computed through Eq (13), given that the inertial tensor can be assumed constant for each small time step of the simulation (Δt = 10 ms). The angular velocity vector of the eye, ω (and hence, its instantaneous axis of rotation,n, Eq (5)) is obtained from the angular acceleration by time integration. The resulting new orientation of the eye, q (or, alternatively, its rotation vector, r; Eq (4)), is subsequently computed by integrating the coordinate-velocity _ q, which was obtained from q(t − Δt), and ω from Eq (10). When the new orientation, q(t), is calculated, it is used to update the coordinates of the insertion points on the eye Q i and to update the tensor of the moment of inertia by I(R) = RI 0 R T (with I 0 the moment of inertia tensor in the eye's primary reference frame, and R the rotation matrix associated with the quaternion q; see Fig 6). For a perfect symmetric sphere I 0 would be given by I sphere ¼ 2 5 ma 2 , with a its radius, and m the mass, and invariant to eye orientation. In our robotic prototype, it is given by: A simplified diagram of the simulator is shown in Fig 6. The simulator was implemented in Matlab's Simulink toolbox (version 2018b). It is important to note that the equations governing the simulator are nonlinear. This is already obvious from the nonlinear updates of the cranial insertion points (Eqs (14) and (15)). However, also the pulling directions of the elastics change as the eye rotates, causing the gains between inputs and outputs to depend on the operating point. Moreover, the 3D kinematic relations include vector products which are noncommutative, and contain trigonometric dependencies on the angles. This all means that the model can only be well approximated by a linear system around a given operating point, and outside that region is expected to produce large errors.
To check whether the model was mechanically capable to cover the full range of 3D eye movements, we simulated a large feedforward set of 21 rotational steps for each of the three motors over a range of [-35, +35] deg (i.e., a total of 21 3 = 9216 combinations), to evaluate the 3D oculomotor range of the model. Fig 7 shows the result of this simulation. Clearly, the horizontal, vertical, and torsional ranges of the model eye were quite comparable, and close to [-50, +50] deg in all directions.
Linear system identification. The system identification procedure [53,54] consisted of the following steps. As input for the identification we used a 180 s long Pseudo-Random Binary Sequence (PRBS) signal, uncorrelated between the three motor inputs. The first 120 s were used to train the system-identification algorithm (typically yielding >98% accuracy for all components), while the remaining 60 s of the data were used for testing the model predictions and validating the quality of the fit. We constructed a state-space model from our system, which represents the relation between its inputs and outputs through a set of n th order differential equations, comprised within a single first order n × n matrix difference equation. We consider that the dominant dynamics is given by the eye mechanical inertia (the motor electric and mechanical time constants are comparatively much smaller), so the state vector has intrinsically 6 dimensions (3 orientations and 3 angular velocities). In this way, the discrete dynamical system is represented by the following two matrix equations, with n = 6: where t is the discrete time variable. To find the optimal parameters, the matrices [A, B, C, E], of the model, we applied the subspace method ( [55]; implemented in Matlab's "n4sid.m" function), which is known to provide numerically stable state-space models for multivariable dynamical systems [54,56].
To quantify the quality of the identification, we calculated the normalized root-mean squared error (NRMSE) fitness value, which indicates how well the predicted model response matches the non-trained target data for the nonlinear simulation: where r is the true validation data output of the nonlinear system,r is the output of the identified linear system and � r is the mean of r.
Even though the predictions on the validation set for the linearized fit of Eq (22) were not perfect, the result was generally judged satisfactory. The torsional and horizontal components exhibited NRMSE values above 75% for the torsional component, and 90% or higher for the horizontal component, which is excellent. The vertical component, however, lagged somewhat behind with a best correlation of 40%. It should be noted, however, that the linearization will always produce errors, since the simulator is highly non-linear, especially for the strongly cross-coupled torsional/vertical muscle pairs. Indeed, correlations for r x and r y increased substantially, when the via points for SO and IO were positioned less laterally (not shown).
Optimal control. We next used the identified system, Eq (22), to find the optimal control commands u i that minimize a total cost functional. Several cost functions for the gaze control system have been proposed in past studies to account for the nonlinear main-sequence properties of Eq (1). An optimal control strategy that minimizes saccade duration and maximizes saccade accuracy has proved successful (speed-accuracy trade off; e.g. [50]). However, models on optimal gaze control have so far been confined to one dimension only (horizontal saccades), and therefore miss other important aspects of gaze control, like component cross-coupling, and Listing's law. Here, we evaluated the contribution of different cost terms on overall saccade performance.
In general, the optimal open-loop control problem is represented as the solution to the following constrained optimization problem: Our aim is to identify the simplest cost function (with the smallest set of cost terms) that best accounts for the following properties: • Control of 3D saccadic movements, such that all axes describing eye orientations lie in a plane (Listing's plane), according to Eq 7 • Nonlinear dynamic properties of saccades (their main-sequence relations, including skewness, straight trajectories, and the associated component stretching of velocity profiles).
Costs. The first three cost terms to be considered have also been proposed in earlier studies. They acount for (i) minimization of the gaze-direction error (i.e., maximize saccade accuracy), (ii) minimize the total energy spent by the input signals, and (iii) minimize saccade duration, D.
Accuracy. The accuracy cost minimizes the mean-squared difference between the target direction, ðr y ;r z Þ, and the actual 2D gaze direction at the end of the saccade, (r y,D , r z,D ), without putting any additional constraints on the amount of cyclo-torsion, r x . Moreover, as the saccade reaches its goal, the velocity and acceleration should be zero. Thus, the classical twodimensional accuracy cost is written as: Clearly, this cost term is essential for any optimization procedure, as without it the system's performance will not be bound to any particular input-output relationship. For some solvers, it may be better to include zero velocity and acceleration constraints in the cost function to penalize the norms of the velocity and acceleration vectors. Energy. The energy expenditure by the input signals was taken to be proportional to the total squared input velocity vector across time. The latter is proportional to the instantaneous changes in the input vector, du t = u t+1 − u t . In this way, the cumulative energy cost was expressed by the length of the time-series difference vector: with Δ a matrix consisting of blocks of 3x3 identities, such that DU ¼ : Duration. The duration term penalizes a long saccade duration, D. To implement a duration cost, we took the hyperbolic discount of reward at D, proposed by [57]: with β is the rate at which the value of reward decays.
Listing's law. Because the first three costs do not put any constraints on the 3D behaviour of the eye movement, we also considered costs that might enforce a behaviour according to Listing's Law (Eq (7)). We considered two costs that explicitly penalize deviations from Listing's plane: one that only incorporates the final 3D orientation of the eye, regardless the trajectory followed, and one that constrains the entire eye-movement trajectory to Listing's plane.
Thus, the cost that penalizes deviations from Listing's plane at saccade offset is determined by: while constraining the entire eye movement trajectory to Listing's plane was associated with the following cost: with R x = [r x,0 , r x,1 , � � �, r x,D ].
Force. Note that the functionals of Eqs (29) and (30) assume an a-priori representation of Listing's plane in the system. As such, the primary position is explicitly assumed to correspond to the straight-ahead eye orientation. However, if Listing's plane were to be an emerging property of the optimization, such an assumption should preferably not be imposed a priori. In the default mechanical model of Fig 5, and Table 1, the origins of the extraocular muscles at the back of the orbit were arranged approximately symmetrically around the horizontal plane of the eye (z = 0). However, if these origins were to be jointly shifted upward or downward, the orientation of the primary position (the normal vector to Listing's Plane) might change accordingly. This possibility would suggest that the orientation of Listing's Plane could be related to an (a)symmetry in the pulling directions around the horizontal plane of the eye with respect to straight ahead. We thus included a cost that would minimize the total tension exerted by the extraocular muscles on the globe at the final eye orientation (i.e., fixation at t = D), and hence, would minimize the total effort required from the neural inputs when the eyes are not moving. As an empirical approximation, we related the total force and eye orientation by the following quadratic form: This quadratic relation resulted from simulations with the model across the system's oculomotor range for different configurations of the origins and eye orientations. To estimate the 3 × 3 matrix H F , we gathered data across a large range of 3D eye orientations, like in Fig 7. The weight of the force cost was taken much smaller than the weight for the accuracy cost, i.e. λ F � λ A , because the objective was that the force minimization should only differentiate between saccades ending in the same vertical and horizontal desired orientations, and not deteriorate the accuracy of the final gaze-direction components r y,D and r z,D .
Although, in principle, 32 different cost functions can be constructed from any combination of the accuracy cost with the other five, we considered the six cost functionals in Table 2 as the most relevant candidates for the optimal control.
The optimal saccade commands U, and saccade durations, D (at 10 ms resolution), were found by solving the quadratic problems with linear constraints, outlined in Eq (24) and Table 2, through applying Matlab's "quadprog.m" optimization algorithm. Fig 8 (left) shows an example of the optimization for the J AED cost as function of the movement duration for a 10 deg horizontal saccade. In the simulations, we ran 1500 saccades in random directions. The saccade vectors were drawn from a Gaussian distribution for the vertical (r y ) and horizontal (r z ) components with a standard deviation of 15 deg. The first saccade started at the (assumed) primary position, r 0 = 0, and any subsequent saccade started from the new end position in a random direction and amplitude (Fig 8, right). On this set of randomly directed saccades, we implemented the different optimization strategies outlined in Table 2.

Results
We characterized the 3D behaviour of the model by analyzing the saccade trajectories, kinematics, and velocity profiles. To quantify the distribution of rotation vectors, we considered all data points, R, and analyzed their deviations from the best-fit plane. We analyzed the behavior for two different simulator configurations: one, in which the insertion points, P i were symmetrical with respect to the eye's equator (z = 0, like in Fig 5), and one in which we systematically varied all insertions together along the z-direction (the asymmetrical case).
For the symmetrical simulator, the primary position is expected to coincide with the straight-ahead orientation, so that Listing's plane is simply defined by r x = 0 for all eye orientations. To compare the different optimal control strategies, we analyzed the dispersion Table 2. Six different cost functionals, all optimized to control saccades by tuning the weights, λ α for the contributing terms. Each total cost includes the accuracy cost, but differs in the contribution from other kinematic terms. AD: accuracy and duration; AE: accuracy and energy; AED: accuracy, energy and duration; AEDL 1 : AED with Listing's law at saccade offset only; AEDL 2 : AED with Listing's law for the entire trajectory. The weights are the (two to four) free parameters for the optimal control.
(standard deviation) of the torsional component of the eye orientations with respect to the model's Listing's Plane: s r x .
To quantify the straightness of the resulting saccades (that is, to what extent they are generated as a single-axis rotation, and follow the rules of component stretching, described https://doi.org/10.1371/journal.pcbi.1008975.g008 Table 3. Results for the different optimization strategies specified in Table 2. n: normal to the best-fit plane through the rotation-vector data, and absolute angle between n and the straight-ahead direction, P = [1, 0, 0]; Width LP: standard deviation of the data around the best-fit plane. MS: nonlinear main-sequence data fit on all saccade data: V 0 : asymptotic velocity (deg/s); α: angular constant (in deg); r CS : component stretching correlation of v pk,ΔH=8 o(F) vs. cos(F) (cf. Fig 11B). in the Introduction) we calculated the correlations between the vertical and horizontal velocity profiles, represented by � r v y ;v z , with v ¼ _ r. For a straight saccade, this value should be close to one. Finally, we analyzed the main-sequence properties of the saccades: the relations between duration, peak velocity and amplitude of saccades (Fig 4), and the stretching of the vector components for oblique saccades. Because the temporal resolution of the simulations was set at 10 ms, we did not calculate a measure for the skewness of velocity profiles, as it would be too coarse. We analyzed the dynamic properties for the entire set of 1500 saccades, as well as for three selected sets of 12 saccades made in three different directions (horizontal, oblique (at a 45 deg angle) and vertical), all starting from the primary position.
In what follows, we present the results for the optimal control of J AED (no constraint on ocular torsion), of J AEDL 1 (constraining the final eye orientation to Listing's Plane), and of J Force for nine asymmetric configurations (without putting an explicit constraint on ocular torsion). The results for the other optimization strategies (AE, AD, and AEDL 2 ) are summarized in Table 3, and in S2-S5 Figs. Listing's plane analysis. As can be seen in Fig 9A and 9B, the model eye looked all across the visual field (about 50˚range both in horizontal and vertical gaze directions in the yz plane). The best-fit plane, r x = ar y + br z , through the 3D eye-movement trajectories had a normal vector ofn ¼ ½0:955; À 0:266; À 0:129�, with the data scattering around the plane with a standard deviation of σ x = 0.043 half-radians (4.88 deg). The best-fit plane is therefore rotated with respect to the frontal (yz) plane, along the y-and z-axes. This can be appreciated by the slight forward tilt of the data cloud in the (xz)-view (panel B), in which the torsional scatter is nearly twice as large (10.95 deg; cf. Fig 3C). Clearly, the scatter around the best-fit plane is much smaller than the torsional oculomotor range of the model, illustrated in Fig 7. The axis and angle that will rotate the normal to the plane,n (which is expressed in the eyecentered reference frame), onto the true primary direction, P = (1, 0, 0) of Listing's reference frame, is given byn P ¼n � P ¼ ½0; À 0:437; 0:900�, and r P ¼ arccos ðn � PÞ ¼ 17:2 deg. Applying the associated rotation matrix, R P , to the rotation-vector data of Fig 9 then yields the eye orientations expressed in Listing's frame of reference:

Minimizing AED
The result of this coordinate transformation is shown in Fig 10 for the (x, z) view. Trajectory analysis. The simulations also showed that saccade trajectories were approximately straight, as the correlations between the horizontal and vertical velocity profiles had their mode very close to 1.0 (histogram, Fig 9D). This property emerges from the simultaneous optimization of the duration and energy costs. A straight trajectory thus ensures a movement with the shortest possible duration, at the lowest energy cost.
Main sequence analysis. The main-sequence plot of all saccade data is presented in Fig  9C. Note that there is a considerable variability in the peak vectorial eye velocities, which is due to the fact that saccades started all across the oculomotor range, and were pooled for all movement directions. Because of the mechanical properties of the eye muscles, it may be expected that the peak eye velocity will be position-and direction-dependent. To illustrate the direction-dependence on the saccade dynamics we ran the model once more to elicit saccades from straight-ahead in three selected directions: purely horizontal (red) dots, purely vertical (blue) and oblique (45 deg re. horizontal; green). The dashed curves show the best-fit exponential functions, described by V pk = V 0 (1 − exp(−ΔR/α)).
Note that the horizontal saccades were by far the fastest, while the vertical saccades were much slower. While almost all the torque applied to the eye by the horizontal muscles is converted in horizontal motion, the vertical and oblique muscles divide their torques along the vertical and cyclo-torsional components (see Table 1). Furthermore, because of the pulleys, the

PLOS COMPUTATIONAL BIOLOGY
force vectors applied to the eye by the vertical and oblique muscles are less orthogonal to the moment arm than the horizontal recti. Together, these factors reduce, on average, the effective torque, hence the peak velocity, in the vertical direction. The velocities of the oblique saccades tended to fall between these two extremes.
In line with the nonlinear main-sequence properties, and the relatively straight saccade trajectories, the model saccades also exhibited strong component cross-coupling. Fig 11 shows how the shorter component of an oblique saccade (here kept fixed at 8 deg in the horizontal direction) is systematically stretched in duration (and lowered in its peak velocity) as the orthogonal vertical component increases in size from 0 to 30 deg (i.e., saccade directions from 0 (horizontal) to nearly 80 deg (upward)). According to the common-source vectorial pulsegenerator model of [47], the small-component's peak velocity declines as the cosine of the saccade direction angle with the horizontal meridian (dashed line). The J AED model's result is quite close to this prediction. Note, however, that the common vectorial drive, as proposed by the common-source model, is not implemented in the hardware of our mechanical analogue,  as the three driving inputs, u = [u 1 , u 2 , u 1 ] were taken essentially independent. Hence, the observed cross-coupling between the eye-movement components is an emergent property that resulted from the constraints imposed by the optimal control. Zero torsion at saccade offset 3D kinematics and dynamics. Although the saccade dynamics resulting from J AED minimization closely follow experimental data, the 3D orientations of the eye still allowed for a wide range of torsional values. Thus, this optimization criterion by itself is not sufficient to explain the emergence of Listing's law. To check whether expanding the accuracy criterion for gaze, which is essentially 2D, to a 3D analogue by requiring that the saccade should have (close to) zero torsion at the end of the movement, we reran the optimization by including the endpoint Listing cost in the total cost function: J AEDL 1 . The result is shown in Fig 12A and 12B. Now, the saccades follow Listing's law, as the standard deviation of all eye orientations around the best-fit plane (which now has its normal vector very close to the P = [1, 0, 0] direction), has reduced to only 1.7 deg. This value is close to what is observed in monkey and human data [39] (cf. Fig 3).
As can be appreciated, adding the zero-torsion cost did not influence the shape of the saccade-velocity profiles, the straightness of their trajectories, or their main-sequence properties (panel C). There was nearly full component stretching in the oblique saccades (panel D). Thus, the AED minimization by itself forces the eye to behave as a single-axis rotation, and the constraint on the final ocular torsion then suffices to keep the entire eye trajectory near Listing's plane. Including a cost on eye torsion would suggest that the saccadic system has an internal representation of its torsional position state. As argued in the Introduction, experimental evidence seems to supports this idea. For example, any deviation from Listing's plane, even as  small as one to two degrees, is typically corrected by the next saccade [18]. Moreover, in more natural 3D gaze-control behaviors, like in eye-head coordination, data suggest that the system appears to plan a gaze trajectory for which the ocular torsion is zero when the entire gaze shift is over [32]. Such a strategy would be in line with a cost that includes eye torsion at the end of the saccade.

Minimizing fixation-force
Symmetrical case: Δz = 0. Yet, explicitly forcing the eye to have zero torsion at the end of the saccade also presupposes that the orientation of Listing's plane coincides with the straightahead direction. This, however, is not supported by experimental data (Fig 3). Thus, we

PLOS COMPUTATIONAL BIOLOGY
wondered how the orientation of Listing's plane would depend on the pulling directions of the extra-ocular muscles. We reasoned that also the forces exerted by the muscles on the eye when it is not moving, are associated with a certain cost. At steady fixation, the brain might want to minimize the total force required to maintain that eye orientation.
In Fig 13 we show the results when including a force cost at fixation for the mechanical prototype shown in Fig 5 and Table 1. Interestingly, the inclusion of force minimization in the optimal control led to good results for all aspects of saccade control. The scattering around Listing's plane resulted to be very small: a width of only 1.55 deg, which is even better than for the J AEDL 1 cost functional. The trajectories were approximately straight, with near-complete  Table 1; column P z 0 ). To further illustrate the appropriate saccade dynamics, Fig 14 shows that for all (slow and fast) saccades, regardless their trajectories and starting points, a tight relation emerged for the saccade amplitude vs. the product of its peak velocity with its duration. The simulated data show a larger variability than is typically observed in the human data (Fig 4C), which is mainly due to the discrete nature of the saccade-duration estimates (here taken at a 10 ms resolution).
Asymmetrical cases: Δz 6 ¼ 0. To verify the influence of the muscular geometry on the orientation of Listing's plane, we varied the relative position of the muscles' cranial insertion points with respect to the equatorial plane of the eye. We shifted all insertion points either upwards or downwards with respect to the horizontal ocular midplane, by Δz, ranging from  Table 3 summarizes the results for all optimization strategies. Although the J AE control strategy (no duration constraint) seems to result in a 3D behavior constrained by Listing's law, it does not generate normal saccade dynamics. Typically, the main sequence resulted to be close to linear, as evidenced by the large angular constant of the saturating exponential curve. The saccade durations for this strategy were all D = 220 ms, as no cost was associated to the saccade speed. As a result, saccades were slow, in order to minimize the total energy (proportional to kΔUk 2 ) spent by the inputs. Moreover, although saccades were straight, they did not display any component stretching, as the velocity profiles for the horizontal components of oblique saccades did not change with saccade direction (r CS = -0.04, Table 3; see also Supporting Information, S2-S4 Figs). This is to be expected from a linear main sequence: in that case, component stretching is a simple consequence of linear vector decomposition.

Results for all optimizations
Abnormal dynamics were also obtained for the J AD functional (no energy constraint): again the saccade durations were constant, but now at D = 100 ms, and the model produced very high saccade speeds (up to 1000 deg/s), as the amount of energy spent by the input signals was not penalized. Component cross-coupling in oblique saccades was observed, but in a completely opposite way as seen in Fig 11: as the vertical component increased in size, the peak velocity of the constant-size horizontal component also increased (rather than decreased) with the saccade direction (r CS = -0.98, Table 3). To keep the final amplitude at 8 deg, however, the horizontal component initially made a large overshoot, and then returned to the final target location during the remaining vertical movement. As a result, the saccade trajectories were In contrast, imposing torsional constraints on the saccade endpoints, or on the entire saccade trajectories, as well as minimization of the fixation force for all muscle geometries led to both Listing's law, and to normal nonlinear saturating saccade dynamics. In all cases, saccades were nearly straight, and they displayed appropriate component cross-coupling as expected for single-axis rotations.

Summary
This work accounts for the dynamic and kinematic behaviors of saccadic eye movements and eye orientations in three dimensions, by simulating a simplified mechanical model of the eye,

PLOS COMPUTATIONAL BIOLOGY
driven by three independent motor signals, and subjected to different optimal control strategies. Our model explains the major properties of saccade kinematics and dynamics in 3D with a minimum of prior constraints on its biomechanical properties and neural control: directiondependent nonlinear main-sequence behavior, component cross-coupling, 3D single-axis rotations of saccades, straight trajectories in Listing's plane, and Listing's law.
In our model, all six extraocular muscles were given the same elasticity (and kept constant), the viscous force on the eye was taken proportional to eye angular velocity, and the eye's inertial tensor depended slightly on instantaneous 3D eye orientation. The direction of the torque vector, exerted by each tendon, was determined by the shortest path from the muscle's insertion point on the globe, Q i , to its craniocentric via point, X i (through passive pulleys for the vertical and oblique muscles). Below, we will argue that these simplifications were not essential to obtain the resulting eye-movement properties. Our model provides novel suggestions for the potential control principles underlying the emergence of Listing's Law and the mainsequence in our robotic prototype, and it may serve as a starting point for more realistic biomimetic implementations in humanoid robotic systems.
Although previous studies have applied optimal control principles to explain human (saccadic) oculomotor behavior, they have so far been confined to horizontal saccades only [13,[50][51][52][57][58][59]. These studies suggested that the nonlinear main-sequence behavior of saccades (Fig 4) may result from a neural control strategy that aims to optimize speed-accuracy tradeoff, rather than from mere saturation of neuronal peak-firing rates at the level of the brainstem burst generator as proposed by [45]. Indeed, more recent single-unit recordings from the midbrain SC provided supporting evidence for such an optimal control by revealing synchronous firing properties of Superior Colliculus cells in the neural population that encodes the vector of the saccade goal, and its trajectory [48,49].

Emerging properties
Our 3D mechanical model accounts for a number of emerging saccade properties, that have so far not been explained by previous models: i. the nonlinear main-sequence holds for saccades in all directions, and its direction-and initial position-dependence results from the geometric arrangement of the six eye muscles. For example, purely vertical saccades from the straight-ahead initial orientation require a precise coordination of the vertical recti and oblique eye muscles, which have their optimal pulling directions not aligned with the cardinal (vertical and cyclo-torsional) Cartesian axes (see Table 1), and had less-effective, non-orthogonal, moment arms. As a result, vertical saccades of the model were slower than horizontal saccades, which involved an optimal orthogonal activation of the horizontal recti (e.g., Fig 9).
ii. Despite the direction-dependence of the main sequence, the model produced a considerable amount of cross-coupling between the components of the instantaneous eye-velocity vector. This resulted in trajectories in which both velocity components had very similar durations, and were roughly scaled versions of each other (Fig 11A).
iii. Because of this strong component cross-coupling, saccade trajectories were approximately straight in all directions. As a result, the smaller component of the saccade vector was stretched in time, and reached a lower peak velocity than when executed in isolation. This effect could be well described by the cosine (horizontal component, as in Fig 11B) and sine (vertical component; not shown) projections of the vectorial velocity vector, as predicted by the common-source model [47]. This conceptual model assumes that saccades are driven by a central vectorial pulse generator, which represents the vectorial velocity of the eye, and that the horizontal/vertical/torsional saccade components are jointly derived from this common drive by vector decomposition. As our mechanical model is driven by three essentially independent motors (Fig 5), the apparent common vectorial velocity drive is implicitly encoded by their outputs as a result of the optimal control constraints.
iv. By including the minimization of static total force on the eye, Listing's law emerged without any additional assumptions regarding the eye's cyclo-torsional state, or with the inclusion of precisely positioned pulleys for the horizontal and vertical/torsional muscle pairs. Moreover, this force cost could also account for the orientation of Listing's plane with respect to the head, and hence, for the direction of the primary position with respect to the oculomotor range (Fig 3). Although suggestions in this direction have been made before in the literature [16,24,27,60], our model explicitly demonstrates this potential relationship for the first time.

Saccade dynamics
The dynamic properties of the saccades require the simultaneous optimization of energy, duration and accuracy (J AED minimization; Fig 8). Putting a constraint on the saccade duration in the form of a hyperbolic discount leads to the saturation of saccadic peak eye velocity for large amplitudes, and the affine relation between saccade duration and amplitude. Adding a penalty on the total energy consumption of the 3D input drive is cause for the strong cross-coupling between the components. This cross-coupling ensures straight saccades (and thereby the shortest path taken by the eye), and can be effectively envisioned to promote a single-axis rotation of the globe from starting to end orientation. Dropping either the penalty on saccade duration (J AE ) or energy (J AD ) disrupts the main-sequence and cross-coupling properties of the saccades (Table 3; Supporting Information, S2-S4 Figs). However, it should be noted that the current simulations did not incorporate the effect of signal-dependent (multiplicative) noise in the input driving signals. It can be demonstrated that speed-accuracy trade-off in the presence of signal-dependent noise will show up in a similar quadratic way in the total cost as the energy term in the current formulation (Eq 26), so that the energy term in the cost functional may be dropped when signal-dependent noise is included in the model ( [58,61]; unpublished results from our lab).

3D saccade kinematics
The classical cost functional that jointly minimizes the mean accuracy of (2D) target acquisition, input energy expenditure, and movement duration (J AED ) could not readily account for the 3D behavior of saccades, as the torsional excursions of the model eye resulted to be too large, at about ±10 deg. Thus, to obtain a behavior that would be better in line with Listing's law, the accuracy constraint had to be extended by either explicitly penalizing the amount of ocular cyclo-torsion at saccade offset, or (but perhaps trivially) by constraining the entire trajectory to Listing's plane. Because the AED criterion by itself results in single-axis rotations of the eye, putting an additional constraint on the final eye orientation in 3D would thus ensure that the entire trajectory remains in Listing's plane, by virtue of Eq 12. This extra constraint hardly affected the dynamic behavior of the saccades. The requirement that the entire saccade trajectory should stay in Listing's plane, however, led to much slower saccades than the other cost functionals of Table 2, as the asymptotic velocity was only 440 deg/s. Inclusion of a Listing cost only at saccade offset generated saccades that could reach velocities of nearly 600 deg/s (Table 3; see also Supporting Information, S5 Fig).
Yet, the explicit requirement of zero cyclo-torsion presupposes that the primary direction is known, and aligned with the straight-ahead direction, i.e., the center of the oculomotor range. Experimental results, however, indicate that the primary direction may be located more eccentric, in the upward viewing direction (see Fig 3, for an example). Although the precise reason for this off-center direction is unclear, our model simulations with the asymmetrical muscle insertions indicate that the geometrical arrangement of the extra-ocular muscles could be directly related to this phenomenon. Interestingly, by adding a static force cost to the overall classical AED cost functional, Listing's law became an emerging property of the system, but now with the primary direction systematically depending on the muscle geometry (Figs 15 and  16). Note that the human eye generates about 3-4 saccades/s, with an average duration of 80-100 ms, which means that 60-70% of the time the eye is not moving, but kept at the eye orientation obtained after the previous saccade. It would therefore make sense for the system to minimize the total effort (i.e., "optimize the skill") needed to keep the eye in any particular orientation.
Our simulations thus suggest that Listing's law (Fig 3) and the main-sequence behavior of saccades (Fig 4) could be part of a joint optimization strategy of the oculomotor system that aims to minimize movement duration and fixation effort when collecting sensory evidence about the environment. In dealing with this problem, the system has to overcome a sluggish (overdamped), anisotropic plant with quite complex static and dynamic mechanical properties [25-27, 35, 38, 62], as well as the noncommutative nature of sequential rotations [3,6,16,19,37,39].
Note that Donders' law and Listing's law could also help to optimize visual-sensory functions, like stereoscopic depth perception and the perception of (binocular) visual line orientations seen from different eye positions [21,33,34]. The pulley system of the horizontal recti [25,26] might then also function to facilitate joint binocular control. If so, their plane of action, and even the location of the primary position, could then perhaps reflect an adaptive response of the oculomotor system to the statistics of natural binocular visual inputs [63].

Model simplifications
Clearly, the mechanical model of Fig 5 is a coarse simplification of the true complexity of the human eye: the actual properties of the muscles were modelled as simple linear springs with a fixed elasticity, taken to be equal for all six tendons (cf., [27,64]). Moreover, the paths taken by these muscles to exert their forces on the globe were not constrained by (static or dynamic) pulleys [25,26,29,38]. These mechanical simplifications for sure have an influence on the exact motor commands needed to control the plant. This is already demonstrated by the effects of the different geometrical arrangements of the muscles with respect to the eye's equatorial plane (Fig 15), but also by adding head-fixed pulleys for the horizontal recti in our simulator (results not shown). Yet, the major saccade properties: their main-sequence behavior, straight trajectories, component cross-coupling and Listing's law, all emerge for these different geometries, albeit in slightly different ways (e.g., as a tilt of LP). The concept is illustrated in the schematic of Fig 17. More realistic models of (nonlinear) extraocular eye muscles [27,64], and incorporation of a pulley mechanism [24][25][26], will obviously change the model eye's 3D transfer characteristic, but the learning algorithm that uses the eye's measured outputs (3D orientation and angular velocity) in combination with the set of optimization constraints, will result in appropriate (albeit different) commands to cope with more complex plant properties. The relative RMS power of the two terms (calculated over the same 300 ms as in C) as function of saccade amplitude. Its value increases with amplitude up to about 0.02% for the largest saccades. As a result, the latter term was ignored in our simulations. (TIF) S2 Fig. (A,B) 3D oculomotor behavior and (C,D) saccade dynamics, resulting from minimizing the AD cost. Although the 3D eye orientations are close to Listing's law, the saccade dynamics (main sequence and component cross-coupling) are abnormal. The saccade peak velocity (C) doesn't seem to saturate, as the exponential fit yields a non-physiological asymptote near 10 4 deg/s, with a angular constant close to the oculomotor range. (D) In the fixedcomponent oblique saccade test (as in Fig 9), the 8 deg horizontal component develops a large overshoot (with all durations the same, at about 100 ms), which is followed by a long return phase. This results in goal-directed, but curved saccade trajectories. See also Table 3 for numerical details. (TIF) S3 Fig. (A,B) 3D oculomotor behavior and (C,D) saccade dynamics, resulting from minimizing the AE cost. Also for this control strategy the saccades are abnormal: (C) the peak velocity increases linearly with saccade amplitude (but saccades remain much slower than for the J AD strategy, with the fitted asymptote at around 650 deg/s), and absence of component stretching (D), as all horizontal saccade components have the same, slow, velocity profile (all Fig 17. The total model encompasses two linked (nonlinear) components: the sequence of motor commands, U t , and the 3D eye plant, H. The plant's output (described by 3D eye orientation, r t , and angular velocity, ω t ) is used to tune the cost function, J Force (r, ω), that drives the motor outputs. The system's behavior follows the dynamics and kinematics of real saccades (insets, right). Any change in the plant's parameters will require corresponding changes in the motor outputs, u t . r D : 2D target location. The thin feedback arrow represents the training phase during which eye movements recurrently train the motor outputs to minimize the imposed cost.  Fig 11), indicating that saccade trajectories are curved. See also Table 3. (TIF) S5 Fig. A,B) 3D oculomotor behavior and (C,D) saccade dynamics, resulting from minimizing the J AEDL 2 cost. As expected, saccades follow Listing's law (width of the plane is 2.25 deg), but peak velocities (C) of the saccades are much lower than for the strategies constraining AED and AED with final eye position in LP. (D) Saccade trajectories are straight, as the correlations between horizontal and vertical velocity profiles are high (see also