**Abstract**:

**Abstract**
Optokinetic nystagmus (OKN) is a fundamental gaze-stabilizing response found in almost all vertebrates, in which eye movements attempt to compensate for the optic flow caused by self-motion. It is an alternating sequence of slow compensatory eye movements made in the direction of stimulus motion and fast eye movements made predominantly in the opposite direction. The timing and amplitude of these slow phases (SPs) and quick phases (QPs) are remarkably variable, and the cause of this variability is poorly understood. In this study principal components analysis was performed on OKN data to illustrate that the variability in correlation matrices across individuals and recording sessions reflected changes in the noise in the system while the linear relationships between variables remained predominantly the same. Three components were found that could explain the variance in OKN data, and only variables from within a single cycle contributed highly to any given component. A linear stochastic model of OKN was developed based on these results that describes OKN as a triple first order Markov process, with three sources of noise affecting SP velocity, the QP trigger, and QP amplitude. This model was used to predict the degree of signal dependent noise in the system, the duration of the transient state of SP velocity, and an apparent undershoot bias to the QP target location.

*accumulator models*and found that the best fitting distribution to QP interval histograms was given by the ratio of two correlated random variables (Waddington & Harris, 2012). We concluded that this was because SP duration was determined by the ratio of a variable SP amplitude threshold and SP velocity that varies between cycles.

*SD*= 5) years, participated in the study and had no neurological, ophthalmological, or vestibular impairments.

^{2}in either direction. Then, during a forward pass of the data, eye velocity was recorded and the peak velocity was determined at the time when velocity first began to decrease in magnitude and remained at a lower magnitude for 4 ms. The start and end point of the QP was then determined by both a (respectively) backward and forward pass of the data from the time of peak velocity to the time when eye velocity returned to a value between 0°/s and stimulus speed for 2 ms. This allowed us to collect QPs that were made in the direction of optic flow as well as QPs made in the opposite direction. All eye movements were reviewed in a customized interactive graphical interface. Blinks were detected manually and intervals containing blinks were marked and removed from analysis.

*m*SP-QP cycles, where

*m*ranged from 64 to 442. Six measurements were taken from each OKN cycle:

*x*,

_{i}*S*,

_{i}*V*,

_{i}*T*,

_{i}*y*and

_{i},*Q*, (

_{i}*i*= 1,

*m*) according to the scheme in Figure 1. Eye velocity changed between adjacent SPs, but could also vary during a SP (i.e., SPs could be nonlinear). We did not analyze the additional variance introduced by this intra-SP variability and instead measured SP velocity as the difference in eye position divided by the SP duration such that

*V*(

_{i}=*y*)

_{i}– x_{i}*/T*.

_{i}*n*-dimensional cloud of data points and rotating the line to minimize the squared orthogonal distance of each point to that line. This line passes through the maximum variation of the data and is termed the first principal component. Successive components are found by repeating this procedure to explain the remaining variance, with the added constraint that each component must be orthogonal to all preceding components. Mathematically, PCA is performed by eigenvalue decomposition, where the

*n*×

*n*input correlation (or covariance) matrix is diagonalized to yield

*n*eigenvalues and

*n*eigenvectors, where

*n*is the number of dimensions analysed in the data. In our investigation PCA was performed using the correlation (rather than covariance) matrices as it is based on standardized data, so the variables measured on different scales do not influence the results of the analysis disproportionately (Blunch, 2008).

*n*is an overestimate of the maximum degrees of freedom (DoF) because of the constraints imposed by the geometric relationships among the variables (Figure 1): Thus there is only a maximum of

*d*dimensions in the data, where

*d*is determined by the number of cycles of OKN and the number of variables in each cycle that are analyzed. We discarded those

*n*−

*d*components that do not explain any significant variance in the data. During our investigation we found that discarded components had eigenvalues that were either negligible (< 1 × 10

^{−5}) due to numerical round-off error or had values less than one (0.3 ± 0.2) due to the residual error from the linear approximation of Equation 3 implicit in PCA. This meant that the components were able to explain approximately 95% of the variability across all data sets. The error introduced by the linear approximation of Equation 3 could not be overcome by performing a reciprocal transformation of

*T*as this introduces a new nonlinear relationship between

_{i}*S*and 1/

_{i}*T*. However, PCA was performed using 1/

_{i}*T*as a variable in place of

_{i}*T*and we found qualitatively similar results.

_{i}*n – d*components we performed factor rotation using the varimax strategy to obtain orthogonal (uncorrelated) rotated components. Varimax rotation is a rotation of the component axes that maximizes the sum of the variances of the squared loadings of a component. This results in components that have high loadings for a few variables and low loadings for the other variables and aids in making the underlying structure of the loading patterns more clear. Oblique rotated components were also explored, as suggested by Costello and Osborne (2005) (using the promax strategy), but no qualitative differences were found in our results.

*d*components from each trial were sorted into categories according to their loading pattern. This was initially done by eye, but later numerical heuristics were developed to speed up the procedure. A principal component can be rotated such that it faces in the opposite direction but remains in the same dimension, causing all the loadings on a component to have the opposite sign. It was necessary to flip the sign of all loadings in these mirrored components so that they could be sorted in to the correct category. Each loading pattern was then displayed as a line plot of loading value against the original variable, and components placed in the same category were plotted on the same axis. This sorting was exhaustive and we found that all the loading patterns fell obviously into only three qualitatively different categories (see Results).

*x*, the partial autocorrelation of lag

_{t}*k*is the autocorrelation between

*x*and

_{t}*x*that is not already accounted for by lags one up to and including

_{t+k}*k*− 1. The sample partial auto-correlation function was calculated with the Forecasting options in SPSS 19.0 and the sample autocorrelation function was calculated using the built-in Matlab function

*corrcoeff*.

*x*at each lag as a function of

_{t}*x*

_{t}_{-1}plus some unexplained noise, whereas a first order moving average model explains the change in

*x*at each cycle as a function of the unexplained noise in the previous cycle. These identifying techniques were used in conjunction with the results of PCA to interpret the update dynamics of OKN variables.

_{t}*x*,

_{i}*S*,

_{i}*V*,

_{i}*T*,

_{i}*y*, and

_{i}*Q*were created. The model parameters were based on robust regression analysis (Matlab function

_{i}*robustfit*) of the relevant variables for each trial. These simulated sequences were then subjected to exactly the same correlation and PCA procedures as the original data.

*F*= 18.2,

*p*= 0.001) that saturated at higher stimulus speeds, producing some values of SP gain lower than 0.5 (Figure 2A). However, even when SP gain was close to unity we observed that SP velocity fluctuated over time to values much higher and lower than the mean SP velocity. We found that at stimulus speeds of 10°/s, 20°/s, 30°/s, and 40°/s retinal slip reached values greater than 5°/s in 6%, 48%, 79%, and 91% of SPs, respectively. This was unexpected as it indicated that even at moderate stimulus speeds a large proportion of SPs were being generated with retinal slip that should dramatically reduce visual acuity.

*F*= 10.1,

*p*< 0.001) with the mean asymptote at approximately 300 ms for speeds greater than 20°/s (Figure 2A). The distribution of SP duration was positively skewed and all but one were significantly different to Gaussian (Lilliefors test,

*p <*0.016 for 39 trials).

*F*= 13.5,

*p*< 0.001) and caused mean QP amplitude to become more negative (

*F*= 17.7,

*p*< 0.001) (Figure 2B). Mean SP amplitude and mean QP amplitude had a strong tendency to compensate for each other, so that over a trial eye position did not wander to the limit of gaze. However, the distributions of SP and QP amplitude during an experimental trial were often different.

*F*= 0.6,

*p*= 0.50).

*x*,

_{i}*S*,

_{i}*y*,

_{i}*Q*, and

_{i}*x*

_{i+}_{1}. Figure 4 illustrates, step by step, PCA performed on the correlation matrices generated from Participant 1 (left panel) and Participant 2 (right panel), viewing an OKN stimulus at 10°/s.

^{−6}, and 2.37 × 10

^{−6}. Due to the constraints imposed by Equations 1–3 there are at most three DoF in these data, so two components were discarded before factor rotation. The rotated loading patterns of the three retained components are illustrated as line plots of loading values against the OKN variable in Figure 4C, with each component placed in descending eigenvalue order from left to right. The largest eigenvalue component has negative loading from

*y*, positive loading from

_{i}*Q*, and negative loading from

_{i}*x*, thus representing QPs with a more negative start position having more positive amplitude and a more negative end position. The second component represents QPs with more negative amplitude having a more negative end position, and the third component represents SPs with a more negative start position having more positive amplitude.

_{i}^{−6}, and 3.63 × 10

^{−7}, and the rotated loading patterns of the three retained components are illustrated in Figure 4F. Despite the differences in the correlation matrices we observed similar loading patterns in the components but expressed in different eigenvalue order.

*Q*-component, the

*S*-component, and the

*V*-component. These noise processes extended in time predominantly only across adjacent cycles, implying a first order Markov process.

*noise*(corresponding to the

*S*-,

*Q*-, and

*V*-components) were normally distributed. Second, we assumed parsimoniously that the process could be explained by a first order Markov process. That is, the state of the system is conditional only on the immediate previous cycle and not explicitly on earlier cycles. Third, we assumed the model parameters associated with the eigenvectors (

*a*,

*b*,

*c*, and

*d*; see below) would be constant across trials. Fourth, we assumed the model parameters associated with the eigenvalues (

*e*,

*ε*(

_{s}*i*),

*ε*(

_{q}*i*), and

*ε*(

_{v}*i*); see below) would be free parameters and allowed to vary between trials.

*S*-component there was a high positive loading from

*S*, positive loading from

_{i}*V*, and negative loading from

_{i}*x*(Figure 5). This loading pattern indicated that SP amplitude becomes more positive (increases in magnitude) with more positive SP velocity and a more negative SP start position. This correlation could also be observed in regression of

_{i}*S*against

_{i}*x*and

_{i}*V*with a mean

_{i}*R*

^{2}= 0.25 (

*p*≤ 6 × 10

^{−5}for all trials). We express these relationships in the form: where

*a*and

*b*are constants and

*ε*(

_{s}*i*) is a normal random process with mean

*ŝ*and standard deviation

*σ*. We assume that the values of

_{s}*ŝ*and

*σ*remain constant during each trial. In the

_{s}*S*-component there were also positive loadings from

*T*and

_{i}*y*, sometimes high and sometimes low, but these loadings appear to be forced by the constraints in Equation 3 and Equation 2, respectively.

_{i}*Q*-component with a positive loading from

*Q*and a negative loading from

_{i}*y*(Figure 5) indicating that QP amplitude becomes more positive (decreases in magnitude) with a more negative QP start position. Evidence of

_{i}*V*loading on the

_{i}*Q*-component was not as clear as in the S-component, as it was sometimes positive and sometimes negative and occurred in association with loading from

*S*. However, regression of

_{i}*Q*against

_{i}*V*and

_{i}*y*revealed a linear trend with a mean

_{i}*R*

^{2}= 0.34 (

*p*≤ 2.2 × 10

^{−5}for all trials). We express these relationships in the form where

*c*and

*d*are constants and

*ε*(

_{q}*i*) is a normal random process with mean

*qˆ*and standard deviation

*σ*. We assume that the values of

_{q}*qˆ*and

*σ*remain constant during each trial. We interpret the other loadings in the

_{q}*Q*-component to be forced by the constraints expressed in Equations 1–3.

*V*-component of each cycle there was a high positive loading from

*V*and corresponding negative loading from

_{i}*T*with a small or absent loading from

_{i}*S*(Figure 5). This indicates that the

_{i}*V*-component describes whether SP velocity is fast or slow independent of its amplitude. This loading pattern implies that the

*V*-component is self-contained and does not depend directly on the other variables. However, from Figure 3A we see that

*V*is autocorrelated and the regression of

_{i}*V*

_{i}_{+1}against

*V*reveals a clear linear relationship with a mean

_{i}*R*

^{2}= 0.28 (

*p*≤ 3.3 × 10

^{−4}for 35 of 40 trials). As the principal components reflect uncorrelated noise processes, we express this relationship in the form of an autoregressive rather than a moving average model where

*e*is a free parameter, and

*ε*(

_{v}*i*) is a normal random process with mean

*vˆ*and standard deviation

*σ*. We assume that the values of

_{v}*e*,

*vˆ*, and

*σ*remain constant during each trial.

_{v}*σ*,

_{s}*σ*, and

_{q}*σ*) were calculated from the residuals of regression.

_{v}*a*,

*b*,

*c*, and

*d*were −0.250, 0.158, −0.478, and −0.166, respectively, and were considered constants for the purposes of our analysis. The mean values of the noise processes (

*ŝ*,

*qˆ*, and

*vˆ*) were then recalculated for each trial, keeping the model parameters

*a*,

*b*,

*c*, and

*d*fixed.

*a*,

*b*,

*c*, and

*d*), and the free parameters (

*e*,

*ε*,

_{s}*ε*, and

_{q}*ε*). The simulated data sets were created iteratively using Equations 1–6 with the additional constraints where

_{v}*z*is a small positive value to mimic the QP dead zone, which we set to 1°.

*e*or by

*vˆ*(Equation 6). A simplified closed-loop linear control system that could generate OKN is illustrated in Figure 7A, and the update dynamics of SP velocity in this system are given by where

*g*

_{1},

*g*

_{2}, and

*g*

_{3}are the gain elements of retinal slip, efference copy, and the forward gain of SP velocity, respectively.

*e*and

*vˆ*relate to the gain elements in the proposed control system: and Thus the model predicted that

*e*and

*vˆ*must be related, as they share the term

*g*

_{3}

*g*

_{1}. Rearranging Equation 13 and substituting into Equation 12 illustrates that: and this relationship can be observed in a scatter plot of the values of

*e*against −

*vˆ*/

*V*from each trial (Figure 7B;

_{S}*R*

^{2}= 0.72). Orthogonal regression performed on these data gives a slope of 1.02 and a

*y*-intercept of 0.87 that indicated the approximate value of the product of the gain elements

*g*

_{3}and

*g*

_{2}.

*V*-component may be a consequence of either

*g*

_{1}or

*g*

_{3}(or both) varying, as these gain elements are effectively multiplied by the stimulus speed to arrive at the noise process

*ε*.

_{v}*S*- and

*Q*-components we compared the values of σ

*and σ*

_{s}*against the mean magnitude of SPs and QPs for each trial. Orthogonal regression illustrated a linear relationship between the standard deviation of the*

_{q}*S*-component and the mean SP magnitude (Figure 8A;

*R*

^{2}= 0.75), and a linear relationship between the standard deviation of the

*Q*-component and the mean QP magnitude (Figure 8B;

*R*

^{2}= 0.50). These relationships can be expressed in the form and where

*k*= 0.36,

_{s}*k*= 0.38, and

_{q}*u*and

_{s}*u*were not significantly different to zero. It is quite well-known that there is proportional noise in saccades to stationary targets, but the constant of proportionality is usually much lower at ∼ 0.1 (van Beers, 2008). It appears that in OKN the constant of proportionality is not only much larger but is also observed in the amplitude of SPs.

_{q}*V*-component, regression of

*σ*was performed against the values of stimulus speed (

_{v}*V*), mean SP velocity (

_{S}*V¯*), and mean retinal slip. There was a weak relationship between

*σ*and

_{v}*V¯*(

*R*

^{2}= 0.29), a moderate relationship between

*σ*and retinal slip (

_{v}*R*

^{2}= 0.46), and a strong relationship between

*σ*and

_{v}*V*(

_{S}*R*

^{2}= 0.69). We express the relationship between the standard deviation of the

*V*-component and the stimulus speed (Figure 8C;

*R*

^{2}= 0.69) in the form: where

*k*= 0.12 and

_{v}*u*was not significantly different to zero. This agrees with our hypothesis that a variable gain element in the OKN control system would cause noise in SP velocity that is proportional to the stimulus speed.

_{v}*x*

_{i}_{+1},

*y*

_{i}_{+1},

*S*

_{i}_{+1}and

*Q*

_{i}_{+1}, given by: Here we can see that each of these variables depends on its value only in the previous cycle and also the SP velocity in the previous cycle. In other words, the dynamics are first order Markov and driven by SP velocity, which is itself another first order Markov process. The variables

*x*

_{i}_{+1},

*y*

_{i}_{+1},

*S*

_{i}_{+1}, and

*Q*

_{i}_{+1}are all correlated with each other through a common dependence on SP velocity and the three noise sources

*ε*(

_{s}*i*),

*ε*(

_{q}*i*), and

*ε*(

_{v}*i*).

*x*

_{i+}_{1},

*y*

_{i+}_{1},

*S*

_{i+}_{1}, and

*Q*

_{i+}_{1}all have the same transient rise time determined by the quantity (1 +

*a*)(1 +

*c*). Empirically, this value is approximately 0.39 indicating that the process is stable and will eventually reach steady state. For SP velocity the transient rise time depends on

*e*and with a range of values from 0.08 to 0.90 was also always stable.

*vˆ*and

*e*calculated from each trial. To reach 63% of the steady state SP velocity from an initial speed of 0°/s took on average 2.3 cycles (

*SD*= 2.0). Typical values of SP duration would equate this transient response to a period less than 1 s. Thus, steady state is reached quickly, and this result is in good agreement with the study by Abadi, Howard, Ohmi, and Lee (2005) who showed that steady state SP velocity was reached in approximately two OKN cycles. We also simulated the update dynamics of the SP start position and found from an initial position of 0° that the start position took an average 2.4 cycles (

*SD*= 1.6) to reach 63% of the steady state value, demonstrating that the transient response for contraversion was also extremely fast.

*V¯*→ ∞.

*Q*= 0 to give: where

_{i}*y*is the desired target location. Estimating

_{T}*V*with

_{i}*V¯*and

*ε*(

_{q}*i*) with

*qˆ*, we found the average

*y*for each trial. Mean values of

_{T}*y*were −6.8°, −8.3°, −8.9°, and −10.2° for trials with a stimulus speed of 10°/s, 20°/s, 30°/s, and 40°/s, respectively, illustrating a main effect of stimulus speed (

_{T}*F*= 5.9,

*p*= 0.003). These results indicated that QPs were being targeted in to the negative field and not straight ahead. Calculating values of

*y¯*and

*x¯*from Equations 19 and 20 we found the average amplitude of QPs and the average distance to the target location. We found that QPs undershoot this target location with mean values of −2.1°, −2.3°, −3.1°, and −4.3° for trials with a stimulus speed of 10°/s, 20°/s, 30°/s, and 40°/s, respectively, demonstrating a main effect of stimulus speed (

*F*= 15.9,

*p*< 0.001). We considered this undershoot bias as a value of saccadic gain (QP amplitude ÷ QP target amplitude) and found it was approximately 0.64, with no main effect of stimulus speed (

*F*= 2.1,

*p*= 0.122). These results implied that QPs have a tendency to undershoot their target location, possibly due to adaptive changes in the value of saccadic gain or as a result of not fully compensating for the movement away from the target location during the SP.

*S*-,

*Q*-, and

*V*-components (Equations 4–6). Although mathematically simple, the behavior of the three interacting noise processes is non-trivial.

*V*-component is related to stimulus speed, in accordance with Kolarik, Margrain, and Freeman (2009) who also found that variation in velocity between SPs increased linearly with stimulus speed. Stimulus speed, retinal slip, and mean SP velocity are dependent on each other (retinal slip = stimulus speed – eye velocity). However, only retinal slip and eye velocity by way of efference copy (Sperry, 1950) are readily available to the system, whereas stimulus speed must be reconstructed internally. The error term in our

*V*-component then appears to reflect an internally generated noise. We propose that the reason for this signal dependent noise is a variable element in the forward loop of the OKN control pathway (either the overall reflex gain or the neural integrator), which is effectively multiplied by stimulus speed.

*Q*-component is linearly related to the mean QP magnitude. This result is not entirely unexpected as it is well-known that larger saccades have larger errors to static visual targets. What is surprising is that the constant of proportionality is approximately 0.38, which is much more than typical static saccades. The reason for this is unknown but may reflect targeting error when saccades are made to a position without any obvious visual target or simply error when making saccades to a moving target. The standard deviation of the

*S*-component is also linearly related to the mean SP magnitude, indicating the presence of proportional noise in determining the threshold at which to trigger a QP and was unexpected.

*Q*- and

*S*-components control mean eye position in a rather simple way. QP amplitude is negatively correlated with the start position of the QP (or end position of the previous SP) and will therefore tend to correct for the end position of the previous SP. Similarly, SP amplitude is negatively correlated with SP start position (or the end of the previous QP) and will tend to correct for end position of the previous QP. The

*Q*- and

*S*-components cooperate to maintain mean eye position, via rapid start position feedback.

*b*estimated from our data was 0.16, indicating a refractory period of 160 ms. However, this conclusion should be considered with caution, as SPs clearly can occur with durations under this value. The QP amplitude is also dependent on both the QP start position and SP velocity. Similar relationships have been observed in the threshold for triggering the QP of vestibular nystagmus in humans (Lau, Honrubia, & Baloh, 1978) and rabbits (Lau & Honrubia, 1986). This implies that these relationships have biological relevance in other species and even other forms of nystagmus, and future studies might investigate how the model parameters change between foveate and afoveate species or between optokinetic and vestibular nystagmus.

*T*=

_{i}*S*, and assuming normality of the error terms this variable would be given by a positively skewed distribution. In a previous study we investigated the goodness of fit of six different probability density functions to the histograms of SP duration and found that the distribution of two correlated and left truncated (at zero) normal variables gave a significantly better fit to the data than the other models tested (Waddington & Harris, 2012). We propose that this is due to the approximately linear trajectory of SPs and that QPIDs are given by the ratio of a variable SP amplitude threshold and SP velocity that varies between cycles.

_{i}/ V_{i}*Perception*

*,*34, 1015–1024. [PubMed]

*Biological Cybernetics*

*,*75, 1–9. [CrossRef]

*Chaos, Solitons and Fractals*

*,*9, 343–361. [CrossRef]

*Biological Cybernetics*

*,*66, 203–216. [CrossRef] [PubMed]

*Journal of Neurophysiology*

*,*93, 1510–1522. [CrossRef] [PubMed]

*Introduction to structural equation modelling using SPSS and AMOS*. Thousand Oaks, CA: Sage Publications, Ltd.

*Time series analysis: Forecasting and control*. San Francisco, CA: Holden-Day.

*Vision Research*

*,*12, 293–304. [CrossRef] [PubMed]

*Vision Research*

*,*22, 479–484. [CrossRef] [PubMed]

*Ophthalmic Research*

*,*25, 91–93. [CrossRef] [PubMed]

*Contemporary Ocular Motor and Vestibular Research*(pp. 185–187). Stuttgart, Germany: Georg Thieme.

*Vision Research*

*,*14, 1053–1058. [CrossRef] [PubMed]

*Biological Cybernetics*

*,*28, 209–221. [CrossRef] [PubMed]

*Journal of Physiology*

*,*270, 321–344. [CrossRef] [PubMed]

*Annals of the New York Academy of Sciences*

*,*374, 421–433. [CrossRef] [PubMed]

*Signal Processing*, 36, 287–314. [CrossRef]

*Practical Assessment Research & Evaluation*

*,*10, 1–9. Retrieved from http://pareonline.net/getvn.asp?v=10&n=7.

*Journal of Neurophysiology*

*,*103, 543–556. [CrossRef] [PubMed]

*Experimental Brain Research*

*,*81, 46–52. [CrossRef] [PubMed]

*British Journal of Ophthalmology*

*,*85, 1477–1483. [CrossRef] [PubMed]

*Annals of the New York Academy of Sciences*

*,*956, 445–448. [CrossRef] [PubMed]

*Investigative Ophthalmology and Visual Science*

*,*44, 3833–3841, http://www.iovs.org/content/44/9/3833. [PubMed] [Article] [CrossRef] [PubMed]

*Visual Neuroscience*

*,*5, 107–122. [CrossRef] [PubMed]

*Nonlinear Dynamics*

*,*44, 367–380. [CrossRef]

*Annals of the New York Academy of Science*

*,*545, 180–186. [CrossRef]

*Acta Oto-Laryngologica*

*,*65, 441–448. [CrossRef] [PubMed]

*Frontiers in Bioscience*

*,*13, 1899–1916. [CrossRef] [PubMed]

*European Journal of Neuroscience*

*,*8, 92–105. [CrossRef] [PubMed]

*Journal of Vision*, 9(13):1, 1–11, http://www.journalofvision.org/content/9/13/1, doi:10.1167/9.13.1. [PubMed] [Article] [CrossRef] [PubMed]

*Markov processes for stochastic modeling*. London: Chapman & Hall.

*Experimental Brain Research*

*,*201, 271–282. [CrossRef] [PubMed]

*Journal of Comparative Physiology A*

*,*185, 341–352. [CrossRef]

*Vestibular mechanisms in health and disease*(pp. 37–44). New York: Academic.

*Journal of Comparative Physiology A*

*,*160, 585–592. [CrossRef]

*Journal of Neurophysiology*

*,*10, 153–158.

*Eye movement research: Mechanisms, processes and applications*(pp. 47–62). Amsterdam: Elsevier.

*European Journal of Neuroscience*

*,*10, 811–822. [CrossRef] [PubMed]

*Handbook of physiology. The nervous system. Motor control*(pp. 1275–1320). Bethesda: American Physiological Society.

*IEEE Transactions on Biomedical Engineering*

*,*39, 1319–1321. [CrossRef] [PubMed]

*Biological Cybernetics*

*,*76, 237–250. [CrossRef]

*Journal of Comparative and Physiological Psychology*

*,*43, 482–489. [CrossRef] [PubMed]

*Neuroreport*

*,*11, 1987–1990. [CrossRef] [PubMed]

*Biological Cybernetics*

*,*87, 68–78. [CrossRef] [PubMed]

*PLoS ONE*3

*,*e2070. doi: 10.1371/journal.pone.0002070.

*Experimental Brain Research*. Advanced online publication. doi: 10.1007/s00221-012-3297-z.

*Vision Research*

*,*2, 69–80. [CrossRef]

*Acta Oto-Laryngologica*

*,*511, 104–108. [CrossRef]