Abstract
Our movements are variable, but the origin of this variability is poorly understood. We examined the sources of variability in human saccadic eye movements. In two experiments, we measured the spatiotemporal variability in saccade trajectories as a function of movement direction and amplitude. One of our new observations is that the variability in movement direction is smaller for purely horizontal and vertical saccades than for saccades in oblique directions. We also found that saccade amplitude, duration, and peak velocity are all correlated with one another. To determine the origin of the observed variability, we estimated the noise in motor commands from the observed spatiotemporal variability, while taking into account the variability resulting from uncertainty in localization of the target. This analysis revealed that uncertainty in target localization is the major source of variability in saccade endpoints, whereas noise in the magnitude of the motor commands explains a slightly smaller fraction. In addition, there is temporal variability such that saccades with a longer than average duration have a smaller than average peak velocity. This noise model has a large generality because it correctly predicts the variability in other data sets, which contain saccades starting from very different initial locations. Because the temporal noise most likely originates in movement planning, and the motor command noise in movement execution, we conclude that uncertainty in sensory signals and noise in movement planning and execution all contribute to the variability in saccade trajectories. These results are important for understanding how the brain controls movement.
Introduction
The movements that we produce are variable. This is apparent in sports such as golf, darts, and basketball in which it is crucial to have good control over one's movements. Even after many years of training, highly skilled sportsmen are unable to produce exactly the same movement several times in succession. From a neuroscientific point of view, movement variability is important because it can help us understand how the nervous system controls movement. First, movement variability forms part of the behavior that theories of motor control should explain. Second, because variability places a limit on how successful movements can be, it has been suggested that movements are planned such that the variability is minimized (Harris and Wolpert, 1998). For both of these arguments, it is crucial to know the sources of movement variability.
In theory, movement variability could originate from (1) the sensory information defining the target of a movement, (2) the central planning of a movement, and/or (3) the peripheral processes of movement execution that take place after motor commands have been planned. Surprisingly, different studies aimed at identifying the sources of movement variability reached very different conclusions. Whereas Gordon et al. (1994) and Churchland et al. (2006) concluded that the variability in arm movements originates mostly in central movement planning, van Beers et al. (2004) demonstrated that movement execution is the major source. Moreover, Osborne et al. (2005, 2007) showed that sensory information is the major source of variability in smooth pursuit eye movements. Although the major source could be different for movements of different body parts, the conflicting results for arm movements demonstrate how poorly the origin of movement variability is understood. An explanation for the conflicting conclusions could be that each study focused on one particular source while neglecting other sources that could have produced parts of the observed variability.
To get a better insight into the sources of movement variability, we examined the variability in human saccadic eye movements. Saccades are particularly well suited for this purpose, because, unlike most other movements, we have no voluntary control over their kinematic properties such as speed and duration (Becker, 1989). As a result, saccades are stereotyped movements whose kinematic properties are determined by the amplitude, a phenomenon known as the main sequence (Bahill et al., 1975). Establishing the sources of variability is therefore not hampered by variability in the subjects' voluntary control. Another advantage of saccades is that they are too quick to be corrected based on sensory feedback received during the movement (Becker, 1989).
Our approach is to analyze the full pattern of spatiotemporal variability of human saccade trajectories to targets in different directions and at different distances. The sources of the observed variability are then inferred via computational modeling of various types of noise in the sensory signals and motor commands. The analysis of the full pattern of spatiotemporal variability is crucial in this approach because it allows us to separate the contributions from different sources.
Materials and Methods
Subjects
Five subjects (two females, three males; 24–36 years of age) participated in this study. The author was one of them; the other subjects were unaware of the purposes of the study. The results of the author were not systematically different from those of the naive subjects. Four subjects had normal vision and one had corrected-to-normal vision, by wearing spectacles. This study is part of a research program that has been approved by the local ethics committee.
Apparatus
Stimuli were presented on a 40 × 30 cm computer screen with a resolution of 1024 × 768 pixels that was viewed from a distance of 67 cm. For each subject, the point midway between the eyes was aligned with the center of the screen. The two-dimensional orientation of both eyes was recorded at 250 Hz with an EyeLink system (SR Research, Mississauga, Ontario, Canada). The EyeLink system was used rather than the more precise scleral search coils because recent studies identified two effects of search coils that make them unsuitable for the present study. First, search coils on the eyes make saccades slower (by ∼8%) and last longer (by ∼8%) than in the natural situation in which no coils are present (Frens and van der Geest, 2002), and second, they increase the SD in the endpoints by some 20% (Smeets and Hooge, 2003). The video-based EyeLink system does not have these problems but measures natural saccade trajectories with a, for our study, sufficient accuracy and precision (van der Geest and Frens, 2002). Note that modern video-based eye trackers that can also record the torsional component of eye orientation, have insufficient accuracy for our purposes (Houben et al., 2006).
To minimize movement of the EyeLink cameras relative to the eyes, the subjects' heads were immobilized by a combination of a bite-board and a forehead support, and the EyeLink cameras were not mounted on the supplied headband, but, to avoid headband slippage, they were rigidly connected to the bite-board. Each subject chose her or his “natural” upright head orientation.
Procedure
In experiment 1, we tested saccades with different amplitudes, whereas experiment 2 examined saccades in different directions. Both experiments had a blocked design. In each block, subjects changed their gaze 75 times back and forth between two targets. A blocked design was used rather than randomly appearing targets to avoid that slow changes in saccade performance during an experimental session (Jürgens et al., 1981) would distort the estimated variability. One target, the starting position, was always straight ahead of the right eye; the position of the other target varied across blocks. Both targets were continuously present throughout the entire block. Stationary targets were used rather than jumping targets because saccades to stationary targets are more accurate (they undershoot less) and they are not preceded by smooth, anticipatory eye movements in the direction of the target (Lemij and Collewijn, 1989).
Stimuli consisted of black rings (0.33° outer diameter; 0.17° inner diameter) shown on a white background. The small inner disk was the actual target that subjects were asked to fixate.
A calibration block preceded each experimental block. In a calibration block, 25 targets, organized in a 5 × 5 grid covering the whole screen, were shown sequentially. Each target was shown for 2 s, giving subjects ample time to achieve accurate fixation. The eye orientation recordings of a calibration block were used to calibrate the EyeLink output for the experimental block that immediately followed it (for details of the calibration method, see below, Analysis).
Experiment 1.
Nine targets were used exactly to the left of the starting position at distances of 2, 4, 6, 8, 10, 12, 14, 16, and 18°. The time to initiate a saccade was specified by a metronome that produced beeps at a rate of 40 beeps per minute. Each experimental block therefore lasted 3 min and 45 s and followed a calibration block that took 50 s. Subjects thus had their heads immobilized for periods of ∼4.5 min. Between blocks, they could move their head and relax as long as they wanted but not <30 s. Blocks were presented in a random order, and all blocks were tested in a single session.
Experiment 2.
Targets were placed in 24 equally spaced directions at a fixed distance of 9° from the starting position. The metronome produced 50 beeps per minute and blocks were again presented in a random order. Each subject took part in two sessions of 12 blocks each that were run on different days. Both sessions took place ∼50 d before experiment 1. All other details were identical with those of experiment 1.
Analysis
Throughout this study, eye orientation was expressed as rotation vectors (Haslwanter, 1995). A rotation vector is a three-dimensional vector that represents the single-axis rotation that is required to bring the eye from some reference orientation (also known as primary position or primary orientation) to its current orientation. The direction of a rotation vector defines the rotation axis and its length equals one-half the tangent of the rotation angle. Because most readers will be more familiar with degrees to express the magnitude of rotations, results will be reported in degrees.
We represented the full three-dimensional eye orientation because this allowed us to estimate motor commands (see below, Estimation of motor commands). However, we recorded only the horizontal and vertical components. To obtain the third component (torsion), we assumed Listing's law to hold and we assumed Listing's plane (Haslwanter, 1995) to be oriented perpendicular to the right eye's line of sight when fixating the starting position. Control simulations were performed to estimate the effects of possible departures from these assumptions (see Results).
Calibration.
We did not use the standard EyeLink calibration procedure but instead used the following, more accurate method. During both the calibration and the experimental blocks, we saved the raw pupil coordinates of EyeLink. In the output of the calibration blocks, we selected 100 ms intervals during which the output was stable at the ends of the 2 s intervals for which each target was presented. We assumed that in this interval the subject fixated the target. Two-dimensional quadratic regressions were then made to find the relationship between the mean EyeLink output during these 100 ms intervals and the eye orientations required to fixate the targets. This produced calibrations with typical root mean square values of 0.3°.
Despite our attempts to minimize movement of the head relative to the cameras, the experimental data still displayed some drift. This was probably attributable to small head movements allowed by the soft tissues in the forehead. Even though this drift was slow (<1 min arc/s), we compensated for it to avoid misestimation of saccadic variability. We determined for each saccade from the starting position to the target the mean eye orientation in the 20 ms interval immediately before the saccade. This mean orientation was then subtracted from the entire trajectory of the saccade. This method was justified by plotting, for entire experimental blocks, all eye orientations immediately after the primary saccade and also the eye orientations 600 ms after the onset of the primary saccade (i.e., after possible corrective saccades were made). Whereas the former endpoint distributions were often not centered on the target and displayed a relatively large amount of variability, the latter distributions were much smaller and were well centered on the target.
Kinematic analysis.
Eye orientation data were filtered with a second-order, zero phase-lag Butterworth filter with a 125 Hz cutoff frequency. Cubic spline interpolation was then used to upsample the data to 1000 Hz. This did not change the values for the frames for which data were recorded, but only interpolated smoothly in between. This enabled us to estimate movement duration at a finer timescale, and it improved the accuracy of the peak velocity estimates.
Angular velocity Ω, a three-dimensional vector, was determined as follows (Haslwanter, 1995):
where r is the eye orientation (expressed as a rotation vector), ṙ is its temporal derivative, and × denotes the cross product. Temporal derivatives were calculated numerically using a three-point central-difference method.
Saccade onset and offset were determined using a velocity threshold of 0.5 s−1 (expressed as rotation vectors), which corresponds to 29°/s. Saccade onset occurred at the time of the last frame before the norm of the angular velocity vector exceeded this threshold; saccade offset occurred at the time of the first frame after the velocity fell below the threshold. The difference between saccade offset and onset times defined the movement duration M. Peak velocity vp was defined as the peak value of the norm of the angular velocity vector.
Analysis of variability.
We aimed to identify the sources of the spatiotemporal variability in human saccades by computational modeling of different sources of variability. The comparison of observed and modeled variability is only valid when both are computed from sets of (observed or modeled) saccades that have no suspicious outliers, because outliers can strongly influence measures of variability. The sets of observed saccades appeared to have outliers, possibly related to variations in the subjects' levels of attention and alertness. We rejected saccades whose starting position, peak velocity, time of peak velocity, movement duration, maximum perpendicular deviation (the maximum deviation of the saccade trajectory from a straight line between the start and end orientation), or end position differed >3 SDs from the mean of all saccades of that subject to that target. In addition, occasional saccades that had a tail in the velocity profile with a velocity <2 s−1 (∼115°/s) for >30 ms were also rejected. The calculation of variability was based on, on average, 67.8 saccades per target per subject in experiment 1, and 63.6 in experiment 2.
To analyze the variability in saccade endpoints, we determined for each subject for each target the covariance matrix of all (not rejected) endpoints. This variability was visualized by 95% confidence ellipses. Following van Opstal and van Gisbergen (1989), we decomposed the endpoint variability into two components. SD ςR represents the variability in the direction defined by the straight line between the starting location and the mean of all endpoints (“the variability in movement amplitude”), whereas ςφ represents the variability in the direction orthogonal to that (“the variability in movement direction”). The variability in all other quantities was expressed as an SD and/or as a coefficient of variation (CV) (the SD divided by the mean).
We also analyzed correlations between different variables. The results show that movement duration, amplitude, and peak velocity are all correlated with one another. In that situation, it is not correct to consider simple correlation coefficients between each pair of variables because these fail to take into account the interactions with the other variable. We therefore calculated partial correlation coefficients among these three variables.
Estimation of motor commands
Our model (see below) requires us to estimate the motor commands that were sent to the extraocular muscles to generate the observed saccades. Here, we explain how we estimated these motor commands.
We first derived mean trajectories for the movements that each subject made toward each target. For that purpose, we temporally rescaled all trajectories toward a target to their mean duration and resampled them at 1 ms using cubic spline interpolation. The mean trajectories were then obtained by averaging these resampled trajectories.
We next calculated the total torque, T, that was generated by the muscles, from the mean trajectories. Following Tweed and Vilis (1987), Schnabolk and Raphan (1994), and Raphan (1998), we assumed that the muscle torques must counteract the inertial, viscous, and elastic forces of the oculomotor plant as follows:
The first term on the right-hand side represents the inertia of the globe and equals the product of the eye's moment of inertia J and the angular acceleration. The second term represents viscosity, and is proportional to the coefficient of viscosity B and the angular velocity. The last term describes elasticity, which is proportional to stiffness K and rotation angle φ from the reference position, and has the direction of rotation axis n̂. We used a moment of inertia of J = 2.00 × 10−7 kg · m2 (Robinson, 1964). The viscosity and stiffness were derived assuming time constants of the overdamped human oculomotor plant of 224 and 13 ms (Robinson et al., 1986): B = 1.63 × 10−5 kg · m2 · s−1 and K = 6.87 × 10−5 kg · m2 · s−2.
We next calculated how this total torque was produced by the muscles. We treated antagonistic muscle pairs as single ideal muscles that can produce positive and negative torques, and whose insertions on the globe are orthogonal to each other (Quaia and Optican, 1998). The net torque generated by each muscle pair can then be found if the pulling directions of all muscle pairs are known. These pulling directions vary with the eye orientation because of the effects of muscle pulleys (Demer et al., 1995) and/or orbital fat (Schutte et al., 2006):
Here, τ is a three-dimensional vector that contains the torques generated by the three muscle pairs and RM is the matrix for a rotation of angle ½φ about axis n̂ (Raphan, 1998).
The final step is to calculate motor commands from muscle torques. Physiologically, the best definition of a motor command would be the entire set of signals sent to all extraocular muscle fibers. This, however, is not a useful definition here because we do not know all of these signals. Instead, we will use the equivalent of the sum of all the signals sent to a muscle as the definition of a motor command. Motor commands can then be calculated from torques by modeling the muscle as a first-order low-pass filter with a time constant of 10 ms (Harris and Wolpert, 1998; Tanaka et al., 2006). A low-pass filter is a simplification of the true complexity of extraocular muscles, but it captures the temporal filtering properties that are important for propagation of noise.
In terms of the well known pulse-step pattern of saccadic innervation, the commands we calculated correspond to the innervation up to the end of the pulse component. The step component that is still produced afterward to keep the eye in its position was not modeled because we considered the variability in the saccades only, not in possible slow movement afterward. We confirmed that the shape of our estimated motor commands resembles their recorded shape (van Gisbergen et al., 1982).
The major shortcoming of our procedure is that it cannot represent coactivation of antagonistic muscles. Our definition of noise, however, allowed us to include the effects of coactivation on the saccadic variability (see below, Model).
Model
Individual noise sources.
We mentioned three possible sources of saccade variability in the introduction: sensory signals, movement planning, and movement execution. Noise in planning and execution will result in noise in the motor commands that are sent to the muscles. Our approach allows us to estimate this motor command noise, but it does not allow us to identify the source of this noise upstream of the motoneurons. As a result, we are unable to determine unambiguously which part of the overall variability is caused by movement planning and which by movement execution. Accordingly, we will for now treat these two sources together, and refer to them as “motor noise.” We defer the decomposition of motor noise into planning noise and execution noise to Discussion.
In addition to sensory signals and motor noise, noise in the EyeLink system must also have contributed to the observed variability and should therefore be taken into account when modeling the observed variability. Hence, our model includes three different noise sources. We will term them “sensory noise,” “motor noise,” and “measurement noise,” and now discuss each of them.
The effect of measurement noise was estimated from the variability measured during stable fixation. Expressed as SDs, the precision was 0.04° horizontally and 0.03° vertically for all subjects. This gave rise to an uncertainty in the velocity estimate of 25°/s (SD). A simulation of the process of determining movement onset and offset from noisy EyeLink data recorded at 250 Hz and upsampled to 1000 Hz, suggested that, for a saccade of a certain duration, this process adds an uncertainty in the movement duration estimate of 1.6 ms (SD) and it produces a correlation coefficient of 0.45 between amplitude and duration.
Saccades cannot be executed more precisely than their target is localized. Localization by the peripheral retina is limited by several factors such as the size and density of photoreceptors and the precision with which retinal information is conveyed to and processed in central visual areas. Our term “sensory noise” refers to the combined effect of all of these factors. We estimated how precise the targets in our experiments were localized from a large number of psychophysical studies. Three factors determine this precision. First, localization of peripheral targets is more precise in the ςφ direction than in the ςR direction (Yap et al., 1987; White et al., 1992; Westheimer, 2005). Second, localization precision in the ςφ direction displays a strong “oblique effect”: precision is much better along the horizontal and vertical meridians than in oblique directions (White et al., 1992; Westheimer, 2001). Third, localization precision varies across the visual field in the same way as visual acuity does (Yap et al., 1987; Levi and Klein, 1990; Aitsebaomo and Bedell, 1992; White et al., 1992; Westheimer, 2001, 2005). We quantified the first two factors using the values reported by White et al. (1992): ςR/ςφ = 6 for targets on the horizontal and vertical meridians, and ςR/ςφ = 2 for targets in other directions. The way in which visual acuity varies across the visual field was estimated from the data of Wertheim (1894). Visual acuity decreases with eccentricity. When expressed as a CV (i.e., as an SD divided by the eccentricity), it is approximately constant for the larger eccentricities in our experiment, but it is larger for small eccentricities. A fit to Wertheim's data produced the following: CVacuity(e) = CVacuity,0(1.95exp(−e/2.65) + 1), where e is the eccentricity in degrees, and CVacuity,0 is the CV at large eccentricities. The data of Wertheim (1894) further suggest that CVacuity,0 depends on direction. It is smallest along the horizontal axis, it is ∼1.30 times larger upward, and ∼1.56 times larger downward.
We used the average of the results of White et al. (1992) and Aitsebaomo and Bedell (1992) as the CV of localization in the ςR direction at large eccentricities along the horizontal meridian: 0.045. To account for the directional dependence of the precision, we applied the following scaling to the localization precision CV in both the ςφ and ςR directions:
where CVhor represents the precision on the horizontal meridian, and angle θ represents the target direction (with θ = 0° for rightward, θ = 90° for upward, etc.). Factor V was 0.30 for targets above the horizontal axis and 0.56 for targets below that axis.
Saccades could have been planned to any position of the target disk (0.17° diameter). Assuming an equal probability for any position on this disk, this finite size led to an additional variance of
We will now describe the “motor noise” that we added to the estimated motor commands. Direct recordings of motor command signals in monkeys during saccades (Hu et al., 2007) and in cat (Gómez et al., 1986) and goldfish (Pastor et al., 1991) during fixation have shown that their firing exhibits “signal-dependent noise” (SDN). The properties of the motor unit pool of a muscle, such as recruitment, are such that the force produced by a muscle also has SDN (Jones et al., 2002; Hamilton et al., 2004). SDN is zero-mean, white Gaussian noise in the magnitude of the signal with an SD proportional to the magnitude of the signal, and it has been used extensively in modeling studies (Harris and Wolpert, 1998; Tanaka et al., 2004, 2006). We added SDN to the motor commands at each time-step of the mean trajectories. The SD ςSDN of this noise was defined as ςSDN2 = kSDN2u2, where u is the motor command and kSDN defines the level of the noise. SDN in the three muscle pairs was assumed to be independent.
As described above, we modeled muscle pairs as single ideal muscles that can produce positive and negative torques. This is a problem for modeling motor noise because antagonistic muscles are generally activated simultaneously (Miller and Robins, 1992). As a result of such coactivation, the torques generated by the two muscles will partially cancel, but the variances therein will add up. Consequently, the net torque will not have SDN, even when all motoneurons contributing to this torque individually do have SDN. Because no quantitative measures of coactivation are available, we included “constant noise” (CN) to account for its effect. A combination of SDN and CN can give an appropriate description of the noise in the torque. CN was modeled in the same way as SDN, but with an SD independent of the motor command: ςCN2 = kCN2. Constant and signal-dependent noise were assumed to be independent, giving rise to an SD of motor command noise of
Not all six extraocular muscles are equally strong. This is relevant because stronger muscles are less noisy than weaker muscles when both produce the same torque (Hamilton et al., 2004). The critical factor determining the noise level of a muscle is its number of motor units. The CV of the torque produced by a muscle is proportional to its number of motor units raised to the power of −1/2 (Hamilton et al., 2004). We therefore scaled the levels of SDN and CN according to this relationship for each muscle pair. Unfortunately, accurate estimates of the numbers of motor units for the human extraocular muscles are not available. However, two studies (Bors, 1926; Torre, 1953) suggest that the innervation ratio (the mean number of muscle fibers per motoneuron) is the same (∼5.5) for different extraocular muscles. Assuming the same innervation ratio for all muscles, we could estimate the numbers of motor units from the reported numbers of muscle fibers (Bors, 1926; Kato, 1938; Torre, 1953; Oh et al., 2001). After summing the numbers of muscles belonging to the same pair together, we arrived at the following ratios of motor unit numbers: horizontal:vertical:oblique = 1:0.77:0.50. We therefore scaled the levels of SDN and CN of the vertical and oblique muscle pairs by 0.77−0.5 = 1.14 and 0.50−0.5 = 1.41, respectively.
Finally, we added “temporal noise” (TN) to account for variability in movement duration. Because the results of experiment 1 show that the CV of movement duration does not vary with saccade amplitude (see Results), we defined the level of TN, kTN, as the CV of movement duration. Temporal noise was defined as a simultaneous scaling of the movement duration and the velocity profile: when duration was scaled by a factor 1 + α, velocity was scaled by a factor 1/(1 + α). This definition of temporal noise does not lead to variability in saccade endpoints; it produces variability in the temporal aspects of the saccade only. The same scaling was applied simultaneously to all muscle pairs. Temporal noise does not reflect noise that is added directly to the motor commands but it produces motor command variability in an indirect way. In the simulations, movement duration was varied by varying the time step in the calculation of a saccade trajectory from a motor command. All time steps within a single simulated movement were scaled by the same factor, but this factor was varied across simulated movements to obtain the desired CV of duration.
Combining the effects of all noise sources.
We combined the effects of all types of noise into a total variability that we could compare with the observed variability. The variability caused by the three types of motor noise was determined simultaneously in simulations in which all three noise types were present (see below, Fitting the model). We then added the effects of measurement and sensory noise.
The total endpoint variability was found by adding the covariance matrices attributable to motor, sensory, and measurement noise together (for each subject and target separately). Similarly, the SDs in amplitude, peak velocity, and movement duration were found by adding the variances attributable to all sources together. Sensory noise leads to variability in not only endpoints, but also in duration and peak velocity, because a movement that is planned to have a larger amplitude will, according to the main sequence, also have a longer duration and a larger peak velocity. To estimate these variabilities, we estimated each subject's main sequence relationships from the data of experiment 1 by fitting the observed duration and peak velocity to the function A(1 + exp(−amp/B)), where amp is the amplitude of the movement, and A and B are fit parameters. Local linearization of these functions then allowed us to transform SDs in amplitude into SDs in duration and peak velocity.
To determine predicted correlations between each pair of amplitude, duration, and peak velocity, we used the general equation for the correlation coefficient ρ(X,Y) for the correlation between X and Y as follows:
where Cov(X,Y) is the covariance between X and Y, and Var(X) is the variance in X. The preceding paragraph explained how we obtained the variances in the denominator. To determine the covariance in the numerator attributable to sensory noise, we again linearly approximated the main sequence relationships. This resulted in correlation coefficients of unity, so that each covariance equaled the product of two SDs. Measurement noise was assumed to produce a correlation coefficient of 0.45 between amplitude and duration (see above, Individual noise sources), and zero correlations in the other cases. This led to the following equations for the predicted correlation coefficients:
Here, the subscripts “mot,” “sens,” and “meas” refer to motor, sensory, and measurement noise, respectively. As for the experimental data, partial correlations were calculated from these simple correlations to account for the interrelationship between the three variables.
Fitting the model
We generated 500 sets of noisy motor commands for each target for each subject, and calculated the resulting movement trajectories. The endpoint of a simulated movement was defined as the position immediately after the last motor commands had been sent (other methods, such as using a velocity threshold produced very similar results). The variability in these movements was then combined with the variability attributable to sensory and measurement noise as described above.
The model had two free parameters: the levels of SDN and CN (kSDN and kCN). The level of TN was derived from the data of experiment 1, such that the model reproduced the observed mean CV of movement duration (i.e., kTN = 0.068). To estimate the free parameters, we fitted the model to the data from both experiments by optimizing the log likelihood of the observed endpoints (van Beers et al., 2004). This involved the following steps. First, for a given set of parameters, we simulated 500 movements for each target and each subject. We then determined the means and the covariance matrices of the simulated endpoints. This defined the model prediction. We next computed the log likelihood for each observed endpoint given the model prediction, and added these log likelihood scores together. This sum quantifies how well the observed endpoints match the predictions. The sum was optimized using an unconstrained nonlinear optimization algorithm (fminsearch in MatLab; The Mathworks, Natick, MA) that determined the two parameters that produced the best likelihood score. The same parameters were used for all subjects and targets.
Results
Experiment 1
Mean main sequence parameters
The top row of Figure 1 shows the mean amplitude, duration, and peak velocity as a function of the target amplitude. The first plot shows that the subjects on average made saccades of the required amplitude. The duration and peak velocity plots display the well known main sequence characteristics. Duration increases approximately linearly with amplitude, whereas peak velocity increases at a decreasing rate.
Main sequence parameters in experiment 1: mean and variability. A–F, Mean values for amplitude (A), duration (B), and peak velocity (C), and CVs of ςR and ςφ (D), duration (E), and peak velocity (F) as a function of target amplitude in experiment 1. In D, the CV of ςφ is the SD in the ςφ direction, divided by the mean amplitude. The black curves represent the mean of all subjects, and the gray areas represent 95% confidence intervals of the intersubject variability (1.96 times the SE).
Variability in main sequence parameters
All the endpoints, and their 95% confidence ellipses, of a representative subject are shown in Figure 2. The target is in some cases not located within this ellipse, which indicates that constant errors were made. Constant errors were also found for the other subjects, but their size and direction differed across subjects, and no general pattern emerged. We did not analyze the constant errors further. The confidence ellipses of the endpoints are approximately aligned with the movement direction. In other words, the variability in movement amplitude (ςR) is larger than the variability orthogonal to that (ςφ), in agreement with the results of van Opstal and van Gisbergen (1989). Endpoint variability increases with target amplitude.
Raw data in experiment 1. All endpoints, 95% confidence ellipses, and mean trajectories of a representative subject (HL) in experiment 1. The results for different targets have been offset vertically to avoid overlap.
The variability of the saccade endpoints, averaged over all subjects, is shown in Figure 1D. Two CVs are plotted; these are defined as ςR and ςφ divided by the mean amplitude. Both CVs decrease with amplitude, and both decrease in approximately the same way. This means that the ellipse shape does not change with amplitude, which was also found by van Opstal and van Gisbergen (1989). The CV of duration (Fig. 1E) is approximately constant, but the CV of peak velocity (Fig. 1F) decreases somewhat with amplitude. Interestingly, for most amplitudes, the CV in the amplitude is smaller than those in duration and peak velocity.
Correlations between main sequence parameters
The observation that the variability in amplitude tends to be smaller than that in duration and peak velocity can easily be explained by variations in movement speed. Movements of the same amplitude but with different durations will also have different peak velocities (assuming they have similar velocity profiles). This will result in a negative correlation between duration and peak velocity. Figure 3C shows that across different targets (denoted by different colors), peak velocity increases with duration. However, within sets of movements to the same target, we find the expected effect that peak velocity decreases with duration. This was also found by Jürgens et al. (1981) and Becker (1989). Figure 3F shows that the partial correlation coefficient (pcc) between duration and peak velocity, averaged over all subjects, is strongly negative for all amplitudes.
Correlations between main sequence parameters in experiment 1. A–C, Plots of raw data of duration as a function of amplitude (A), peak velocity as a function of amplitude (B), and peak velocity as a function of duration (C) for one representative subject (HL) in experiment 1. Confidence ellipses have been plotted for the sake of clarity; these have not been used in the analysis. D–F, Shown are pcc values between amplitude and duration (D), amplitude and peak velocity (E), and duration and peak velocity (F), as a function of target amplitude, averaged over all subjects. The format is the same as in Figure 1.
Figure 3A shows that across targets, movement duration increases with amplitude. This is one of the well known main sequence characteristics that was already shown in Figure 1B. Figure 3A also shows that duration and amplitude are also positive correlated within sets of saccades to the same target. Thus, saccades to a certain target that had a longer than average duration also tended to have a larger than average amplitude. This is confirmed by Figure 3D, which shows the pcc between duration and amplitude as a function of target amplitude, averaged over all subjects. The pcc decreases with amplitude but remains strongly positive for all amplitudes tested.
Across targets, peak velocity increases with amplitude (Fig. 3B). This is another part of the main sequence that was also shown in Figure 1C. Peak velocity increases with amplitude also within sets of saccades to the same target. The pcc between amplitude and peak velocity is strongly positive for all amplitudes, although its magnitude decreases with amplitude (Fig. 3E). Thus, amplitude, duration, and peak velocity are all correlated with one another.
Experiment 2
Mean main sequence parameters
The top row of Figure 4 shows plots of the mean amplitude, duration, and peak velocity as a function of the target direction. The amplitude generally agreed with the target amplitude of 9°. Saccade duration increased with the magnitude of the vertical component of the saccade, especially for downward saccades. Duration was minimal for purely horizontal saccades (mean duration, 48 ms), and it was significantly longer (p < 0.003) for purely upward (mean, 59 ms) and downward (mean, 61 ms) saccades. This pattern was found for all subjects and it agrees with previous studies (Oohira et al., 1983; Collewijn et al., 1988a,b; Kubo et al., 1991).
Main sequence parameters in experiment 2: mean and variability. A–F, Shown are plots of mean amplitude (A), duration (B), and peak velocity (C), and CVs of ςR and ςφ (D), duration (E), and peak velocity (F) as a function of target direction in experiment 2. The format is the same as in Figure 1.
Variability in main sequence parameters
All saccade endpoints and the mean trajectories of one representative subject are shown in Figure 5. As in experiment 1, confidence ellipses of the endpoints are approximately aligned with the movement direction. This figure also suggests that the ellipse shape may vary with the target direction. Relatively elongated ellipses are found for the horizontal targets and for the purely downward target. Ellipses in other directions tend to be rounder. Figure 4D confirms that, across all subjects, ςφ has minima for targets along or near the cardinal axes (one-tailed t test comparing the values for the four cardinal directions against the other 20 directions, p < 0.0003). There are no such minima for ςR; this SD varies smoothly with direction and increases with the magnitude of the vertical component of the saccade, especially for downward saccades. Peak velocity (Fig. 4C) is approximately the same for all directions.
Raw data in experiment 2. All endpoints, 95% confidence ellipses, and mean trajectories of a representative subject (JR) in experiment 2.
The variability of the duration and peak velocity is shown in Figure 4, E and F. As in experiment 1, the mean CV of movement amplitude was smaller (0.065) than those for duration and peak velocity (0.104 and 0.100, respectively).
Correlations between main sequence parameters
As in experiment 1, the pcc values between amplitude and duration and between amplitude and peak velocity are strongly positive, whereas that between duration and peak velocity is strongly negative (Fig. 6). None of these pcc values varies with direction and the values are similar to those for horizontal saccades of comparable amplitudes in experiment 1.
Correlations between main sequence parameters in experiment 2. A–C, Shown are partial correlation coefficients between amplitude and duration (A), amplitude and peak velocity (B), and duration and peak velocity (C) as a function of target direction in experiment 2, averaged over all subjects. The format is the same as in Figure 1.
Model predictions
Predictions of individual noise types
The aim of this study was to determine the sources of variability in saccades. We did this by determining which combination of measurement, sensory, and motor noise can explain the observed variability. We will first demonstrate the variability resulting from the individual noise sources. The contribution of measurement noise is very small and will not be shown.
The properties of the assumed sensory noise have been described in detail (see Materials and Methods, Model). The resulting variability is shown in Figure 7 in orange for what appeared to be very informative measures: the endpoint variability in experiments 1 and 2, and the three pcc values between amplitude, duration, and peak velocity in experiment 1. Endpoint variability is represented by 95% confidence ellipses centered on the mean endpoint location. Given sensory noise alone, the three pcc values between pairs of amplitude, duration, and peak velocity are unity. This is because a target that is localized further away will, according to the main sequence, provoke a saccade with a longer duration and a higher peak velocity.
Variability resulting from individual noise sources. A, B, Endpoint variability for subject HL in experiment 1 (A) and subject JR in experiment 2 (B), represented by 95% confidence levels. C–E, Shown are partial correlation coefficients between amplitude and duration (C), amplitude and peak velocity (D), and duration and peak velocity (E) in experiment 1. Shown are the following: observed data (black) and predictions by sensory noise only (orange), SDN only (red; kSDN = 0.172), CN only (green; kCN = 1.37 × 10−5 kg m2 s−2), and TN only (blue; kTN = 0.068). For each type of noise, the assumed noise level is the level in the model containing all types of noise. The levels are thus not optimal for the situations shown here.
The variability resulting from SDN in the motor commands is shown in Figure 7 in red (for kSDN = 0.172; not a best fit for this situation). This variability was calculated by setting the levels of all other noise types to zero. The predicted endpoint variability varies strongly with movement direction. Very elongated ellipses are produced along the cardinal axes. This can be understood from the extraocular muscles that are used. For a purely horizontal saccade, for instance, the horizontal muscles are activated strongly, resulting in noisy motor commands for these muscles. The vertical and oblique muscles are hardly activated and will therefore have little noise. This results in a large variability in movement extent, and virtually no orthogonal variability. In contrast, saccades in oblique directions are produced by substantial activation of both horizontal and vertical muscles. Both muscle pairs will thus have appreciable noise, which results in rounder confidence ellipses. The size of the confidence ellipses increases with movement amplitude because saccades with larger amplitudes require larger motor commands (that have more noise), and they take longer so that noise is added over a longer time. The pcc values involving duration are zero because SDN does not lead to variations in movement duration. The pcc between amplitude and peak velocity is strongly positive. This is because fluctuations of the motor command above its desired value will lead to both a higher peak velocity and a larger amplitude.
The variability resulting from CN in the motor commands is shown in Figure 7 in green (for kCN = 1.37 × 10−5 kg m2 s−2). The resulting endpoint variability hardly varies with movement direction. This is because CN does not depend on the motor commands. The variability is slightly larger in the vertical than in the horizontal direction because the vertical muscles are slightly noisier. The endpoint variability increases with movement amplitude because saccades with larger amplitudes have longer durations so that more noise is added. The three pcc values resulting from CN are quite similar to those caused by SDN. The correlations involving duration are zero because CN does not lead to variations in movement duration. The pcc between amplitude and peak velocity is positive but not as large as for SDN. This difference can be understood from the different natures of SDN and CN. For SDN, the fluctuations in the motor commands will be larger when the motor command is larger. Such large fluctuations will lead to stronger correlations between amplitude and peak velocity than the smaller fluctuations that result from CN.
The variability caused by TN is shown in Figure 7 in blue (for kTN = 0.068). The confidence ellipses are very small because temporal noise is defined in such a way that it does not produce variability in endpoints. The small variability shown is the result of inaccuracies in the numerical differentiation and integration in the simulations. The pcc between duration and peak velocity is strongly negative. This reflects how the temporal noise was defined: as variations of speed such that a saccade with a longer duration has a lower peak velocity. The two pcc values involving amplitude also take extreme values, but these are less important because they are caused by the above-mentioned inaccuracies in numerical differentiation and integration.
The full model
Figure 7 makes clear that none of the considered noise types in isolation can explain the observed variability. We therefore examined whether a combination of these noise types can do this. For that purpose, we optimized the log likelihood of the observed endpoints by finding the optimal values of the levels of SDN and CN, which were assumed to be the same for all subjects. The optimal values were as follows: kSDN = 0.172 and kCN = 1.37 × 10−5 kg m2 s−2.
The variability predicted by the model with these parameters is illustrated in Figure 8. The confidence ellipses produced by the model in Figure 8, A and B, are quite similar to those of the data, although they are somewhat larger. Overall levels of variability varied across subjects; the examples shown are from subjects with a smaller than average variability. For other subjects, the observed variability could be larger than that predicted. The other panels show that the model captures all aspects of the observed spatiotemporal variability very well: the 95% confidence areas of observations and predictions overlap almost everywhere. The main exception in experiment 1 is that the predicted magnitudes of the three pcc values for the smallest amplitudes are smaller than the actual ones. The most important result of experiment 2 (the direction dependence of ςR and ςφ, with minima for ςφ in the cardinal directions, and a ςφ that increases with the magnitude of the vertical component of the saccade) is predicted very well. Although the predicted CV in duration is slightly too small, especially for downward saccades, and the predicted magnitudes of the three pcc values are slightly too large, the overall results of experiment 2 are also well reproduced by the model. The underprediction of the CV in duration for downward saccades is related to the observation that this CV is somewhat larger for saccades in these directions, whereas we assumed it to be the same (kTN = 0.068) for all directions. It is unclear why the CV is larger for downward saccades. It is possible that duration is intrinsically more variable in these directions, but other factors, such as a reduced reliability of the eye tracker for downward gaze, could also play a role. We therefore decided not to model this effect, and kept the definition of temporal noise as simple as possible (i.e., as a constant CV in duration that holds for all directions).
Variability predicted by the full model. A, B, The endpoint variability predicted for subject HL in experiment 1 (A) and subject JR in experiment 2 (B). Results for experiment 1 occur in C–H: CVs of ςR and ςφ (C), CV of duration (D), CV of peak velocity (E), pcc between amplitude and duration (F), pcc between amplitude and peak velocity (G), pcc between duration and peak velocity (H). Results for experiment 2 occur in I–N: CVs of ςR and ςφ (I), CV of duration (J), CV of peak velocity (K), pcc between amplitude and duration (L), pcc between amplitude and peak velocity (M), and pcc between duration and peak velocity (N). The values predicted by the model are indicated by black confidence ellipses in A and B and by red curves in the other panels.
To find out how much each noise source contributes to the total variability, we calculated the endpoint variability for each noise source in isolation. In fact, the confidence ellipses in Figure 7 represent these individual contributions. Figure 9 shows the contributions as a percentage of the total endpoint variance. These contributions vary somewhat with movement amplitude and direction, but sensory noise is always the major source of endpoint variability. On average, sensory noise is responsible for 57% of the variance, whereas SDN and CN both contribute 21%. The contributions of TN and measurement noise are <1%.
Contributions of individual noise sources. A, B, Contributions to the total endpoint variance as a function of movement amplitude in experiment 1 (A) and movement direction in experiment 2 (B). “Total motor noise” is the sum of the contributions of SDN, CN, and TN.
Is it necessary to include all noise types to reproduce the observed variability? For some noise types, it is obvious that they have to be included. Sensory noise, for instance, has such a profound influence on saccade variability that not including this would be unrealistic. The importance of TN can be understood from the observed negative correlation between duration and peak velocity. Without TN, it is not possible to obtain a negative correlation because all the other noise types produce a positive or zero correlation (Fig. 7). The need to include CN is evident from the observed ςφ for saccades in cardinal directions. All other noise types lead to a ςφ that is much smaller than that observed (Fig. 7). CN is the only type of noise that generates appreciable values in these directions. The only noise source whose essentiality is not immediately obvious is SDN. We included SDN because physiological studies (Gómez et al., 1986; Pastor et al., 1991; Hu et al., 2007) suggest that oculomotor commands have SDN.
To examine whether SDN is an essential part of the model, we compared the predictions of the model just described to the predictions of an alternative model that includes all noise types except SDN. This alternative model was optimized in the same way as the full model. An objective way to compare how well these models predict the endpoint variability is to compare the likelihoods of the observed endpoints. However, a direct comparison of these likelihoods is not fair because the model with SDN has one free parameter more than the model without SDN. We therefore calculated Akaike's information criterion (AIC) (Akaike, 1973) for both models. AIC quantifies the goodness of fit of a statistical model based on its likelihood while taking into account the number of free parameters of the model. The lower the AIC, the better the model. For the full model, the AIC is −174.879, whereas it is −174.444 for the model with no SDN. Hence, including SDN significantly improved the quality of the fit, and the improvement is large (it amounts to the expected effect of >200 noninformative free parameters). This implies that all noise types included are essential.
Model validation
We validated the model by examining how well it predicts the variability in other data sets. The first data set consists of the trajectories of the left eye in experiments 1 and 2 for the saccades back from the target to the starting position. The variability in this data set was independent from that in the main data set because it consists of different saccades, made by different eyes. Importantly, the model was fitted to data from saccades that started from the primary position, whereas this data set consists of saccades that started from secondary and tertiary positions. A good prediction of the variability in these data would thus enlarge the generality of the model. We did not fit the model to this new data set but used the same parameter values as for the main data set.
The variability in these data was analyzed in the same way as that in the main data set. The results are shown in Figure 10, A–N. All aspects of the observed variability depend in approximately the same way on movement direction and amplitude as in the main data set. Thus, starting from a different position does not strongly influence the saccade variability. This is also what the model predicts (Fig. 10A–N). The agreement between observed and predicted variability is comparable with that of the main data set.
Model validation: observed and predicted variability in other data sets. A–N, Saccade variability for saccades of the left eye from the target back to the staring position. In B, the trajectories have been shifted such that the initial positions of all saccades coincide. Otherwise, the format is the same as in Figure 8. O–T, Saccade variability for saccades of the left eye from the starting position to the target in experiment 2.
We applied the model also to the left eye saccades from the starting position to the target in experiment 2, because most of these saccades moved this eye from a secondary to a tertiary position. Such saccades may be special because the eye orientation could tilt out of Listing's plane (Raphan, 1998). Figure 10, O–T, shows that the variability in these saccades is similar to that for saccades starting from the primary position (Fig. 8I–N), and this is also predicted by the model. We evaluate the effect of possible torsional blips in these data in the next section.
Control simulations
We made a large number of assumptions in the data analysis and in the formulation of the model. We have run control simulations to estimate the effects of possible inaccuracies of these assumptions.
Torsional variability
We assumed that Listing's law was obeyed during the recorded saccades. This allowed us to determine the full three-dimensional eye orientation from the recorded two-dimensional gaze directions. However, Listing's law is not obeyed perfectly, but the torsional component of eye orientation varies with an SD of 0.5–1.0° (Straumann et al., 1991; Desouza et al., 1997). Torsional variability can influence the predicted variability in the horizontal and vertical components because the motor commands for the horizontal and vertical muscles depend on the eye's torsion. We added variability to the torsion in control simulations. The same mean trajectories were used as before, but we added random torsion to each of the 500 simulated trajectories for each subject and target, while keeping the gaze direction the same. Random torsion values, drawn from zero-mean Gaussian distributions with a fixed SD, were added to the initial and final eye orientation of each simulated trajectory. The torsion at intermediate frames was assumed to follow the displacement profile of the gaze direction. We then fitted the model in the same way as the main model by optimizing the likelihood of the two-dimensional (horizontal–vertical) endpoints.
Table 1 shows the fit parameters for torsion SDs of 0.0 (the main model), 0.5, and 1.0°. The differences in the fit parameters are very small, at most several percent. The predictions for the quantities shown in Figure 8 completely overlap those of the main model and are therefore not shown. This suggests that the effects of assuming zero torsional variability are very small.
Optimal fit parameters for the main model and for simulations with variability added to the torsional component of the modeled saccades
Note that, as a direct result of motor noise, our model itself predicts nonzero torsional variability. How does this predicted variability compare with the reported variability? For the saccades in experiment 1, the predicted torsion SD increases with target amplitude. For an 18° saccade, the SD is ∼0.4°. This is only slightly smaller than the reported SD of 0.5–1.0° (Straumann et al., 1991; Desouza et al., 1997). Those values were in fact determined from experiments in which subjects produced many large saccades with amplitudes of up to 60°. Extrapolation of the predicted SD to such large amplitudes results in values within the reported range. This suggests that our model could well explain the reported torsional variability.
Torsional blips
So-called torsional blips (Tweed et al., 1994) form another possible violation of Listing's law. Torsional blips are torsional deviations from Listing's plane (the plane that contains all possible eye orientations obeying Listing's law) that can be induced by a saccade. Consistent torsional blips have been reported for horizontal saccades, with the peak amplitude occurring shortly before the end of the saccade (Straumann et al., 1995). The amplitude of these blips varied with saccade amplitude and eye elevation. Based on these data, we optimized the model again, this time assuming torsional blips with peak amplitudes of 10% of the horizontal saccade component, with intorsion occurring for abduction and extorsion for adduction. Table 2 shows that these blips influenced the fit parameters by several percent. However, the quality of the fit (as shown in Fig. 8) was not affected. Another simulation showed that torsional blips that depend on the vertical component of the saccade do not influence the fit parameters or the quality of the fit. Most importantly, the predictions for the data set that could be expected to have the largest torsional blips (the left eye saccades from the starting position to the target in experiment 2) (see above, Model validation), were not affected by the inclusion of torsional blips either. These results suggest that the estimated noise levels depend at most weakly on the presence of torsional blips, whereas the quality of the fit is not affected.
Optimal fit parameters for the main model and for a simulation with torsional blips
Orientation of Listing's plane
Listing's plane of the right eye was assumed to be orthogonal to the direction of gaze when this eye looked straight ahead. However, the actual orientation can be different. Listing's plane can have a nonzero pitch angle (i.e., being rotated about the interocular axis) and a nonzero yaw angle (i.e., being rotated about a vertical axis), where the angles are defined relative to the assumed orientation. Several studies (Ferman et al., 1987; Tweed and Vilis, 1990) have shown that when the head is in an upright position that is experienced as “natural,” the pitch angle is on average approximately zero. Pitch angles of individual subjects, however, can differ substantially from this mean, but the vast majority differs only several degrees (Haslwanter et al., 1994). The yaw angle is primarily determined by vergence (Mok et al., 1992; Bruno and van den Berg, 1997). When the vergence angle approximates zero (i.e., when the viewing distance is large), the yaw angle is close to zero. When the vergence increases (i.e., when the viewing distance decreases), Listing's plane rotates temporally at a rate that varies across individuals but averages at ∼0.72° for each degree of vergence per eye (Mikhael et al., 1995). For our viewing distance of 67 cm, the expected temporal rotation will be ∼2°.
We optimized the model again assuming different orientations of Listing's plane. We examined forward and backward rotations of 2.5°, and temporal rotations of 2 and 4°. Table 3 shows that the fit parameters hardly vary with orientation. Moreover, the predictions for the quantities shown in Figure 8 are in all cases virtually identical with those of the main model (data not shown). These results suggest that the estimated noise levels and the quality of the fit are as good as insensitive to possibly incorrect assumptions about the orientation of Listing's plane in our subjects.
Optimal fit parameters for the main model (top row) and for simulations in which different orientations of Listing's plane were assumed
Direction of muscle torques
We assumed that the directions of the torques produced by the extraocular muscles vary with eye orientation in the way proposed by Raphan (1998). It is, however, unclear how well this describes the actual directions during saccades. To estimate the possible effect of this assumption, we optimized an alternative model that was at the other extreme: muscle torque direction was assumed to be independent of eye orientation (as in Tweed and Vilis, 1987; Schnabolk and Raphan, 1994). Surprisingly, both fit parameters, and also the predictions for the quantities in Figure 8, were virtually identical with those of the main model (the largest difference occurred for kCN: 0.15%). This means that the matrix proposed by Raphan (1998) has virtually no effect on the predicted variability. Although surprising at first sight, this can easily be understood: this matrix mainly changes the torsional commands but those for the horizontal and vertical muscles are hardly affected.
Level of sensory noise
The level of sensory noise was not a free parameter because the uncertainty in target localization was taken from literature values (see Materials and Methods, Model). However, the actual localization uncertainty in our subjects could have been different. We repeated the optimization for different levels of sensory noise. This was done by scaling all the SDs and CVs that define the localization uncertainty by factors 0.9 and 1.1 relative to the baseline values, while keeping the anisotropy of localization uncertainty constant. In these simulations, the optimal level of TN (kTN) varied also because different levels of sensory noise lead to different amounts of variability in saccade duration. Table 4 shows that the level of sensory noise has a rather large effect on the fit parameters. All parameters decrease with increasing level of sensory noise. Both kCN and kTN do so weakly, whereas kSDN decreases more strongly. No notable changes occurred for the predictions shown in Figure 8. Hence, the level of sensory noise influences the fit parameters but not the quality of the fit.
Optimal fit parameters for different levels of sensory noise
Discussion
We measured the spatiotemporal variability in the trajectories of human saccades as a function of movement amplitude and direction. Our results confirm many previous findings of isolated aspects of this variability (Jürgens et al., 1981; Becker, 1989; van Opstal and van Gisbergen, 1989; Smeets and Hooge, 2003). One of our new observations is that the variability in movement direction is smaller for purely horizontal and vertical saccades than for saccades in oblique directions. We also found that saccade amplitude, duration, and peak velocity are all correlated with one another. All of the results together allowed us to determine the levels of different types of noise in a model that we used to identify the sources of saccade variability. This modeling suggests that uncertainty in localization of the target is the major source of variability in saccade endpoints, with noise in motor commands explaining a slightly smaller fraction. On top of this, there is temporal variability such that saccades with a longer than average duration have a smaller than average peak velocity, without affecting the saccade amplitude. The model has a large generality. It was fitted to a data set that contained only saccades starting from the primary position, but it predicts the variability in saccades starting from secondary and tertiary positions equally well.
Validity of the model
Although our model accounts for most aspects of the observed variability (Fig. 8), the predictions differed from the observations for some aspects. This occurred mainly for saccades with small amplitudes (2 and 4°), especially for the correlations between amplitude, duration, and peak velocity. This could well be related to the fact that these saccades had relatively short durations (22 and 31 ms on average) compared with the sampling rate of the eye tracker (one frame every 4 ms). This made the determination of the saccade onset and offset relatively inaccurate, and it may have led to underestimations of the variability in duration and overestimations of the magnitudes of the correlations in the data.
We performed control simulations to determine how the model predictions depend on a number of assumptions that we made regarding the torsion component of eye orientation during saccades. The simulations showed that these factors may have influenced the estimated levels of the various noise types weakly, but they did not affect the quality of the fit. The assumed direction of muscle torques, resulting from the effects of muscle pulleys and/or orbital fat, affected the results even less. This suggests that the effects of possible noise in the active control of muscle pulleys (Demer et al., 2000) will also be negligible. The only factor that influenced the fit parameters more strongly is the level of sensory noise. We used literature values about the precision with which visual targets are localized, but this precision can vary across subjects (Aitsebaomo and Bedell, 1992; White et al., 1992). Indeed, the observed saccade variability varied across our subjects (the raw data in Figs. 2 and 5 were from subjects with a lower than average variability). We initially tried to fit the model for each subject individually, with the level of sensory noise as an additional free parameter, but the individual data sets were too small to obtain reliable fits. We therefore assumed that for each noise type the noise levels were the same for all subjects. Although this will have led to suboptimal predictions for individual subjects, the noise levels determined this way can be considered as representative population estimates.
The sources of variability
Our model contains different types of noise: sensory noise and three types of motor noise: SDN, CN, and TN. Aitsebaomo and Bedell (1992) previously showed that both sensory noise and motor noise contribute to saccadic variability. Our results extend this by quantifying their contributions. Sensory noise explains ∼57% of the variance in saccade endpoints, whereas SDN and CN each explain ∼21%. It is not easy to interpret the absolute levels of SDN and CN that we estimated, because these depend on the time step used in the simulations (they scale with the time step raised to the power of −1/2). The level of TN, 0.068, is easier to interpret because it is the coefficient of variation in movement duration that remains after subtracting the variability in duration resulting from uncertainty in target localization. Thus, the intrinsic SD in duration is ∼7% of the mean duration. We modeled TN as a simple scaling of the duration and velocity of the saccade trajectory, such that the amplitude remained unaffected. In addition to our data, such temporal noise can also explain other phenomena such as the remarkably small variability in the product of peak velocity and duration (Smit et al., 1987), and compensation for variations in saccade speed (Quaia et al., 2000).
What are the sources of the various types of motor noise? SDN could have been expected to be the major type of motor noise because the noise in oculomotor neurons is approximately signal dependent (Gómez et al., 1986; Pastor et al., 1991; Hu et al., 2007). However, our results show that CN is equally important. As explained in Materials and Methods, coactivation of antagonistic muscles (Miller and Robins, 1992) is a likely explanation for CN. However, it is not possible to relate the level of CN to coactivation quantitatively, because no quantitative data on coactivation are available.
If the noise in individual muscles is indeed SDN, where does this come from? It is possible that noise in the magnitude of motor commands originates in the motoneurons themselves, but it seems more likely that at least a part of the noise arises upstream, in the brainstem regions in which the motor commands are generated. Noise in these areas will inevitably lead to variations in motor commands. Noisy variations in the representation of the saccade vector in the motor map of the superior colliculus, for instance, could lead to variability in saccade endpoints (van Opstal and van Gisbergen, 1989). Noise in motor commands could therefore be related to noise in the brainstem areas involved in saccade command generation.
However, the effects of noise in these structures are thought to be counteracted by a local feedback loop in the saccade command generating circuit (Robinson, 1975). This local feedback loop contains a neural integrator that mathematically integrates the activity emitted by the pulse generator. It compares the integrated activity to the activity that would be required for a saccade of the desired amplitude. The pulse is terminated when the difference between the two becomes zero. Consequently, the desired amplitude will be achieved regardless of fluctuations in the pulse activity. When the pulse happens to be higher or lower than desired, its duration is adjusted automatically. This most likely forms the basis of what we called temporal noise.
In Introduction, we mentioned three possible sources of saccade variability: sensory signals, movement planning, and movement execution. The preceding paragraph suggests that movement planning is the likely source of the temporal noise in our model.
A perfect feedback loop would almost completely counteract the consequences of noise in the pulse generator. Therefore, any noise in the magnitude of the motor commands that affects the endpoint of the saccade, should have arisen downstream of the feedback loop (i.e., in movement execution). However, it is possible that the feedback loop does not operate perfectly. In that case, noise in the magnitude of the motor commands could also originate in movement planning. Because it is not known how well the feedback loop works, we conclude that the SDN and CN in our model arise in movement execution and possibly also partly in movement planning, but we are unable to determine their exact individual contributions.
Conclusions
We conclude that the variability in human saccades is caused by a combination of uncertainty in target localization and noise in movement planning and execution. Uncertainty in target localization is the major source of endpoint variability, but the combined effect of movement planning and execution is only slightly smaller. This stresses the importance of carefully considering all possible sources when studying the variability in a motor system. It is unlikely that a single source can explain all variability, although the relative contributions may vary across motor systems and movement types.
Because, for saccades, the consequences of noise in movement planning and execution depend on the movement trajectory, our results suggest that the variability will depend on the desired trajectory. This is important because it suggests that, for a given saccade target, there may be an optimal trajectory for which the endpoint variance and therefore the expected movement error are minimal (Harris and Wolpert, 1998). Importantly, SDN in motor commands is a relatively small source of this variance. This is a serious departure from the commonly made assumption (Harris and Wolpert, 1998; Tanaka et al., 2004, 2006) that endpoint variance results from signal-dependent noise only. It remains to be determined which trajectories minimize the consequences of the actual noise.
Footnotes
-
This work was supported by The Netherlands Organization for Scientific Research Grant 451-02-013. I thank Jos van der Geest for the use of his setup and software, Steve Goldberg for pointing us to literature on motor unit numbers, and Jeroen Smeets, Eli Brenner, Maarten Frens, and Beerend Winkelman for fruitful discussions.
- Correspondence should be addressed to Robert J. van Beers, Department of Physics of Man, Helmholtz Institute, Utrecht University, Princetonplein 5, 3584 CC Utrecht, The Netherlands. r.j.vanbeers{at}phys.uu.nl
