Free
Article  |   November 2012
Human optokinetic nystagmus: A stochastic analysis
Author Affiliations
  • Jonathan Waddington
    Centre for Robotics and Neural Systems and Plymouth Cognition Institute, Plymouth University, Plymouth, UK
    jonathan.waddington@plymouth.ac.uk
  • Christopher M. Harris
    Centre for Robotics and Neural Systems and Plymouth Cognition Institute, Plymouth University, Plymouth, UK
    C.M.Harris@plymouth.ac.uk
Journal of Vision November 2012, Vol.12, 5. doi:10.1167/12.12.5
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Jonathan Waddington, Christopher M. Harris; Human optokinetic nystagmus: A stochastic analysis. Journal of Vision 2012;12(12):5. doi: 10.1167/12.12.5.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract
Abstract
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.

Introduction
Optokinetic nystagmus (OKN) is a fundamental oculomotor response elicited in nature by the image motion across the retina as an animal moves through the environment. It is essential for maintaining optimal visual acuity and can be observed in almost all vertebrates (Huang & Neuhauss, 2008; Walls, 1962) and some invertebrates (Land, 1999) with mobile head or eyes. Although archaic, OKN differs among species depending on the repertoire of oculomotor subsystems and patterns of optic flow. In humans the response is dominated by OKN with fast dynamics, termed early OKN and the ocular following response (Cohen, Henn, Raphan, & Dennett, 1981; Cohen, Matsuo, & Raphan, 1977; Gellman, Carl, & Miles, 1990; Miles, 1995, 1998; Miles, Kawano, & Optican, 1986). 
Human OKN can be elicited by instructing a subject to maintain gaze on a region in space whilst viewing a large moving visual scene (Honrubia, Downey, Mitchell, & Ward, 1968), either around the subject (rotational OKN) or in the fronto-parallel plane (translational OKN). The response is an alternating sequence of slow phases (SPs) and quick phases (QPs). QPs are fast movements with velocity profiles similar to saccades (Garbutt et al., 2003; Garbutt, Harwood, & Harris, 2001; Kaminiarz, Königs, & Bremmer, 2009). They are usually made in the opposite direction to optic flow and tend to direct the eyes into a more eccentric position. SPs are slow movements made in the direction of stimulus motion with a gain (SP speed ÷ stimulus speed) less than unity, which decreases with stimulus speed (Fletcher, Hain, & Zee, 1990) so complete retinal stabilization is seldom achieved. SPs tend to bring the eye position to a more central location, but the timing and amplitude of both SPs and QPs are highly variable (Anastasio, 1996; Balaban & Ariel, 1992; Carpenter, 1993, 1994; Cheng & Outerbridge, 1974; Trillenberg, Zee, & Shelhamer, 2002). This variability is intrinsic and implies either an embedded stochastic process or complicated deterministic behavior manifesting as chaos. Low fractional correlation dimensions have been reported, implying chaotic behavior (Shelhamer, 1992, 1996). However, stochasticity affects the correlation dimension (Argyris, Andreadis, Pavlos, & Athanasiou, 1998), and purely stochastic models of nystagmus can also demonstrate low fractional dimensions (Harris & Berry, 2006). 
A complete model of the system has remained elusive due to the essentially distinct neural pathways that generate QPs and SPs and the different approaches used to investigate the two subsystems. The neural circuitry that generates SPs is relatively well understood in terms of the underlying control system (Robinson, 1981) and those structures in the brain that are known to play an important role such as the pretectum complex (Hoffmann, 1988; Ilg & Hoffmann, 1996). Although much is known about the neural circuitry involved in generating saccades, the mechanism that determines when a QP is triggered and its target position is not well understood. However, the timing of QPs is often modeled as a stochastic signal accumulating information (or integrating neural activity) to a threshold, at which time a QP is generated (Anastasio, 1996; Balaban & Ariel, 1992; Carpenter, 1993, 1994; Trillenberg et al., 2002). Recently, we compared six of these 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. 
Despite a variable gaze position from cycle to cycle, average position across cycles tends to remain steady. This indicates some adaptive mechanism tends to keep gaze in the desired location and implies that SP and QP variables must be correlated. Correlations between the amplitude and velocity of SPs (Watanabe, Ohmura, Shojaku, & Mizukoshi, 1994) and between the start position of QPs and QP amplitude (Balaban & Ariel, 1992) have been observed. However, a stochastic model of the system has been difficult to ascertain due to variability in the strength of correlations across individuals and experimental trials. 
In this study principal components analysis (PCA) was performed on eye movement data recorded during a sustained period of OKN stimulation to investigate the underlying structure of the correlation matrices of OKN variables. Each component extracted with PCA was made up of a weighted linear sum of the analyzed variables (the eigenvector) and explained a certain percentage of the variance in the data (proportional to the eigenvalue). The results of our analysis indicated that individual differences in the correlation matrices were due to changes in the eigenvalues of components while the eigenvectors remained predominantly the same. This implied that the amount of noise in the system was changing between recording sessions while the linear relationships between variables remained relatively constant. 
The similarity of the components between participants and trials allowed us to develop a reliable stochastic model of OKN eye movements using only simple linear stochastic equations, with dependencies among adjacent cycles and three sources of noise affecting SP velocity, QP triggering, and QP amplitude. The most surprising result is that as SP velocity wanders in a Markov fashion above and below average SP velocity (in our data and in the model). In many cases retinal slip reaches values above 5°/s which should dramatically reduce visual acuity. 
Materials and methods
Ethics statement
All protocols were approved by The University of Plymouth Faculty of Science Human Research Ethics Committee. Participants gave informed written consent and were made aware of their right to withdraw at any time. 
Participants
Ten healthy adults (six female and four male), mean age 25 (SD = 5) years, participated in the study and had no neurological, ophthalmological, or vestibular impairments. 
Procedure
Participants sat in a chair 1 m from the middle of a flat white screen that covered as much of the field of view as reasonable to approximate a hemispherical display (viewable area of 1.57 m by 1.17 m landscape, subtending 76° by 61°). The OKN stimulus was rear projected (EPSON EMP-500; Seiko Epson Corp., Japan) onto the screen at a frame rate of 60 Hz. The room was kept completely dark except for the projection onto the screen during stimulation. The participant's head was constrained using a chin rest. Eye movements were measured using a binocular head-mounted limbus tracker (Skalar IRIS Infrared Light Eye Tracker; Skalar Medical BV, Netherlands) that recorded horizontal movement with a resolution of 3 minarc at a sampling rate of 1 kHz. The eye tracker was calibrated for each participant at the start of each procedure by recording the voltage output during fixation of 40 targets at different positions on the screen. The voltage was linearly proportional to eye position within ±20° of the center of the display and linear regression was used to generate a calibration scale and offset. 
Eye movement measurements were recorded on computer (vsgEyetrace v.3.0.beta software for Windows; Cambridge Research Systems, UK) and stored on the hard disk of a computer for offline analysis. 
Stimulus
Translational OKN was elicited with a flat vertical square wave grating, comprised of alternating black and white vertical stripes moving horizontally at a fixed tangential speed. Each recording session was comprised of a pseudorandom sequence of four trials, each with a different stimulus speed equivalent to 10°/s, 20°/s, 30°/s, or 40°/s at the center of the screen. We note that translational optic flow generates different rotational stimulus speeds depending on where the subject looks on the screen. As an example, in our procedure the participant sits at a distance of 1 m and the rotational stimulus speed appears to be approximately 6% slower when the participant looks 20° away from center. Each stimulus presentation was recorded over 160 s so that we could collect a long time series of OKN data to reliably perform the proposed analyses. As participants often find it uncomfortable to view very high spatial frequency stimuli for long periods we chose to use a stimulus with a spatial frequency of 0.1 cycles/°, which generated a reliable OKN response but minimized the likelihood that participants would purposefully look away from (or through) the screen to alleviate discomfort. Participants were asked to stare at the center of the screen rather than follow the moving lines to evoke ”stare” OKN rather than “look” OKN, and attention was maintained by giving brief verbal feedback about the amount of time lapsed during the stimulus presentation approximately every 10 s. Participants were also given a break for one minute between each trial to alleviate discomfort, tiredness, and to minimize the effects of any optokinetic after-nystagmus. 
Data analysis
Movement in the direction of the stimulus motion was defined as positive and movement in the opposite direction as negative. Position was defined relative to the center of the stimulus display and is positive on the side that the stimulus is drifting towards and negative on the opposite side. 
All programs and algorithms for analyzing data were developed and created in MATLAB (MATLAB; Mathworks, USA). Each eye was calibrated separately, and the average was computed to yield a cyclopean eye position. Eye velocity was derived from the eye position using a central difference algorithm and zero-phase digital filtered with an 80 Hz Butterworth filter. Eye acceleration was derived from the filtered eye velocity data using a central difference algorithm. 
Possible QPs were detected when eye acceleration was greater than 1000°/s2 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. 
After blinks were extracted, each recorded trial contained m SP-QP cycles, where m ranged from 64 to 442. Six measurements were taken from each OKN cycle: xi, Si, Vi, Ti, yi, and Qi, (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 Vi = (yi – xi)/Ti
Figure 1
 
Geometric representation of OKN and the characteristic variables. The ith cycle in a series contains one SP followed by one QP and is defined by six variables: SP start position xi, SP amplitude Si, average SP velocity Vi, SP duration Ti, QP start position yi, and QP amplitude Qi.
Figure 1
 
Geometric representation of OKN and the characteristic variables. The ith cycle in a series contains one SP followed by one QP and is defined by six variables: SP start position xi, SP amplitude Si, average SP velocity Vi, SP duration Ti, QP start position yi, and QP amplitude Qi.
We define an OKN cycle as one SP followed by one QP. Variables from adjacent OKN cycles were grouped to create a vector that was used to generate a correlation matrix for each trial (40 trials, 10 Subjects × 4 Speeds). Vectors that included cycles separated by blinks were not included in generating the correlation matrix. Each correlation matrix was then analyzed using PCA. 
Principal component analysis (PCA)
PCA can be thought of as passing a line through the centroid of an 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). 
In our application, 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 nd 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 Ti as this introduces a new nonlinear relationship between Si and 1/Ti. However, PCA was performed using 1/Ti as a variable in place of Ti and we found qualitatively similar results. 
After discarding 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. 
After factor rotation, similar loading patterns were observed among the trials but with different eigenvalues. The 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). 
Independent component analysis (ICA) is a popular method used in signal processing that is similar to PCA. However, performing ICA finds the underlying components by maximizing the statistical independence of the estimated components rather than maximizing the proportion of explained variance. While this appears to make ICA a more powerful tool than PCA, to identify all of the underlying components, no more than one of them may be Gaussian (Comon, 1994). This makes it inapplicable in the case where the underlying signals represent multiple sources of Gaussian noise, and for this reason PCA was performed rather than ICA on the OKN data. 
Modeling
The Box-Jenkins (ARIMA—autoregressive integrated moving average models) approach is a common method used to analyze time series data and evaluate models based on their forecasting ability. This method assumes that values in the time series occur at equally spaced intervals, do not contain missing data, and the variance in the fluctuations over time are constant (Box & Jenkins, 1976). Whereas the duration between cycles of OKN is exceptionally variable, blinks create gaps in the time series, and there is known to be signal dependent noise in the ocular motor system (e.g., in the amplitude of saccades). Although ARIMA methods allow for statistical testing of the validity of time series models, the statistical significance can be interpreted incorrectly when the assumptions of the model are not met. In our investigation we only use the methods specified for identifying the order of an ARIMA model, such as examining the sample autocorrelation and partial autocorrelation function of variables to characterize the temporal characteristics of OKN variables. 
The autocorrelation function describes the correlation between the same variable at different points in time as a function of the time difference. Given a time series xt, the partial autocorrelation of lag k is the autocorrelation between xt and xt+k that is not already accounted for by lags one up to and including 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
When identifying the order of an ARIMA model, exponential decay of the autocorrelation function is one indication that the data can be modeled with an autoregressive model. Then, if the partial autocorrelation function has a sharp cut-off at a lag of one, the time series is identified as first order. A first order autoregressive model explains the change in xt at each lag as a function of xt-1 plus some unexplained noise, whereas a first order moving average model explains the change in xt 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. 
To model the relationships of OKN variables within cycles, we considered possible linear statistical relationships that could give the principal components observed in the data. It was not possible to use the rotated loadings directly, as the constraints in Equations 13 forced correlations among some pairs of variables, so it was necessary to interpret the loading patterns and constraints to yield a statistical model. Other possible linear relationships were explored, but they could not reproduce the distributions and relationships observed in the data. To test the model, 40 simulated sequences of the six OKN variables: xi, Si, Vi, Ti, yi, and Qi were created. The model parameters were based on robust regression analysis (Matlab function 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. 
Results
Descriptive statistics of OKN variables
Increasing stimulus speed resulted in an increase in the mean SP velocity (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. 
Figure 2
 
Descriptive statistics of OKN variables. Effect of OKN stimulus speed on (A) mean SP velocity (filled circles) and median SP duration (empty circles) and (B) mean SP amplitude (filled circles) and mean QP amplitude (empty circles). Each data point represents the global mean (or median) values averaged across 10 participants and error bars represent ± standard deviation between participants. Data points are shifted on the x-axis slightly for legibility. Representative histograms of (C) SP amplitude illustrating that SPs are not made in the direction opposite to stimulus motion and (D) QP amplitude illustrating that some QPs are made in the direction of stimulus motion and the existence of a dead zone close to 0°. Data for histograms obtained from Participant 1, stimulus speed 10°/s.
Figure 2
 
Descriptive statistics of OKN variables. Effect of OKN stimulus speed on (A) mean SP velocity (filled circles) and median SP duration (empty circles) and (B) mean SP amplitude (filled circles) and mean QP amplitude (empty circles). Each data point represents the global mean (or median) values averaged across 10 participants and error bars represent ± standard deviation between participants. Data points are shifted on the x-axis slightly for legibility. Representative histograms of (C) SP amplitude illustrating that SPs are not made in the direction opposite to stimulus motion and (D) QP amplitude illustrating that some QPs are made in the direction of stimulus motion and the existence of a dead zone close to 0°. Data for histograms obtained from Participant 1, stimulus speed 10°/s.
Increasing stimulus speed also resulted in a decrease in SP duration (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). 
Increasing stimulus speed caused mean SP amplitude to become more positive (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. 
SPs were always made in the direction of stimulus motion and SP amplitude had a unimodal distribution (Figure 2C). QPs were usually made in the negative direction, although approximately 5% were made in the positive direction and usually the variance of QP amplitude was greater than the variance of SP amplitude (36 of 40 trials). QP amplitude was sometimes bimodally distributed, with a dip around 0° extending approximately 0.8° in both directions (Figure 2D). However, some very small amplitude QPs did occasionally occur, and in trials where there were very few QPs made in the positive direction it was not possible to identify a dead zone. 
Over the course of a trial, QPs tended to compensate for SP eye displacement, thereby maintaining an average eye position. However, mean eye position tended to be held in the negative stimulus field rather than straight ahead (contraversion). The mean gaze position over the course of a trial was on average 3° in to the negative direction with no significant stimulus speed effect (F = 0.6, p = 0.50). 
Correlations between OKN variables and across cycles
Relationships between OKN variables were examined in the form of multiple pair-wise scatter plots. Plotting all 40 trials illustrated some general similarities, implying an underlying structure, but there were often remarkable differences between trials that were not systematic in any way that we could discern, and it was not possible to deduce which relationships were direct (causal) and which were due to indirect dependencies. When relationships were observed between OKN variables they appeared linear except for an approximately hyperbolic relationship between SP velocity and SP duration. Relationships observed between cycles also appeared linear, such as between the SP start position during one cycle and the next. 
The autocorrelograms of OKN variables revealed autocorrelation (Figure 3) in SP velocity and the start and end position of SPs, sometimes extending back as far as five cycles although the decay of autocorrelation was rapid. The partial autocorrelograms for SP velocity (and the start and end position of SPs) illustrated a sharp cut-off at a lag of one, implying that the autocorrelation at a lag of one cycle could explain all the higher order autocorrelations. These results implied that the time series of SP velocity and the start and end position of SPs could be modeled with a first order autoregressive model (see Methods). However, the specific values of the cross-cycle correlations were also remarkably variable between trials. We therefore used PCA to clarify the underlying structure of the correlation matrices within and between cycles of OKN. 
Figure 3
 
Autocorrelograms of OKN variables. A representative sample of autocorrelograms for (A) SP velocity, (B) SP start position, (C) QP start position, (D) SP duration, (E) SP amplitude, and (F) QP amplitude. Bold line: data obtained empirically from Participant 1, stimulus speed 40°/s. Thin line: data simulated using our proposed model (see text) with model parameters estimated from the same empirical dataset.
Figure 3
 
Autocorrelograms of OKN variables. A representative sample of autocorrelograms for (A) SP velocity, (B) SP start position, (C) QP start position, (D) SP duration, (E) SP amplitude, and (F) QP amplitude. Bold line: data obtained empirically from Participant 1, stimulus speed 40°/s. Thin line: data simulated using our proposed model (see text) with model parameters estimated from the same empirical dataset.
PCA performed on one cycle of OKN data
We first performed PCA on the simple 1 × 5 vector of OKN variables: xi, Si, yi, Qi, and xi+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. 
Figure 4
 
Results of PCA performed on a single cycle of OKN. Left panels (A, B, and C) illustrate data from Participant 1 (stimulus speed 10°/s, m = 388 cycles). Right panels (D, E, and F) illustrate data from Participant 2 (stimulus speed 10°/s, m = 201 cycles). Top panels (A and D) illustrate eye position traces during OKN presentation (bold line) and simulated eye position traces from our Markov model (thin line). Middle panels (B and E) illustrate correlation matrices of OKN variables from the corresponding empirical data set. Bottom panels (C and F) illustrate the loading patterns for the three primary components extracted from the correlation matrices, and are plotted in the order of highest to lowest eigenvalue from left to right. Empty circles indicate components extracted from empirical data, filled circles indicate components extracted from simulated data.
Figure 4
 
Results of PCA performed on a single cycle of OKN. Left panels (A, B, and C) illustrate data from Participant 1 (stimulus speed 10°/s, m = 388 cycles). Right panels (D, E, and F) illustrate data from Participant 2 (stimulus speed 10°/s, m = 201 cycles). Top panels (A and D) illustrate eye position traces during OKN presentation (bold line) and simulated eye position traces from our Markov model (thin line). Middle panels (B and E) illustrate correlation matrices of OKN variables from the corresponding empirical data set. Bottom panels (C and F) illustrate the loading patterns for the three primary components extracted from the correlation matrices, and are plotted in the order of highest to lowest eigenvalue from left to right. Empty circles indicate components extracted from empirical data, filled circles indicate components extracted from simulated data.
PCA performed on the correlation matrix generated from Participant 1 (Figure 4B) produced five components with the eigenvalues 2.41, 1.51, 1.08, 5.58 × 10−6, and 2.37 × 10−6. Due to the constraints imposed by Equations 13 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 yi, positive loading from Qi, and negative loading from xi, 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. 
Participant 2 generated a qualitatively different eye position trace to Participant 1 (Figure 4D, cf. Figure 4A), and the correlation matrices were also clearly different (Figure 4E, cf. Figure 4B). PCA performed on the correlation matrix generated by Participant 2 produced components with the eigenvalues 2.50, 1.59, 0.90, 3.09 × 10−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. 
PCA performed on the correlation matrices from all 40 trials illustrated that, despite the differences in the correlation matrices between individuals and stimulus speeds, the loading values always fell into the same three basic patterns but in different eigenvalue order. This indicated that the variability in the correlation matrices was due to differences in the eigenvalues of components (the amount of noise the component represented), but eigenvectors (the linear relationships between the OKN variables) remained similar. 
PCA performed on multiple cycles of OKN
We next performed PCA on the correlation matrices of OKN variables taken from across multiple adjacent cycles in series. The results were similar for two, three, and four cycles, and here the results from across four cycles are presented. The analysis revealed that the eigenvalues of components varied between trials but the loading patterns always fell into 13 distinct categories, corresponding to the maximum 13 DoF. Further, the 13 categories could be sorted into three broad groups that represented similar loading patterns but displaced by one, two, or three cycles, and only OKN variables from within one cycle (or just over one cycle) had high loadings in any given component (Figure 5). The three categories of components indicated that each cycle of OKN constituted three orthogonal (uncorrelated) noise processes, which we call the 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. 
Figure 5
 
Results of PCA performed on four cycles of OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series. Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 5
 
Results of PCA performed on four cycles of OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series. Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Monte-Carlo first order Markov model
Based on the PCA results we produced a model of OKN stochasticity with the same loading patterns as seen in Figure 5. Due to the constraints imposed by Equations 13 it was not possible to use the loading values to generate a model directly, and it was necessary to interpret the loading patterns to determine which correlations we would need to model and which would be forced by the constraints. 
We made four assumptions before developing the model. First, we assumed that the three uncorrelated sources of variability, or 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
In the S-component there was a high positive loading from Si, positive loading from Vi, and negative loading from xi (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 Si against xi and Vi with a mean R2 = 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 σs. We assume that the values of ŝ and σs remain constant during each trial. In the S-component there were also positive loadings from Ti and yi, sometimes high and sometimes low, but these loadings appear to be forced by the constraints in Equation 3 and Equation 2, respectively. 
Q-component
There was an equivalent pattern in the Q-component with a positive loading from Qi and a negative loading from yi (Figure 5) indicating that QP amplitude becomes more positive (decreases in magnitude) with a more negative QP start position. Evidence of Vi loading on the 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 Si. However, regression of Qi against Vi and yi revealed a linear trend with a mean R2 = 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 and standard deviation σq. We assume that the values of and σq remain constant during each trial. We interpret the other loadings in the Q-component to be forced by the constraints expressed in Equations 13
V-component
In the V-component of each cycle there was a high positive loading from Vi and corresponding negative loading from Ti with a small or absent loading from Si (Figure 5). This indicates that the 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 Vi is autocorrelated and the regression of Vi+1 against Vi reveals a clear linear relationship with a mean R2 = 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 and standard deviation σv. We assume that the values of e, , and σv remain constant during each trial. 
Model parameters
Equations 16 describe our stochastic model of the OKN system, containing three constraints and three discrete uncorrelated linear stochastic processes. Equations 1 and 6 describe the update dynamics as a first order Markov process where the variable in one cycle is dependent only on the value of the variable in the previous cycle. We performed multiple linear regressions on the OKN variables recorded from the real data using the weighted least squares method to estimate the values of the model parameters for each trial. Estimates of the standard deviation of the random processes (σs, σq, and σv) were calculated from the residuals of regression. 
The mean values of 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 (ŝ, , and ) were then recalculated for each trial, keeping the model parameters a, b, c, and d fixed. 
Monte Carlo simulations
Monte Carlo simulations of the original 40 data sets were created using the constants (a, b, c, and d), and the free parameters (e, εs, εq, and εv). The simulated data sets were created iteratively using Equations 16 with the additional constraints   where z is a small positive value to mimic the QP dead zone, which we set to 1°. 
PCA was performed on the simulated data in exactly the same way as with empirical data, and the components were sorted using the same numerical heuristics. The results of PCA performed on four cycles of simulated data are illustrated in Figure 6 and match the results of PCA performed on the empirical data (cf. Figure 5) in terms of the number of components extracted and the patterns and values of loadings. 
Figure 6
 
Results of PCA performed on four cycles of simulated OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series, simulated using the Markov model (Equations 19). Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 6
 
Results of PCA performed on four cycles of simulated OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series, simulated using the Markov model (Equations 19). Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Pair-wise scatter plots between OKN variables were created from the simulated data to qualitatively assess how well the model could describe the indirect as well as the direct relationships in the OKN system. Autocorrelograms created from the simulated data illustrated the cross-cycle correlations observed empirically (Figure 3) with only first order update dynamics. Simulated eye position traces were also created and compared to the real data (Figures 4A & 4D), illustrating that the simulation could capture both the random nature of OKN and other gross traits such as the degree of contraversion, the proportion of QPs made in the positive direction, and the mean value of OKN variables. 
The effect of stimulus speed on SP velocity in the Markov model
We considered how stimulus speed affects SP velocity in the Markov model, under the assumption that the dependence was captured either by e or by (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 g1, g2, and g3 are the gain elements of retinal slip, efference copy, and the forward gain of SP velocity, respectively. 
Figure 7
 
Proposed model of how noise in SP velocity is generated. (A) Block diagram illustrating a simple system that could generate OKN. The values g1, g2, g3 are the respective gain elements of retinal slip in the forward loop, efference copy in the feedback loop, and SP velocity in the forward loop. Vs, stimulus speed; Vs – Vi, retinal slip; Vi, SP velocity. (B) Scatter plot relating the values of and e. Data points represent values estimated from each trial by linear regression. Note that under the assumptions of the model εv = g1g3Vs and e = g3(g2g1).
Figure 7
 
Proposed model of how noise in SP velocity is generated. (A) Block diagram illustrating a simple system that could generate OKN. The values g1, g2, g3 are the respective gain elements of retinal slip in the forward loop, efference copy in the feedback loop, and SP velocity in the forward loop. Vs, stimulus speed; Vs – Vi, retinal slip; Vi, SP velocity. (B) Scatter plot relating the values of and e. Data points represent values estimated from each trial by linear regression. Note that under the assumptions of the model εv = g1g3Vs and e = g3(g2g1).
Collecting terms in Equation 10 we can see that: and substituting Equation 11 into Equation 6 illustrates how e and relate to the gain elements in the proposed control system: and Thus the model predicted that e and must be related, as they share the term g3g1. 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 −/VS from each trial (Figure 7B; R2 = 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 g3 and g2
Equation 13 also implied that the noise in the V-component may be a consequence of either g1 or g3 (or both) varying, as these gain elements are effectively multiplied by the stimulus speed to arrive at the noise process εv
Signal dependent noise
We investigated all three noise processes for signal dependence. For the S- and Q-components we compared the values of σs and σq against the mean magnitude of SPs and QPs for each trial. Orthogonal regression illustrated a linear relationship between the standard deviation of the S-component and the mean SP magnitude (Figure 8A; R2 = 0.75), and a linear relationship between the standard deviation of the Q-component and the mean QP magnitude (Figure 8B; R2 = 0.50). These relationships can be expressed in the form and where ks = 0.36, kq = 0.38, and us and uq 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. 
Figure 8
 
Scatter plots illustrating proportional noise. Scatter plots of (A) standard deviation of S-component against SP magnitude (σs = 0.36 || – 0.03), (B) standard deviation of Q-component against QP magnitude (σq = 0.38|| + 0.20), and (C) standard deviation of V-component against stimulus speed (σv = 0.12VS – 0.07). Dashed lines represent lines of orthogonal regression.
Figure 8
 
Scatter plots illustrating proportional noise. Scatter plots of (A) standard deviation of S-component against SP magnitude (σs = 0.36 || – 0.03), (B) standard deviation of Q-component against QP magnitude (σq = 0.38|| + 0.20), and (C) standard deviation of V-component against stimulus speed (σv = 0.12VS – 0.07). Dashed lines represent lines of orthogonal regression.
For the V-component, regression of σv was performed against the values of stimulus speed (VS), mean SP velocity (), and mean retinal slip. There was a weak relationship between σv and (R2 = 0.29), a moderate relationship between σv and retinal slip (R2 = 0.46), and a strong relationship between σv and VS (R2 = 0.69). We express the relationship between the standard deviation of the V-component and the stimulus speed (Figure 8C; R2 = 0.69) in the form: where kv = 0.12 and uv 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. 
Autonomous equations
A clearer picture of the entire process can be obtained by solving Equations 16 for the autonomous update dynamics of xi+1, yi+1, Si+1 and Qi+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 xi+1, yi+1, Si+1, and Qi+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). 
After the OKN stimulus starts the various OKN vectors change over time (the transient state) but will tend towards a statistical steady state where the probability distribution approaches some fixed function independent of the starting conditions (Kijima, 1997). Although the durations of the transient states were not measured directly, they can be easily computed from the model. The variables xi+1, yi+1, Si+1, and Qi+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. 
The update dynamics of SP velocity were simulated using Equation 6 with the values of 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. 
From Equations 3 and 1518, the steady state means for each OKN variable is given by:      Equations 1924 demonstrate how SP velocity affects all measured OKN variables. The model clearly predicts that mean SP and QP amplitude depends linearly on mean SP velocity, whereas mean SP duration is inversely related to mean SP velocity, with the asymptote at approximately 200 ms as → ∞. 
QP targeting
As the participants were asked to stare towards the center of the screen we might expect that the desired target for QPs is straight ahead, but the end position of QPs tends to overshoot the center and drive eye position in to the negative field of view. The data illustrated that the amplitude of QPs depend on their start position, and the Markov model predicted that there is a start position where the OKN system would generate a QP with mean amplitude of zero. We hypothesized that this position is the desired target location of QPs. To estimate this position we can solve Equation 5 when Qi = 0 to give: where yT is the desired target location. Estimating Vi with and εq(i) with , we found the average yT for each trial. Mean values of yT 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 (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 and 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. 
Discussion
In this study of human OKN we have analyzed extensively the statistical relationships between OKN eye movement variables (the start, end and amplitude of SPs and QPs, SP duration, and SP velocity), both within and across adjacent OKN cycles. All of these variables are correlated with each other, but the strength of correlations varied markedly between individuals, and the stochasticity of OKN could not be summarized by a single correlation matrix. We examined the underlying structure of the correlation matrices using PCA and found that they all have similar eigenvectors, which were clustered into three categories. However, their eigenvalues were much more variable, and components were often found in different eigenvalue orders across recording sessions. This indicated that variation in eigenvalues (the explained variance of components) rather than eigenvectors (the linear relationships between variables) were the main reason for the differences between individual correlation matrices (the correlation matrix is a linear sum of its eigenvectors weighted by their eigenvalues). We emphasize that the low number of categories of eigenvectors emerge from the data—and not because of any deliberate dimension-reducing scheme as is often employed in PCA. We conclude that OKN can be explained as a linear stochastic process with three uncorrelated noise sources, which we label the S-, Q-, and V-components (Equations 46). Although mathematically simple, the behavior of the three interacting noise processes is non-trivial. 
The most remarkable finding was that the update dynamics of SP velocity can be described by an autonomous first order Markov process. Although our model predicts an extremely short transient behavior in agreement with Abadi et al. (2005), it appears that SP velocity continues to wander in a haphazard fashion even after the system has reached steady state. When viewing a purely translational stimulus on a tangent screen the apparent rotational stimulus speed will vary depending on gaze position, but during our procedure this equates to a 6% decrease when gaze is directed 20° from the center of the screen. This difference does not seem to be able to account for the degree of variability we have observed in SP velocity, and while high gain can occur (i.e., SP velocity is not limited) it is clear that it is not always maintained. 
While studies have suggested that the resolution of slowly moving targets in the parafovea is better than for static targets (Brown, 1972) the contrast sensitivity for detecting the motion of high spatial frequency stimuli drops rapidly with speeds greater than ∼ 5°/s, and in our data we observe much greater average values of retinal slip than this. However, visual acuity was not tested during this investigation so it is not possible to determine whether participants were attempting to maintain maximal visual contrast. Further, the contrast sensitivity for detecting the motion of low spatial frequency gratings has a peak at stimulus speeds greater than 10°/s (Burr & Ross, 1982), so low spatial frequency stimuli such as those used in this investigation might allow for increased retinal slip. 
We have found that the noise associated with the 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. 
We have also found that the standard deviation of the 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. 
We originally anticipated direct correlations among variables across many OKN cycles, as a long-term adaptive mechanism could explain why average eye position does not randomly wander off to the limit of gaze. We were surprised to find only correlations across adjacent cycles indicating no memory across cycles. With hindsight, we can see that the 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. 
The average position of gaze over the course of each trial was 3° in the direction from which the stimulus originated, rather than directly ahead. An increase in this contraversion has been observed with increasing stimulus speed, when elicited with a rotating full-field patterned curtain (Garbutt, Harwood, & Harris, 2002). An increase in contraversion has also been observed during periods of perceived self-motion (circular vection) that occur during rotational OKN (Thilo, Guerraz, Bronstein, & Gresty, 2000). The lack of a significant stimulus speed effect on contraversion in our data could be the result of using a translational rather than rotational stimulus and the lack of strong linear vection. 
The start of a QP is determined by the saccadic trigger, which could be determined by SP duration (Anastasio, 1996; Carpenter, 1993, 1994), eye position (Chun & Robinson, 1978), or a combination of processes (Balaban & Ariel, 1992). Our results indicate that the threshold is dependent on both the SP start position and SP velocity and cannot be determined solely by position or duration. It is tempting to consider the dependence of SP amplitude on SP velocity as a refractory period during which QPs are not triggered. The mean value of 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. 
A possible target location for QPs has been deduced, but we predict that QPs to this target seem to constantly undershoot with a saccadic gain as low as 0.64. It has been observed that memory guided saccades made to a briefly flashed target during a period of smooth pursuit can only partially compensate for the smooth pursuit eye movement away from the target, and only if given sufficient time (Blohm, Missal, & Lefevre, 2005; Daye, Blohm, & Lefevre, 2010). These previous investigations demonstrated that saccades can fully compensate for the smooth pursuit eye movement when saccade latency is greater than 400 ms, but tend not to compensate when saccade latency is less than 200 ms. Our finding that QPs undershoot the hypothesized target location may indicate that OKN is not able to fully compensate for the eye movement away from the target location during the SP due to relatively short duration SPs. 
The distribution of SP duration is particularly interesting due to its similarity with the distribution of saccade latencies to visual targets, and stochastic models of saccade latency are often used to explain simple models of decision making. For steady-state OKN the Markov model predicts that the probability distribution of SP duration would be a ratio of two random variables Ti = Si / Vi, 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. 
The analysis we have performed here and our model might be applied to understand pathological types of nystagmus such as congenital nystagmus or gaze-dependent nystagmus due to cerebellar dysfunction. Pathological spontaneous nystagmus is often highly nonlinear and it seems unlikely that a linear analysis (such as PCA) would be successful in determining the underlying components of the data. However, nonlinear methods of analysis could be applied (e.g., principle manifolds) and higher order terms (such as acceleration) could also be included in the analysis to investigate how variables that characterize each nystagmus waveform depend on each other. 
In summary, we have examined human OKN cycle by cycle, and this leads us to consider OKN as three first order Markov processes acting in sequence to determine at which point to trigger a QP, the position at the end of the QP, and the new eye velocity at the end of the QP. Surprisingly, the Markovian update dynamics of SP velocity allows retinal slip to reach values where visual acuity should drop dramatically. This brings into question the traditional assumption that OKN is a system that solely minimizes retinal slip and requires further exploration. 
Acknowledgments
This research was funded by an Engineering and Physical Sciences Research Council Standard Research Student Award (DTG: EP/P502675/1). 
Commercial relationships: none. 
Corresponding author: Jonathan Waddington. 
Email: jonathan.waddington@plymouth.ac.uk. 
Address: Centre for Robotics and Neural Systems, Plymouth University, Plymouth, UK. 
References
Abadi R. V. Howard I. P. Ohmi M. Lee E. E. (2005). The effect of central and peripheral field stimulation on the rise time and gain of human optokinetic nystagmus. Perception, 34, 1015–1024. [PubMed]
Anastasio T. J. (1996). A random walk model of fast-phase timing during optokinetic nystagmus. Biological Cybernetics, 75, 1–9. [CrossRef]
Argyris J. Andreadis I. Pavlos G. Athanasiou M. (1998). The influence of noise on the correlation dimension of chaotic attractors. Chaos, Solitons and Fractals, 9, 343–361. [CrossRef]
Balaban C. D. Ariel M. (1992). A “beat-to-beat” interval generator for optokinetic nystagmus. Biological Cybernetics, 66, 203–216. [CrossRef] [PubMed]
Blohm G. Missal M. Lefevre P. (2005). Processing of retinal and extraretinal signals for memory-guided saccades during smooth pursuit. Journal of Neurophysiology, 93, 1510–1522. [CrossRef] [PubMed]
Blunch N. J. (2008). Introduction to structural equation modelling using SPSS and AMOS. Thousand Oaks, CA: Sage Publications, Ltd.
Box G. E. P. Jenkins G. M. (1976). Time series analysis: Forecasting and control. San Francisco, CA: Holden-Day.
Brown B. (1972). Resolution thresholds for moving targets at the fovea an in the peripheral retina. Vision Research, 12, 293–304. [CrossRef] [PubMed]
Burr D. C. Ross J. (1982). Contrast sensitivity at high velocities. Vision Research, 22, 479–484. [CrossRef] [PubMed]
Carpenter R. H. S. (1993). Distribution of quick-phase intervals in optokinetic nystagmus. Ophthalmic Research, 25, 91–93. [CrossRef] [PubMed]
Carpenter R. H. S. (1994). Express optokinetic nystagmus. In Fuchs A. F. Brandt T. Büttner U. Zee D.(Eds.), Contemporary Ocular Motor and Vestibular Research (pp. 185–187). Stuttgart, Germany: Georg Thieme.
Cheng M. Outerbridge J. S. (1974). Inter-saccadic interval analysis of optokinetic nystagmus. Vision Research, 14, 1053–1058. [CrossRef] [PubMed]
Chun K. S. Robinson D. A. (1978). A model of quick phase generation in the vestibuloocular reflex. Biological Cybernetics, 28, 209–221. [CrossRef] [PubMed]
Cohen B. Matsuo V. Raphan T. (1977). Quantitative analysis of the velocity characteristics of optokinetic nystagmus and optokinetic after-nystagmus. Journal of Physiology, 270, 321–344. [CrossRef] [PubMed]
Cohen B. Henn V. Raphan T. Dennett D. (1981). Velocity storage, nystagmus, and visual-vestibular interactions in humans. Annals of the New York Academy of Sciences, 374, 421–433. [CrossRef] [PubMed]
Comon P. (1994). Independent component analysis, a new concept? Signal Processing, 36, 287–314. [CrossRef]
Costello A. B. Osborne J. W. (2005). Best practices in exploratory factor analysis: Four recommendations for getting the most from your analysis. Practical Assessment Research & Evaluation, 10, 1–9. Retrieved from http://pareonline.net/getvn.asp?v=10&n=7.
Daye P. M. Blohm G. Lefevre P. (2010). Saccadic compensation for smooth eye and head movements during head-unrestrained two-dimensional tracking. Journal of Neurophysiology, 103, 543–556. [CrossRef] [PubMed]
Fletcher W. A. Hain T. C. Zee D. S. (1990). Optokinetic nystagmus and afternystagmus in human beings: Relationship to nonlinear processing of information about retinal slip. Experimental Brain Research, 81, 46–52. [CrossRef] [PubMed]
Garbutt S. Harwood M. Harris C. (2001). Comparison of the main sequence of reflexive saccades and the QPs of optokinetic nystagmus. British Journal of Ophthalmology, 85, 1477–1483. [CrossRef] [PubMed]
Garbutt S. Harwood M. Harris C. (2002). Anticompensatory eye position (“contraversion”) in optokinetic nystagmus. Annals of the New York Academy of Sciences, 956, 445–448. [CrossRef] [PubMed]
Garbutt S. Han Y. Kumar A. N. Harwood M. Harris C. M. Leigh R. J. (2003). Vertical optokinetic nystagmus and saccades in normal human subjects. Investigative Ophthalmology and Visual Science, 44, 3833–3841, http://www.iovs.org/content/44/9/3833. [PubMed] [Article] [CrossRef] [PubMed]
Gellman R. S. Carl J. R. Miles F. A. (1990). Short latency ocular-following responses in man. Visual Neuroscience, 5, 107–122. [CrossRef] [PubMed]
Harris C. M. Berry D. L. (2006). A distal model of congenital nystagmus as nonlinear adaptive oscillations. Nonlinear Dynamics, 44, 367–380. [CrossRef]
Hoffmann K. (1988). Responses of single neurons in the pretectum of monkeys to visual stimuli in three-dimensional space. Annals of the New York Academy of Science, 545, 180–186. [CrossRef]
Honrubia V. Downey W. L. Mitchell D. P. Ward P. H. (1968). Experimental studies on optokinetic nystagmus. Acta Oto-Laryngologica, 65, 441–448. [CrossRef] [PubMed]
Huang Y. Y. Neuhauss S. C. F. (2008). The optokinetic response in zebrafish and its applications. Frontiers in Bioscience, 13, 1899–1916. [CrossRef] [PubMed]
Ilg U. Hoffmann K. (1996). Responses of neurons of the nucleus of the optic tract and the dorsal terminal nucleus of the accessor optic tract in the awake monkey. European Journal of Neuroscience, 8, 92–105. [CrossRef] [PubMed]
Kaminiarz A. Königs K. Bremmer F. (2009). Task influences on the dynamic properties of fast eye movements. 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]
Kijima M. (1997). Markov processes for stochastic modeling. London: Chapman & Hall.
Kolarik A. J. Margrain T. H. Freeman T. C. A. (2009). Precision and accuracy of ocular following: Influence of age and type of eye movement. Experimental Brain Research, 201, 271–282. [CrossRef] [PubMed]
Land M. F. (1999). Motion and vision: Why animals move their eyes. Journal of Comparative Physiology A, 185, 341–352. [CrossRef]
Lau C. G. Y. Honrubia V. Baloh R. W. (1978). The pattern of eye movement trajectories during physiological nystagmus in humans. In Hood J. D.(Ed.), Vestibular mechanisms in health and disease (pp. 37–44). New York: Academic.
Lau C. G. Y. Honrubia V. (1986). Fast component threshold for vestibular nystagmus in the rabbit. Journal of Comparative Physiology A, 160, 585–592. [CrossRef]
Miles F. A. Kawano K. Optican L. M. (1986). Short-latency ocular following responses of monkey. I. Dependence on temporospatial properties of visual input. Journal of Neurophysiology, 10, 153–158.
Miles F. A. (1995). The sensing of optic flow by the primate optokinetic system. In Findlay J. M. Kentridge R. W. Walker R.(Eds.), Eye movement research: Mechanisms, processes and applications (pp. 47–62). Amsterdam: Elsevier.
Miles F. A. (1998). The neural processing of 3-D visual information: Evidence from eye movements. European Journal of Neuroscience, 10, 811–822. [CrossRef] [PubMed]
Robinson D. (1981). Control of eye movements. In Brookhart J. M. Mountcastle V. B.(Eds.), Handbook of physiology. The nervous system. Motor control (pp. 1275–1320). Bethesda: American Physiological Society.
Shelhamer M. (1992). Correlation dimension of optokinetic nystagmus as evidence of chaos in the oculomotor system. IEEE Transactions on Biomedical Engineering, 39, 1319–1321. [CrossRef] [PubMed]
Shelhamer M. (1996). On the correlation dimension of optokinetic nystagmus eye movements: Computational variables, filtering, nonstationarity, and surrogate data. Biological Cybernetics, 76, 237–250. [CrossRef]
Sperry R. W. (1950). Neural basis of the spontaneous optokinetic response produced by visual inversion. Journal of Comparative and Physiological Psychology, 43, 482–489. [CrossRef] [PubMed]
Thilo K. V. Guerraz M. Bronstein A. M. Gresty M. A. (2000). Changes in horizontal oculomotor behaviour coincide with a shift in visual motion perception. Neuroreport, 11, 1987–1990. [CrossRef] [PubMed]
Trillenberg P. Zee D. Shelhamer M. (2002). On the distribution of fast-phase intervals in optokinetic and vestibular nystagmus. Biological Cybernetics, 87, 68–78. [CrossRef] [PubMed]
van Beers R. J. (2008). Saccadic eye movements minimize the consequences of motor noise. PLoS ONE 3, e2070. doi: 10.1371/journal.pone.0002070.
Waddington J. Harris C.(2012) The distribution of quick phase interval durations in human optokinetic nystagmus. Experimental Brain Research. Advanced online publication. doi: 10.1007/s00221-012-3297-z.
Walls G. L. (1962). The evolutionary history of eye movements. Vision Research, 2, 69–80. [CrossRef]
Watanabe Y. Ohmura A. Shojaku H. Mizukoshi K. (1994). Optokinetic nystagmus elicited by a random dot pattern and a wide interval stripe pattern in normal subjects. Acta Oto-Laryngologica, 511, 104–108. [CrossRef]
Figure 1
 
Geometric representation of OKN and the characteristic variables. The ith cycle in a series contains one SP followed by one QP and is defined by six variables: SP start position xi, SP amplitude Si, average SP velocity Vi, SP duration Ti, QP start position yi, and QP amplitude Qi.
Figure 1
 
Geometric representation of OKN and the characteristic variables. The ith cycle in a series contains one SP followed by one QP and is defined by six variables: SP start position xi, SP amplitude Si, average SP velocity Vi, SP duration Ti, QP start position yi, and QP amplitude Qi.
Figure 2
 
Descriptive statistics of OKN variables. Effect of OKN stimulus speed on (A) mean SP velocity (filled circles) and median SP duration (empty circles) and (B) mean SP amplitude (filled circles) and mean QP amplitude (empty circles). Each data point represents the global mean (or median) values averaged across 10 participants and error bars represent ± standard deviation between participants. Data points are shifted on the x-axis slightly for legibility. Representative histograms of (C) SP amplitude illustrating that SPs are not made in the direction opposite to stimulus motion and (D) QP amplitude illustrating that some QPs are made in the direction of stimulus motion and the existence of a dead zone close to 0°. Data for histograms obtained from Participant 1, stimulus speed 10°/s.
Figure 2
 
Descriptive statistics of OKN variables. Effect of OKN stimulus speed on (A) mean SP velocity (filled circles) and median SP duration (empty circles) and (B) mean SP amplitude (filled circles) and mean QP amplitude (empty circles). Each data point represents the global mean (or median) values averaged across 10 participants and error bars represent ± standard deviation between participants. Data points are shifted on the x-axis slightly for legibility. Representative histograms of (C) SP amplitude illustrating that SPs are not made in the direction opposite to stimulus motion and (D) QP amplitude illustrating that some QPs are made in the direction of stimulus motion and the existence of a dead zone close to 0°. Data for histograms obtained from Participant 1, stimulus speed 10°/s.
Figure 3
 
Autocorrelograms of OKN variables. A representative sample of autocorrelograms for (A) SP velocity, (B) SP start position, (C) QP start position, (D) SP duration, (E) SP amplitude, and (F) QP amplitude. Bold line: data obtained empirically from Participant 1, stimulus speed 40°/s. Thin line: data simulated using our proposed model (see text) with model parameters estimated from the same empirical dataset.
Figure 3
 
Autocorrelograms of OKN variables. A representative sample of autocorrelograms for (A) SP velocity, (B) SP start position, (C) QP start position, (D) SP duration, (E) SP amplitude, and (F) QP amplitude. Bold line: data obtained empirically from Participant 1, stimulus speed 40°/s. Thin line: data simulated using our proposed model (see text) with model parameters estimated from the same empirical dataset.
Figure 4
 
Results of PCA performed on a single cycle of OKN. Left panels (A, B, and C) illustrate data from Participant 1 (stimulus speed 10°/s, m = 388 cycles). Right panels (D, E, and F) illustrate data from Participant 2 (stimulus speed 10°/s, m = 201 cycles). Top panels (A and D) illustrate eye position traces during OKN presentation (bold line) and simulated eye position traces from our Markov model (thin line). Middle panels (B and E) illustrate correlation matrices of OKN variables from the corresponding empirical data set. Bottom panels (C and F) illustrate the loading patterns for the three primary components extracted from the correlation matrices, and are plotted in the order of highest to lowest eigenvalue from left to right. Empty circles indicate components extracted from empirical data, filled circles indicate components extracted from simulated data.
Figure 4
 
Results of PCA performed on a single cycle of OKN. Left panels (A, B, and C) illustrate data from Participant 1 (stimulus speed 10°/s, m = 388 cycles). Right panels (D, E, and F) illustrate data from Participant 2 (stimulus speed 10°/s, m = 201 cycles). Top panels (A and D) illustrate eye position traces during OKN presentation (bold line) and simulated eye position traces from our Markov model (thin line). Middle panels (B and E) illustrate correlation matrices of OKN variables from the corresponding empirical data set. Bottom panels (C and F) illustrate the loading patterns for the three primary components extracted from the correlation matrices, and are plotted in the order of highest to lowest eigenvalue from left to right. Empty circles indicate components extracted from empirical data, filled circles indicate components extracted from simulated data.
Figure 5
 
Results of PCA performed on four cycles of OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series. Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 5
 
Results of PCA performed on four cycles of OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series. Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 6
 
Results of PCA performed on four cycles of simulated OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series, simulated using the Markov model (Equations 19). Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 6
 
Results of PCA performed on four cycles of simulated OKN in series. Loading patterns of 13 principal components extracted from four cycles of OKN data in series, simulated using the Markov model (Equations 19). Components are classified into three broad groups representing similar components that are shifted by one or more cycles: Q (components 1–5), S (components 6–9), and V (components 10–13).
Figure 7
 
Proposed model of how noise in SP velocity is generated. (A) Block diagram illustrating a simple system that could generate OKN. The values g1, g2, g3 are the respective gain elements of retinal slip in the forward loop, efference copy in the feedback loop, and SP velocity in the forward loop. Vs, stimulus speed; Vs – Vi, retinal slip; Vi, SP velocity. (B) Scatter plot relating the values of and e. Data points represent values estimated from each trial by linear regression. Note that under the assumptions of the model εv = g1g3Vs and e = g3(g2g1).
Figure 7
 
Proposed model of how noise in SP velocity is generated. (A) Block diagram illustrating a simple system that could generate OKN. The values g1, g2, g3 are the respective gain elements of retinal slip in the forward loop, efference copy in the feedback loop, and SP velocity in the forward loop. Vs, stimulus speed; Vs – Vi, retinal slip; Vi, SP velocity. (B) Scatter plot relating the values of and e. Data points represent values estimated from each trial by linear regression. Note that under the assumptions of the model εv = g1g3Vs and e = g3(g2g1).
Figure 8
 
Scatter plots illustrating proportional noise. Scatter plots of (A) standard deviation of S-component against SP magnitude (σs = 0.36 || – 0.03), (B) standard deviation of Q-component against QP magnitude (σq = 0.38|| + 0.20), and (C) standard deviation of V-component against stimulus speed (σv = 0.12VS – 0.07). Dashed lines represent lines of orthogonal regression.
Figure 8
 
Scatter plots illustrating proportional noise. Scatter plots of (A) standard deviation of S-component against SP magnitude (σs = 0.36 || – 0.03), (B) standard deviation of Q-component against QP magnitude (σq = 0.38|| + 0.20), and (C) standard deviation of V-component against stimulus speed (σv = 0.12VS – 0.07). Dashed lines represent lines of orthogonal regression.
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×