Human optokinetic nystagmus: A stochastic analysis

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°/s² 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: x_i, S_i, V_i, T_i, y_i, and 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. 

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 x_i, SP amplitude S_i, average SP velocity V_i, SP duration T_i, QP start position y_i, and QP amplitude Q_i.

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 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⁻⁵) 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_i as this introduces a new nonlinear relationship between S_i and 1/T_i. However, PCA was performed using 1/T_i as a variable in place of T_i 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 x_t, the partial autocorrelation of lag k is the autocorrelation between x_t and x_t+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 x_t at each lag as a function of x_t-1 plus some unexplained noise, whereas a first order moving average model explains the change in x_t 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 1–3 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: x_i, S_i, V_i, T_i, y_i, and Q_i 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.

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.

PCA performed on one cycle of OKN data

We first performed PCA on the simple 1 × 5 vector of OKN variables: x_i, S_i, y_i, Q_i, and 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. 

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⁻⁶, and 2.37 × 10⁻⁶. 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_i, positive loading from Q_i, and negative loading from x_i, 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⁻⁶, and 3.63 × 10⁻⁷, 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).

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 1–3 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 S_i, positive loading from V_i, and negative loading from x_i (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 S_i against x_i and V_i with a mean R² = 0.25 (p ≤ 6 × 10⁻⁵ 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 T_i and y_i, 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 Q_i and a negative loading from y_i (Figure 5) indicating that QP amplitude becomes more positive (decreases in magnitude) with a more negative QP start position. Evidence of V_i 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 S_i. However, regression of Q_i against V_i and y_i revealed a linear trend with a mean R² = 0.34 (p ≤ 2.2 × 10⁻⁵ 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 σ_q. We assume that the values of qˆ 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 1–3. 

V-component

In the V-component of each cycle there was a high positive loading from V_i and corresponding negative loading from T_i with a small or absent loading from S_i (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 V_i is autocorrelated and the regression of V_i+1 against V_i reveals a clear linear relationship with a mean R² = 0.28 (p ≤ 3.3 × 10⁻⁴ 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 σ_v. We assume that the values of e, vˆ, and σ_v remain constant during each trial. 

Model parameters

Equations 1–6 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 (ŝ, qˆ, and vˆ) 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 1–6 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 1–9). 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 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₁, g₂, and g₃ 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 g₁, g₂, g₃ 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. V_s, stimulus speed; V_s – V_i, retinal slip; V_i, SP velocity. (B) Scatter plot relating the values of vˆ and e. Data points represent values estimated from each trial by linear regression. Note that under the assumptions of the model ε_v = g₁g₃V_s and e = g₃(g₂ – g₁).

Collecting terms in Equation 10 we can see that: 

 

and substituting Equation 11 into Equation 6 illustrates how 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₃g₁. 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_S from each trial (Figure 7B; R² = 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₃ and g₂. 

Equation 13 also implied that the noise in the V-component may be a consequence of either g₁ or g₃ (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; R² = 0.75), and a linear relationship between the standard deviation of the Q-component and the mean QP magnitude (Figure 8B; R² = 0.50). These relationships can be expressed in the form 

 

and 

 

where k_s = 0.36, k_q = 0.38, and u_s and u_q 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 |S¯| – 0.03), (B) standard deviation of Q-component against QP magnitude (σ_q = 0.38|Q¯| + 0.20), and (C) standard deviation of V-component against stimulus speed (σ_v = 0.12V_S – 0.07). Dashed lines represent lines of orthogonal regression.

For the V-component, regression of σ_v was performed against the values of stimulus speed (V_S), mean SP velocity (V¯), and mean retinal slip. There was a weak relationship between σ_v and V¯ (R² = 0.29), a moderate relationship between σ_v and retinal slip (R² = 0.46), and a strong relationship between σ_v and V_S (R² = 0.69). We express the relationship between the standard deviation of the V-component and the stimulus speed (Figure 8C; R² = 0.69) in the form: 

  ,

where k_v = 0.12 and u_v 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 1–6 for the autonomous update dynamics of 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). 

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 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. 

The update dynamics of SP velocity were simulated using Equation 6 with the values of 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. 

From Equations 3 and 15–18, the steady state means for each OKN variable is given by: 

 

 

 

 

 

 

 

 

 

 

 

Equations 19–24 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 V¯ → ∞. 

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 Q_i = 0 to give: 

 

where y_T is the desired target location. Estimating V_i with V¯ and ε_q(i) with qˆ, we found the average y_T for each trial. Mean values of y_T 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 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. 

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 4–6). 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 T_i = S_i / V_i, 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: [e-mail]. 

Address: Centre for Robotics and Neural Systems, Plymouth University, Plymouth, UK.