gensymbNot defining

#  Bayesian Dynamical Modeling of Fixational Eye Movements

Lisa Schwetlick1 Sebastian Reich2,3 & Ralf Engbert1,3
1Department of Psychology
2Institute of Mathematics
3Research Focus Cognitive Sciences
University of Potsdam  Germany

###### Abstract

Humans constantly move their eyes, even during visual fixations, where
miniature (or fixational) eye movements are produced involuntarily. Fixational
eye movements are composed of slow components (physiological drift and tremor)
and fast microsaccades. The complex dynamics of physiological drift can be
modeled qualitatively as a statistically self-avoiding random walk [24, SAW
model, see]. In this study, we implement a data assimilation approach for the
SAW model to explain quantitative differences in experimental data obtained
from high-resolution, video-based eye tracking. We present a likelihood
function for the SAW model which allows us apply Bayesian parameter estimation
at the level of individual human participants. Based on the model fits we find
a relationship between the activation predicted by the SAW model and the
occurrence of microsaccades. The latent model activation relative to
microsaccade onsets and offsets using experimental data reveals evidence for a
triggering mechanism for microsaccades. These findings suggest that the SAW
model is capable of capturing individual differences and can serve as a tool
for exploring the relationship between physiological drift and microsaccades
as the two most important components of fixational eye movements. Our results
contribute to the understanding of individual variability in microsaccade
behaviors and the role of fixational eye movements in visual information
processing.

Keywords: Fixational eye movements, visual fixation, physiological drift,
microsaccades, self-avoiding walk, likelihood function, Bayesian data
assimilation

The human eye is never truly at rest. At a macro-level the eyes move in a
sequence of fixations and saccades [54], moving different aspects of the
visual world into the receptor-dense center of the visual field. Despite the
misleading term, even during fixations, the eyes are far from stationary [14,
33, 38, 51]. Microscopic fixational movements constantly shift the visual
input over the receptors in the retina. Three components of fixational eye
movement behavior are distinguished: a slow, meandering physiological drift, a
high frequency tremor, and high velocity microsaccades [1]. The comparatively
small amplitude of tremor movement can not be resolved in current video-based
eye tracking and is therefore neglected in the following.

The function of fixational eye movements, their underlying mechanisms, as well
as the consequences for processing in the visual system have yet to be fully
understood. While fixational eye movements are found ubiquitously across
species and in all primates [32], the generation of fixational eye movements
varies between individuals and is highly characteristic [22, 23, 10, 48]. In
order to model the individually characteristic spatial statistics of
physiological drift and the relationship to microsaccades we begin with a
short discussion of the two movement components.

During physiological drift, the eyes move smoothly in a pattern that resembles
Brownian motion over small time lags [20, 9, 45], i.e., it meanders quasi-
randomly, increasing variance over time. However, a more detailed analysis
suggests that fixational eye movements represent an interesting example of
fractional Brownian motion [42]. A corresponding analysis can be carried out
by computing the mean square displacement (MSD) at different time lags [23,
28]. In Brownian motion, the MSD increases linearly with time lag. In
fixational eye movements, however, a superdiffusive tendency is found over
short time scales ($\lesssim 50$ ms), which is also referred to as
persistence. Over longer time scales ($\gtrsim 100$ ms) physiological drift is
found to be anti-persistent . This behavior can be interpreted by assuming
that fixational eye movements maximize movement over short time scales to
counteract retinal fatigue while reducing variance over longer time scale to
maintain visual fixation on a intended region of interest or object [23].

Microsaccades share their kinematic properties with their larger counterparts,
such as acceleration profile, main sequence (i.e., the linear relationship
between log amplitude and log peak velocity; see [3]), and are generated by
the same neural circuits in the brainstem [50]. Microsaccades are
distinguished mainly by their smaller amplitude, usually thresholded at
$<1\degree$ (see [47] for a review). Microsaccades occur at a rate which is
highly variable between subjects [22, 20]. Attentional mechanisms also
modulate microsaccades in both their rate, e.g., first, a reduced microsaccade
rate following target onset, and then an increased rate [22], and their
direction, e.g., the location of covert attention attracts microsaccades in
endogeneous spatial cueing [27, 22]. The general patterns of interactions of
microsaccade rates and orientations is more complex [20], however, a
computational model based on the SAW model has been proposed [17].

Early accounts of fixational eye movements conceptualized them as the result
of random firing from oculomotor units [16], or else as a nuisance component
that causes blurring, if not corrected by the visual system [44, 9]. More
recent evidence points to fixational eye movement being a necessary and useful
component of visual exploration in counteracting receptor adaptation [51, 38].
First, fixational eye movements prevent visual fading [14, 39] caused by
neural adaptation [12, 38, 39]. Although fading prevention may be achieved by
drift alone, microsaccades are much more effective at restoring vision after
fading has set in [40, 41]. Second, both drift [51, 6] and microsaccades [46]
have been found to facilitate high acuity pattern vision [29]. Specifically,
the performance of an edge detection model can be improved by the addition of
a movement component [52]. In another study, [2] use a Bayesian model of
neurons during early visual processing that simultaneously estimates eye
motion and object shape. The authors also find drift motions benefit high
acuity vision, mainly by averaging over the inhomogeneities in the retinal
receptors and receptor density. Finally, microsaccades and drift have been
found to be both corrective, i.e., moving the eyes back to the intended
fixation position, and exploratory or error-producing, i.e., moving new
details into the center of the visual field [23]. Microsaccades are typically
preceded by a reduction in drift [18, 57]. Thus, both microsaccades and drift
are functionally related and interdependent. It is a research goal of the
current study to analyze the relationship between slow fixational eye movement
components and microsaccades based on quantitative modeling while taking into
account individual differences observed in experimental data.

The SAW model integrates several of the above properties of fixational eye
movement and is biologically plausible [24, 17, 28]. The model describes
physiological drift by assuming a self-avoiding (random) walk (SAW) confined
in a movement potential that limits the movements to reproduce visual
fixation. A random walk on a lattice represents the trajectory of the eye. As
the random walk traverses the lattice, visited locations are activated (see
Figure 1). The activation represents the memory process that keeps track of
the recently visited locations. This generative model successfully reproduces
both the persistent and anti-persistent statistical properties of ocular drift
[24]. An extension of the model [28] implements neurophysiological delays and
thereby matches the characteristic oscillations found in the displacement
autocorrelation. Within the SAW model framework, [24] proposed a mechanism for
for generating microsaccades based on the activation in the SAW model.
However, the model has been a qualitative account so far, as it was not
attempted to quantitatively reproduce experimental data from human observers.

In the present work we implement the SAW model in a likelihood-based framework
in order to enable Bayesian parameter inference from fixational eye-movement
data. We estimate model parameters for individual observers, and find that the
estimated parameters represent individually characteristic spatial statistics
of physiological drift. Based on the quantitative agreement between simulated
and experimental drift movements, we explore the relation to microsaccades. We
use the model’s latent activation to investigate potential mechanisms for
triggering microsaccades. Since fixational eye movements have a strong impact
on the spatiotemporal input that the visual system processes, reproducing
these movements from a computer-implement mathematical model is essential for
a better understanding of visual functioning.

## Results

In a first step we define the SAW model and describe the computation of the
model’s likelihood function. In a second step, we use the likelihood
computation to estimate parameters for individual human participants. Finally,
we use the model to generate data and conduct posterior predictive checks and
exploratory analyses concerning the relationship of drift and microsaccades.

### The model

As a theoretical starting point, the SAW model generates a random walk which
is statistically self-avoiding [26]. The self-avoiding walk is implemented on
an $L\times L$ lattice where nodes are given by $(i,j)$ with
$i,j=1,2,3,...,L$. Each node carries some activation $a_{t}(i,j)$ at time $t$,
which can be interpreted as neural firing rates. Initial activation values are
set to $a_{t=0}(i,j)=10^{-1}$ for all $(i,j)$. At time $t=1,2,3,\dots$, first,
the current activation of each node $(i,j)$ across the field decays, according
to

$a_{t+1}(i,j)=\epsilon\cdot a_{t}(i,j)\;,$ (1)

where $\epsilon=1-(10^{\gamma})$, representing the speed of the process memory
decay. Second, activation is added to the nodes along the walker’s trajectory,
i.e.,

$a_{t+1}(i^{\star},j^{\star})=a_{t}(i^{\star},j^{\star})+1\;.$ (2)

Next, we implement a rule for activating lattice positions
$(i^{\star},j^{\star})$ along the trajectory. We define an ellipse which is
drawn such that the positions at times $t$ and $t+1$ are the foci of the
ellipse. The parameter $\rho$ represents the size of the minor axis of the
ellipse; the numerical value is set to $12$ units (see Fig. 1). The lattice
positions $\vec{v}^{\star}=(i^{\star},j^{\star})$ are defined as all lattice
sites within the ellipse, which are activated according to Eq. (2).

The discretized map $\\{a_{t}(i,j)\\}$ of neural activations can be
interpreted biologically, since grid cells have been found in entorhinal
cortex which keep track of previously visited locations [30].

In principle, the self-avoiding walk can produce persistent motion if
parameters are selected appropriately. However, during visual fixation, human
observers are able to keep the eyes at an intended target. Therefore, the
model implements a movement potential, which confines the random walk and
represents a mechanism of fixation control. As a consequence, the model can
produce anti-persistent motion on the longer time scale. Thus, the self-
avoiding motion in a potential could maintain fixation at the intended
location despite the necessity of refreshing the retinal input. The (time-
independent) confining potential $u$ is centered in the lattice and takes the
form

$u(i,j)=\lambda\frac{\left(\sqrt{(i-\frac{L}{2})^{2}+(j-\frac{L}{2})^{2}}\right)^{\nu}}{L^{\nu}}\;.$
(3)

Within this potential, motion is controlled by the sum of the self-generated
activation $a_{t}(i,j)$ and the potential $u(i,j)$, i.e.,

$q_{t}(i,j)=a_{t}(i,j)+u(i,j)$ (4)

In order to chose the next step of the walker, we compute the probability of
the next eye position as $\pi_{n}$, consisting of $q_{t}(ij)$, self generated
activation and potential, and a time-independent stepping distribution, which
controls the size of the movements, i.e.,

$\pi_{t}(i,j)={q_{t}(i,j)^{-\eta}}\,\exp\left(-\left[\left(\frac{i}{r_{i}}\right)^{\phi}+\left(\frac{j}{r_{j}}\right)^{\phi}\right]\right)\;.$
(5)

From $\pi_{t}(i,j)$ we select the eye position at time $t+1$ using a linear
selection algorithm.

As a result, in each time step follows the sequence of, first, relaxation of
the current activation, second, selection of the next eye position under
consideration of a stepping distribution, and, third, increase of the
activation values along the current trajectory. Our model’s behavior is
controlled by a number of free parameters. In this study we selected a subset
of parameters for estimation, namely $\gamma$ the speed of the relaxation, the
size of the stepping distribution $r_{i}$ and $r_{j}$, the slope of the
stepping distribution $\phi$ and the slope of the potential $\lambda$ (i.e.
$\theta=[\gamma,r_{i},r_{j},\phi,\lambda]$). We set $\rho=12$, $\nu=3$ and
$\eta=1$, in order to constrain the model and to obtain numerically stable
behavior during parameter estimation. The parameters selected for estimation
are of primary interest, as their values are interpretable quantities that may
give insight into the biological plausibility of the model.

${\color[rgb]{0,0,0}\vec{x}_{t}}$${\color[rgb]{0,0,0}\vec{x}_{t-1}}$${\color[rgb]{0,0,0}\rho}$${\color[rgb]{0,0,0}\vec{v}^{*}=(i^{*},j^{*})}$AB
Figure 1:  A Activation of lattice points within distance $\rho$ from the
eye’s trajectory. B Simulated fixational eye movement trajectories,
illustrating the SAW model. Starting from the position marked by the + and
ending at the position marked by the star, the model generates activation
along the trajectory. The movement is constrained by an activation potential
centered around the fixation position. The activation profile is initialized
by simulating some movement causing activation to be visible at locations that
have not been visited. Activation decays gradually over time.

### Likelihood function: sequential computation

The observation that makes the model compatible with a likelihood-based
approach is the fact that the likelihood $L$ at each time for each location on
the lattice is given by the selection map $\pi$. In other words given the
time-ordered data $X_{t}=\\{x_{1},x_{2}\dots,x_{t}\\}$, we can calculate the
probability of observing the walker in position $X$ at time $t$ given the
model and given all previous positions $X_{t-1}$, which corresponds to

${\cal L}_{\mathrm{M}}\left(\theta\mid
X_{t}\right)=P_{\mathrm{M}}\left(x_{t}\mid X_{t-1},\theta\right)\;.$ (6)

The likelihood of a sequence of $n$ events is therefore given by the product
of $n$ conditional probabilities ${\cal L}_{M}(\theta|X_{t})$, which is given
as

$L_{\mathrm{M}}\left(\theta\mid X_{n}\right)=\prod_{t=1}^{n}{\cal
L}_{\mathrm{M}}\left(\theta\mid X_{t}\right)\;.$ (7)

In order to estimate the parameters of the model, we use Bayes’ theorem to
compute the probability of the parameters
$\theta=[\gamma,r_{i},r_{j},\phi,\lambda]$ given the data as

$P(\theta\mid
X_{n})=\frac{L_{\mathrm{M}}(\theta\mid{X_{n}})P(\theta)}{\int_{\Theta}L_{\mathrm{M}}(\theta\mid
X_{n})P(\theta){\mathrm{d}}\theta}\;.$ (8)

A large literature exists to solve likelihood-based parameter estimation
computationally. In order to leverage the full power of the likelihood-based
approach we estimate the full Bayesian posteriors for each parameter using a
differential evolution adaptive metropolis sampling algorithm (DREAM) [37,
56]. More details are provided in the Methods section.

### Parameter estimation results

We estimated the values of the free parameters of the SAW model for each
subject independently based on data. The chosen priors (see Appendix for
details) were relatively uninformative truncated Gaussians. Figure 2 presents
the marginal posteriors for each participant in the study. In the Appendix
(Table A2) we present the point estimates and 98% confidence intervals [35]
for each parameter.

(A)(B) Figure 2:  (A) Marginal posteriors for the five estimated parameters;
the grey lines represent different participants, colored lines highlight three
participants for comparison. For all five parameters, the posteriors converged
to distinct peaks for the different participants. (B) A box plot based on the
sum log likelihood of all trials, sorted by mean log likelihood across
participants.

Biologically plausible models are designed to be grounded in real-world
biological mechanisms and processes. As such, the values of their parameters
reflect the known properties of these mechanisms. It is therefore informative
to investigate the parameter values themselves, as they permit inferences
about concrete aspects of the data generating process. First, the parameters
$r_{i}$ and $r_{j}$ correspond to the widths of the stepping distribution. The
final selection map in the model (Eq. 5), and consequently the resulting
predictions for step sizes, depend on the stepping distribution containing
$r_{i}$ and $r_{j}$, as well as the confining potential and activation. We
define the step size as the distance travelled from one measured experimental
sample to the next, i.e., 2 ms as we are using a sampling rate of 500 Hz. As
shown in Figure 3C the parameters $r_{i}$ and $r_{j}$ are strongly correlated
with the empirical step size distribution.

The parameter $\phi$ represents the slope of the stepping distribution. Higher
values are associated with a stronger tendency to move along the cardinal
directions. The estimated values for $\phi$ range from 0.9 to 1.3, indicating
that the preference for cardinal directions is stronger in some participants
than in others. The model captures individual differences in the stepping
distributions, as demonstrated by the distinct posteriors obtained for $\phi$,
$r_{i}$ and $r_{j}$.

(A)(B)(C) Figure 3:  The fit between step sizes in the model and in the data.
(A) shows the step size distribution, where highlighted lines represent
individual subjects. (B) shows the correlation between simulated and empirical
step sizes. (C) shows the correlation between the observed step size and the
parameter value of $r_{i}$ and $r_{j}$. As expected, the two parameters
correspond directly to the stepping distribution in the data.

Next, the parameter $\gamma$ relates to the speed of memory decay in the
process, where smaller, more strongly negative values cause slower decay,
i.e., longer memory of the visited locations, and larger values closer to 0
cause faster memory decay, i.e., shorter memory. It was bounded in the
estimation to a minimum of $-4$, as smaller value, i.e. less relaxation,
caused the activation in the system to continuously increase, making the model
numerically unstable. We find an average value of $\gamma$ of 3.75 which
indicates that activation decreases by 25% over the course of a trial duration
of 3 s. Parameter $\gamma$ accounts for relatively small individual
variability. For most of the participants, estimates indicate a long memory
for activation that represented the past trajectory in the system.

Finally, parameter $\lambda$ represents the slope of the confining potential.
The shape parameter of the potential function was fixed to 3, fixing the
qualitative form to a steeper gradient than in a quadratic case [24]. The
slope was considered a free parameter for the estimation. As the different
participants’ data lead to different decay and stepping parameters, the slope
of the potential needs to accommodate different resulting values of mean
system activation.

### Posterior predictive checks

Using the estimated parameters (i.e., the posterior parameter distributions),
we simulated artificial data sets of the same size as the experimental test
data sets. We know that the model is, in principle, able to recreate
statistical tendencies of ocular drift found in experimental data on a
qualitative level [24]. Posterior predictive checks serve to show that the
fitted model reproduces the expected statistical tendencies, i.e., turning
angle- and step size distributions, as well as whether individual differences
are reproduced.

Fitted primarily by parameters $r_{i}$ and $r_{j}$, the step size distribution
is represented in Figure 3 A. As defined above, the step size is the spatial
distance between two subsequent samples. In order to condense the distribution
of step sizes into one summary value, we computed the mean of each
distribution. Thus, Figure 3 B shows the correspondence between the mean true
and simulated step sizes. The model fits the step size very precisely and
perfectly captures the differences in the individual preferred step size.

Next, we investigated the absolute and relative turning angles. There exists a
preference in both ocular drift and (micro)saccades for movements in the
cardinal directions, which may be caused by the structure of the ocular
muscles [58]. This fact is captured well by the model (Figure 4A). The
relevant parameter responsible for this effect is the stepping distribution
slope $\phi$, which shapes the stepping distribution to have a stronger or
weaker preference for the cardinal directions. In order to ascertain whether
the differences in individual behavior can also be reproduced, we use the area
under the cumulative density function as a summary statistic (see the Methods
section). The result in Figure 4B shows that the model is capturing individual
differences. Furthermore, there is also a tendency for movements to be in line
with or orthogonal to the previous movement vector. While this tendency is a
lot more pronounced in the experimental data than in the simulated data, the
simulations do show a qualitative reproduction of the trend (Figure 4C).
However, it is likely that the peaks are driven by the preference in absolute
angles, as no mechanism for relative angles was built into the model. We
propose that the difference between the simulated and experimental relative
turning angle distributions reveals the part of the relative turning angle
distribution that must be accounted for by an additional, independent model
mechanism.

Empirical ocular drift data, has been found to be persistent at small time
lags and antipersistent at longer timescales. As shown in Figure 5A, the
persistent trend in our data is not as pronounced for all participants as in
comparable previous studies [23, 28]. The antipersistent period begins around
60 to 80 ms after stimulus onset. The data simulated using the fitted
parameters reflects this behavioral change well. However, it behaves more
randomly than truly antipersistent at short timescales and becomes too
strongly persistent at long timescales. Experience with the model shows that
it can produce truly antipersistent behavior by varying the free parameters.
Possible reasons why this trend was insufficiently captured by the fitted
models are that the MSD is a less dominant tendency compared to other
statistics and that our selection of free model parameters were correlated in
a way that limited the ability to fit this particular tendency. Alternatively,
the self-avoidance of the model may not be the only cause for early
persistence, suggesting a model with an explicit exploration mechanism.
Nonetheless the qualitative change in behavior is clearly present in the
simulated data. Accordingly, the correlation between the experimental data and
the simulated data show that the present model fits only account for a small
part (roughly 10%) of the individual variation (5B and C).

(A)(B)(C)(D) Figure 4:  Turning angle distributions. (A) and (C) show the
absolute and relative turning angles of one sample to the next. Thin lines
represent individual subjects and thicker lines represent the means. Simulated
data is shown in green while experimental data is red. Panels (B) and (D) show
the respective correlations between simulated and experimental data. We
characterize the angle distributions using the area under curve (AUC) of the
empirical cumulative density function (ECDF). For details see the Methods
section. (A)(B)(C) Figure 5:  Persistent and Anti-persistent behavior. (A)
shows the MSD of simulated and experimental data, where thinner lines
represent individual subjects and thicker line represent the averages of
experimental data in red and simulated data in green. Note that the lag is
given in ms, whereas the distance is given in normalized units, in order to
visualize the slopes in comparison to the identity line. (B) and (C) shows the
correlation of simulated and experimental data of linear fits of the Hurst
Exponents for the short and long timescale for individual subjects.

### Investigating microsaccades

Previous work concerning the SAW model has suggested a connection between the
self-avoiding properties of the random walk and microsaccade triggers. A
reduction in movement tends to precede microsaccades [18]. Following this
reasoning, in the model a reduction in movement corresponds to a build-up of
activation in the current position. Thus, we investigated whether the
activation predicted by the model is indeed related to the occurrence of
microsaccades. Specifically, we calculated the activation $q_{t}(i,j)$ at
times relative to microsaccade onsets $t_{\textrm{M}S-on}$ and offsets
$t_{\textrm{M}S-off}$ using the experimental data. Figure 6 shows the average
activation around the time of a microsaccade, which is consistent with the
hypothesis that high levels of activity are related to triggering
microsaccades. In order to better understand the extent of the effect, we
randomized the microsaccade onsets within each subject and computed the same
trajectory on the randomized data. We find that the activation rises more
before a microsaccade and drops more steeply compared to the randomized
controls.

Figure 6:  The average activation in the model around the time of a
microsaccade. The less opaque lines represent randomized controls, where
microsaccade onsets were randomized over the trials within one subject. In the
data we find a distinct rise in activation before the time of a microsaccade
and a drop in activation after.

### Model comparisons

The proposed model comprises three main components: the random walk with a
stepping distribution, the self-activated trace memory, and the potential. In
theory, the combination of both creates an interplay between persistence and
fixation control. In order to better understand the role of each component, we
created 3 control models by removing individual components. Specifically we
investigate

  1. 1.

the full baseline SAW model, that contains all three components

  2. 2.

a model that is a random walk in a potential (W),

  3. 3.

a model that is a random walk without a potential(W-NP),

  4. 4.

a model that is a random walk with self-avoidance but without a potential
(SAW-NP).

Note that the latter two models differ from the first two significantly in
that they are not generative and, therefore, are not biologically plausible.
In the absence of a potential it is still possible to compute the likelihood,
but it is not possible to usefully simulate data from them, as there is
nothing stopping the walker from simply walking away. We include them here,
because they provide a relevant comparison, however, these models must be
treated as substantially different concerning the conclusions they permit.

First, we compare the models in terms of their likelihood (Figure 7A). Each
model was evaluated on the test data set, using the same parameters wherever
applicable. Our findings show that in the generative models, SAW outperforms
the version without self activation (W). Removing the potential and
maintaining the self-avoidance (SAW-NP) also reduces performance. However, we
also find that the non-generative model without either potential or self-
avoidance (W-NP) shows almost identical likelihood to SAW. This suggests that
individually fitted stepping distributions and general random walk behavior by
themselves capture the most predominant features of the data. It is important
to note that the likelihood is a very general measure of model performance-
this model (W-NP) is neither biologically plausible, nor can it capture any
additional statistical properties of the data. Thus, the lesson to be drawn
from this finding is twofold: First, a high likelihood is not always a
guarantee of an appropriate model and, second, that either a potential or
self-avoidance by itself are not beneficial to model performance. Each
component, added individually, actually reduces model likelihood. The fact
that their joint effect reestablishes a similar likelihood while
simultaneously making the model more biologically plausible should be
considered a success.

To this comparison we add another model: the SAW model without individual
parameter fits by subject. We observe that on average, there is a distinct
benefit of fitting individual subjects. However, the averaged parameter model
not only reduces high likelihood values but also the number of very low
likelihood values.

(A)(B) Figure 7:  Comparisons of the different model versions. Panel A shows
the 5 different versions of the model, including different components. We show
that among biologically plausible models that the proposed mechanisms and the
individual fitting procedure confer a likelihood benefit. Among Non generative
models a pure random walk with fitted step sizes, while not biologically
plausible, also achieves a high likelihood score. Panel C shows the model
activation around the time of a saccade. The activation peak emerges only in
the models that contain a confining potential.

Second, using our comparison models, we investigate which components drive the
microsaccade effect. The very clear picture in Figure 7 shows that the
activation peak is present only for the two models with a potential. This
indicates that the effect is not driven, as we supposed, by a build-up of
self-generated activation but rather by the potential. This is consistent with
the idea that the microsaccades that drive the effect are related to control
of the fixation position.

## Discussion

Fixational eye movements display a large degree of randomness and both its
origin and purpose have been much debated in the literature. Mathematical
modeling approaches have contributed insights into the neurophysiological
origin of the movement [16, 5], the desirable or undesirable consequences of
the motion on image processing [52, 2], and the spatiotemporal statistics of
the drift trajectories [9, 24, 49]. We performed Bayesian likelihood-based
parameter inference of a self-avoiding random walk model for fixational drift
at the level of individual observers. The estimation of the parameters
converge to distinct marginal posteriors, and data simulated on the basis of
the fitted models reproduces individually characteristic behavior. In a second
step we propose a relationship between the microsaccade rate and peaks in the
model’s latent activation state. This intuition is confirmed by an
exploratory, data-driven analysis.

### Individual variability

Fixational eye movements are controlled by a complex combination of factors,
such as oculomotor control, attention, and cognition. The specific observed
patterns vary greatly by individual both for measures of ocular drift [10] and
microsaccades [48]. Our results indicate that the individual variability in
drift can at least partly be captured by the parameters of the SAW model. The
average preferred step size is a particularly pertinent example
($r_{i},r_{j}$), but also directional preferences ($\phi$) and potential slope
($\lambda$) are different between subjects. These parametrizations are
sufficient to simulate data which mirrors the characteristic features of
individuals. To our knowledge this is the first paper that models individual
variation of fixational eye movement trajectories.

The variability in drift between individuals has been found to be related to
the individual variability in acuity [11]. Individual differences in
fixational eye movement may therefore be related to a range of factors
including the precise acuity of the eye, the tendency to maintain precise
fixation [10]. Moreover, attentional preferences found in macroscopic eye
movement during facial feature viewing translate to microscopic eye movement
preferences [55].

We find a distinct benefit from using individually fitted data sets over using
a single averaged value for each parameter. This benefit becomes apparent in
better convergence properties of the parameter estimation, as the strong
distinct influences otherwise cause complex, multi-modal posteriors. Moreover,
the individual fits cause a higher overall likelihood and fit of individual
characteristics. However, while model based in averaged parameter values has a
lower likelihood, it also reduces the number of very low likelihood data
points. Thus, individuals who’s data can not be easily fitted, would benefit
from a normative influence of other subjects. This suggests the use of a
hierarchical Bayesian modelling approach in future studies.

### The confining potential

The characteristic Mean Square Displacement (MSD) of fixational drift is
persistent at small time lags and anti-persistent at longer lags [23, 28].
This is consistently the case, when averaging over large amounts of data. On
an individual level, however, this tendency is not always equally pronounced.
While the SAW model reproduces the transition to antipersistent behavior well,
it does not adequately represent the persistent component, with the current
data. As, in principle, self-avoiding random walk models can produce
persistent behavior [24, 49], this may be caused by a number of factors
including the selection of fitted free model parameters, the relatively low
persistence in the present data set, or a dominance of other statistical
tendencies. Alternatively, it is possible that the strength of the persistent
trend, which the SAW model frames as the result of self-avoidance, is in fact
amplified by an additional explicit exploration mechanism which remains to be
identified.

Exploration, or the explicit persistence, of the trajectory is highly
variable, even between trials. The confining potential in the SAW model is
static, representing a fixed intended fixation position. This is most likely a
simplification, as the intention may change over time. Experimentally, we find
a large amount of variation in the cohesion of the drift. In some trials drift
consistently occurs around a specific position. In others it is evident that
two intended fixation positions coincided over the course of the trial. In yet
others, drift consistently maintains its direction away from the starting
point. As the task and stimulus in the experiment was the same for all trials
there is little evidence to explain this variation, aside from random
variation or influence of recent past stimuli. A potential future direction
for the model could be to implement a dynamic confining potential which is
centered around a moving average of a number of recent samples. The limitation
underscores the need for more comprehensive and accurate models to better
capture the complex and individual characteristics of ocular drift behavior.

### The relationship between drift and microsaccades

Early hypotheses suggested that microsaccades serve a corrective function for
ocular drift [15, 13, 43]. Experimentally, however, no reliable correlation
data confirms this hypothesis, as microsaccades can be explorative as well as
corrective. Specifically at shorter time scales microsaccades induce
persistent correlations while at longer time scales they tend to reverse
movement to correct the fixation position [23]. Additionally, studies have
shown that during high-acuity observational tasks, participants naturally
suppress microsaccades without training [8, 59], leading to the conclusion
that microsaccades may serve no useful purpose [34]. However, more recent
research has demonstrated that microsaccades can enhance the visibility
peripheral stimuli [39], facilitate high acuity vision [46, 29] and are
responsive to task demands [31], suggesting a direct link between microsaccade
activity and visual perception.

Thus, the role of microsaccades in fixational eye movement and their
relationship with drift is not fully understood. [18] suggested that
triggering of microsaccades is triggered by to a reduction in drift movement,
i.e., low retinal image slip. This idea was further explored by the suggestion
of a relationship between the self-avoiding random walk model and microsaccade
triggers [24]. Our study provides further evidence supporting this connection.
A build-up of activation in the current position of the SAW model, is
associated with the occurrence of microsaccades. This finding is consistent
with previous work indicating that a decrease in movement precedes
microsaccades [18]. Our data-driven analysis revealed that the activation
rises more before a microsaccade and drops more steeply compared to the
randomized controls, which indicates that high levels of activation are more
likely to trigger microsaccades. By comparing model variations we find that
this trend is primarily related to the potential, indicating that the portion
of microsaccades we capture with our analysis are related to fixation control.
However, the exact mechanism behind this relationship and the potential causal
direction between activation and microsaccades requires further investigation.

### Other trajectory models

Physiological drift, as the slow component of fixational eye movement, is
often modeled as a random walk [9, 36]. Particularly when it is considered
mainly as a component in a model of visual processing [52, e.g.,], this
approximation can yield good results. However, the statistical properties of
the trajectories do differ significantly from simple randomness. To better
capture these aspects a self-avoiding random walk has been proposed [24, 49].
However the number of models that aim to predict fixational movement
trajectories is very limited. The SAW model used in this paper is one example.
Another self-avoiding random walk model was published by [49]. Instead of an
elliptical activation trace [49] implement the self-avoidance by choosing each
step direction from a continuous distribution that is weighted by the density
of recent gaze history in each direction. It achieves a similar result: at
short time scales the model is persistent, avoiding previously visited areas.
The memory of the process is limited and parametrized, allowing the authors
insight into the process memory by estimating parameters. Due to the lack of
an constraining potential, this model does not represent the subdiffusive
component at long time scales.

### Input dependence

Although initially fixational drift was often characterized as noise produced
by the oculomotor units [16], more recent evidence from electrophysiological
recordings in monkeys shows that fixational drift originates higher up in the
chain of command than the oculomotor neurons [5] and is influenced by
attentional processes [55]. Microsaccades, too, are influenced by attention
and preferentially move in the direction of the attended region when there is
covert attention [27, 22]. Thus, drift and microsaccades depend on task and
the features of the fixated target [7]. This interplay of perception and
action is consistent with the idea of active vision [25], even at the scale of
fixational eye movement.

The SAW model is stimulus independent. Other modeling approaches have
investigated the interdependence of fixational eye movements and visual
perception. It can be shown that the visual processing stream is quite capable
of dealing with the hypothesized motion blur caused by the constant
displacement of the stimulus over the receptors [44]. In fact, fixational
drift is beneficial for high acuity vision, presumably because it allows
spatial information to be redistributed into the temporal domain, modulating
the input to individual receptors [11]. Image-computable models of edge
detection can actually be improved by introducing drift [52].

Research investigating the role of motion in visual perception typically use a
random walk and are most likely quite robust to changes in the precise type of
motion. However, the experimentally observed statistical properties of drift
differ from simple random noise. More research is needed to ascertain whether
these properties convey additional benefits to visual processing. In a recent
paper [2] suggested a joined approach to infer movement and stimulus
simultaneously. The model assumes a grid of retinal cells onto which stimulus
patterns are projected and that the visual processing system does not have
access to an efference copy of the movement. Instead they use Bayesian
inference to alternatingly estimate the movement and the stimulus from the
spike rate generated by the retinal cells. The authors conclude that drift is
beneficial for high acuity vision as it helps the system to average over
inhomogeneities in the retinal receptors.

Thus, fixational eye movements play an important role for visual processing.
Integrating the fact that it is both stimulus-dependent and individually
characteristic, suggests that the movement is optimized to account for
individual physiological differences. This is consistent with the finding that
fixational eye movement and visual acuity are related [11]. A future direction
for fixational drift research may be to implement the idea that drift improves
visual acuity in a generative model, to infer the ideal motion to prevent
fading or to enable edge detection. Furthermore, although visual processing
has been found to be quite robust, the development of more accurate models of
fixational eye movement may improve the quality of models of visual
processing.

### Conclusion

In conclusion, our study contributed to fixational eye movement research
through the application of mathematical modelling and Bayesian likelihood-
based parameter inference. Our analyses suggest that self-avoiding random walk
models can effectively capture individual fixational drift behavior, as
evidenced by the convergence of distinct marginal posteriors for each
observer. Furthermore, our data-driven analysis indicates a relationship
between microsaccade rate and peaks in the model’s latent activation state,
providing further insight into the underlying neurophysiological mechanisms of
microsaccade triggering. Overall, our results provide a valuable contribution
to the understanding of fixational eye movements and highlight the importance
of individual differences in this behavior.

## Methods

### The likelihood-based modelling framework

Biologically motivated, mechanistic models allow researchers to test whether
the proposed mechanisms are capable of producing the observed behavior, to
identify which components are essential, and to explore how changes to the
system’s structure alter its output [4]. Historically, the standard approach
for cognitive models involves comparing them to time-independent summary
statistics, e.g., here it may be MSD. The likelihood-based approach offers a
number of advantages. First, it is possible to estimate the model parameters
from the data in a fully Bayesian and statistically rigorous way. The value of
the model is independent of any particular ad-hoc metric, the researcher may
want to investigate [53]. The model likelihood can further be used as a basis
for model comparison. Lastly, using the estimated parameters to simulate data,
it is possible to conduct posterior predictive checks using metrics such as
MSD to investigate whether the data constrains the model in a way that
produces the expected behavior [21]. Thus, likelihood-based parameter allows
compelling conclusions about the underlying mechanisms. Another advantage is
that Bayesian parameter estimation provides a natural way to quantify
uncertainty in the estimates, through the use of posterior distributions and
credible intervals. This can be especially useful in cases where the data are
noisy or the model is complex.

By independently estimating separate parameters for each experimental subject,
it is possible to investigate individual differences. The parameter estimation
yields a separate posterior distribution for each subject. As the parameters
represent interpretable quantities with biological counterparts, the
comparisons of the posteriors of can allow interesting insights. Additionally,
when a model is capable of representing individual differences, it speaks to
the validity of the model and its parametrization.

### Experimental data

The experimental data used for this study were eye movement trajectories
recorded using an Eyelink 2 with a sampling rate of 500Hz. Participants were
seated at a distance of 50cm to the monitor and calibrated using a 9-point
calibration grid. Each trial consisted of a fixation task, where participants
fixated a cross in the center of a white screen for 3 seconds, followed by a
scene viewing task. Here, we use only the enforced fixation data from the
first 3 seconds. Out of 50 recruited participants, 48 completed all 40 trials.
A further 6 were later excluded due to a large number of blinks. As the
experiment included a rigorous online quality control and the calibrate-
ability of subjects varies, 2 participants aborted the experiment. Trials
during which saccades or blinks were detected were repeated immediately. In
order to detect microsaccades we used a velocity-based algorithm [19]. The
data set is publicly available on the Open Science Framework
(www.osf.org/fbuxq)

### Parameter estimation

Here, we used the DREAMZS algorithm [37]. Rooted in the classical Metropolis
(Hastings) Markov Chain Monte Carlo (MCMC) algorithm, DREAMZS iteratively
explores the parameter space by sampling its position and asymptotically
converges to the true posterior distribution of the parameters. The algorithm
includes several (Markov) chains starting at arbitrary (random) positions in
the parameter space. For each chain new positions are chosen by combining
(hence “evolution”) positions of other randomly chosen chains, including their
past states.

In Bayesian parameter estimation, the unknown parameters are treated as random
variables and are assigned a prior probability distribution. This prior
distribution reflects the researcher’s initial beliefs about the likely values
of the parameters based on prior knowledge or experience. We chose relatively
broad truncated Gaussian priors, which did not constrain the estimation very
strongly. The truncated tails were chosen according to experience with the
model to prevent numerical problems in the case of extreme parameter values.

We split this data set into two separate sets: one half (20 subjects) was used
for model development and exploratory analyses. The other half (21 subjects)
was used for the final analyses and model evaluation. Each set contains data
from 27 trials for each subject. We discarded all trials where the movement
during fixation exceeded 1.2 degrees of visual angle. Due to the high
individual variability in the data this criterion excluded 7 subjects, because
too many trials were affected. The 27 trials can be split into training and
test sets, with a 2/3 to 1/3 split. Each trial was 1500 samples long, i.e.
represented a fixation of 3 seconds. This procedure finally yielded a data set
with equal numbers of samples per trial and trials per subject, facilitating
statistical analyses.

## Angle distribution comparisons

In order to compare and correlate the angle distributions of individual
participants as well as between simulated and observed data, it was necessary
to reduce the angle distribution to a single summary value. An area under to
curve (AUC) metric is a common solution to this problem. However, in our case
the distributions were already densities, i.e. their AUC was, by definition,
equal to 1. Instead we compare the AUC of the empirical cumulative density
function (ECDF). The ECDF is obtained by sorting the observations into unique
bins and calculating the cumulative probability for each. By grouping and
computing the ECDF AUC of each group, we can compare the similarity in the
peak height of the angle distributions.

## Acknowledgments

This work was supported by Deutsche Forschungsgemeinschaft (DFG) via
Collaborative Research Center SFB 1294, Project-ID 318763901. We thank Andrea
Miroslava Pomar Robles and Stefan Seelig, who collected the experimental data
set which was used in this study.

## References

  * [1] Robert G. Alexander and Susana Martinez-Conde  “Fixational Eye Movements”  In _Eye Movement Research_ Springer International Publishing, 2019, pp. 73–115  DOI: 10.1007/978-3-030-20085-5_3
  * [2] Alexander G. Anderson, Kavitha Ratnam, Austin Roorda and Bruno A. Olshausen  “High-acuity vision from retinal image motion”  In _Journal of Vision_ 20.7 Association for Research in VisionOphthalmology (ARVO), 2020, pp. 34  DOI: 10.1167/jov.20.7.34
  * [3] A. Bahill, Michael R. Clark and Lawrence Stark  “The main sequence, a tool for studying human eye movements”  In _Mathematical Biosciences_ 24.3 Elsevier, 1975, pp. 191–204  DOI: 10.1016/0025-5564(75)90075-9
  * [4] William Bechtel and Adele Abrahamsen  “Dynamic mechanistic explanation: computational modeling of circadian rhythms as an exemplar for cognitive science”  In _Studies in History and Philosophy of Science Part A_ 41.3 Elsevier BV, 2010, pp. 321–333  DOI: 10.1016/j.shpsa.2010.07.003
  * [5] Nadav Ben-Shushan, Nimrod Shaham, Mati Joshua and Yoram Burak  “Fixational drift is driven by diffusive dynamics in central neural circuitry”  In _Nature Communications_ 13.1 Springer ScienceBusiness Media LLC, 2022  DOI: 10.1038/s41467-022-29201-y
  * [6] Marco Boi, Martina Poletti, Jonathan D. Victor and Michele Rucci  “Consequences of the oculomotorcycle for the dynamics of perception”  In _Current Biology_ 27.9 Elsevier BV, 2017, pp. 1268–1277  DOI: 10.1016/j.cub.2017.03.034
  * [7] Norick R. Bowers, Josselin Gautier, Samantha Lin and Austin Roorda  “Fixational eye movements depend on task and target”  Cold Spring Harbor Laboratory, 2021  DOI: 10.1101/2021.04.14.439841
  * [8] Bruce Bridgeman and Joseph Palca  “The role of microsaccades in high acuity observational tasks”  In _Vision Research_ 20.9 Elsevier BV, 1980, pp. 813–817  DOI: 10.1016/0042-6989(80)90013-9
  * [9] Y. Burak, U. Rokni, M. Meister and H. Sompolinsky  “Bayesian model of dynamic image stabilization in the visual system”  In _Proceedings of the National Academy of Sciences_ 107.45 Proceedings of the National Academy of Sciences, 2010, pp. 19525–19530  DOI: 10.1073/pnas.1006076107
  * [10] C. Cherici, X. Kuang, M. Poletti and M. Rucci  “Precision of sustained fixation in trained and untrainedobservers”  In _Journal of Vision_ 12.6 Association for Research in VisionOphthalmology (ARVO), 2012, pp. 31–31  DOI: 10.1167/12.6.31
  * [11] Alasdair D.. Clarke, Amelia R. Hunt and Anna Hughes  “Building Bayesian Cognitive Models of Visual Foraging”  Center for Open Science, 2021  DOI: 10.31234/osf.io/zacr9
  * [12] D. Coppola and D. Purves  “The extraordinarily rapid disappearance of entopic images.”  In _Proceedings of the National Academy of Sciences_ 93.15 Proceedings of the National Academy of Sciences, 1996, pp. 8001–8004  DOI: 10.1073/pnas.93.15.8001
  * [13] Tom N. Cornsweet  “Determination of the Stimuli for Involuntary Drifts and Saccadic Eye Movements10.1016/0042-6989(76)90156-5”  In _Journal of the Optical Society of America_ 46.11 The Optical Society, 1956, pp. 987  DOI: 10.1364/josa.46.000987
  * [14] R.. Ditchburn, D.. Fender and Stella Mayne  “Vision with controlled movements of the retinal image”  In _The Journal of Physiology_ 145.1 Wiley, 1959, pp. 98–107  DOI: 10.1113/jphysiol.1959.sp006130
  * [15] R.. Ditchburn and B.. Ginsborg  “Vision with a Stabilized Retinal Image”  In _Nature_ 170.4314 Springer ScienceBusiness Media LLC, 1952, pp. 36–37  DOI: 10.1038/170036a0
  * [16] M. Eizenman, P.. Hallett and R.. Frecker  “Power spectra for ocular drift and tremor”  In _Vision Research_ 25.11 Elsevier BV, 1985, pp. 1635–1640  DOI: 10.1016/0042-6989(85)90134-8
  * [17] R. Engbert  “Computational Modeling of Collicular Integration of Perceptual Responses and Attention in Microsaccades”  In _Journal of Neuroscience_ 32.23 Society for Neuroscience, 2012, pp. 8035–8039  DOI: 10.1523/jneurosci.0808-12.2012
  * [18] R. Engbert and K. Mergenthaler  “Microsaccades are triggered by low retinal image slip”  In _Proceedings of the National Academy of Sciences_ 103.18 Proceedings of the National Academy of Sciences, 2006, pp. 7192–7197  DOI: 10.1073/pnas.0509557103
  * [19] R. Engbert, P. Sinn, K. Mergenthaler and Hans A. Trukenbrod  “Microsaccade Toolbox”, 2015
  * [20] Ralf Engbert  “Microsaccades: a microcosm for research on oculomotor control, attention, and visual perception”  In _Visual Perception - Fundamentals of Vision: Low and Mid-Level Processes in Perception_ 154 Elsevier, 2006, pp. 177–192  DOI: 10.1016/s0079-6123(06)54009-9
  * [21] Ralf Engbert  “Dynamical Models in Neurocognitive Psychology”  Springer Nature Publishing, 2021  DOI: 10.1007/978-3-030-67299-7
  * [22] Ralf Engbert and Reinhold Kliegl  “Microsaccades uncover the orientation of covert attention”  In _Vision Research_ 43.9 Elsevier BV, 2003, pp. 1035–1045  DOI: 10.1016/s0042-6989(03)00084-1
  * [23] Ralf Engbert and Reinhold Kliegl  “Microsaccades Keep the Eyes’ Balance During Fixation”  In _Psychological Science_ 15.6 SAGE Publications, 2004, pp. 431–431  DOI: 10.1111/j.0956-7976.2004.00697.x
  * [24] Ralf Engbert, Konstantin Mergenthaler, Petra Sinn and Arkady Pikovsky  “An integrated model of fixational eye movements and microsaccades”  In _Proceedings of the National Academy of Sciences_ 108.39 National Academy of Sciences, 2011, pp. 16149–16150  URL: https://www.pnas.org/content/108/39/E765/1
  * [25] John M. Findlay and Iain D. Gilchrist  “Active Vision: The Psychology of Looking and Seeing”  Oxford, UK: Oxford University Press, 2003  DOI: 10.1093/acprof:oso/9780198524793.001.0001
  * [26] Harald Freund and Peter Grassberger  “The red queen’s walk”  In _Physica A: Statistical Mechanics and its Applications_ 190.3-4 Elsevier BV, 1992, pp. 218–237  DOI: 10.1016/0378-4371(92)90033-m
  * [27] Ziad M. Hafed and James J. Clark  “Microsaccades as an overt measure of covert attention shifts”  In _Vision Research_ 42.22 Elsevier BV, 2002, pp. 2533–2545  DOI: 10.1016/S0042-6989(02)00263-8
  * [28] Carl J.. Herrmann, Ralf Metzler and Ralf Engbert  “A self-avoiding walk with neural delays as a model of fixational eye movements”  In _Scientific Reports_ 7.1 Springer ScienceBusiness Media LLC, 2017  DOI: 10.1038/s41598-017-13489-8
  * [29] Janis Intoy and Michele Rucci  “Finely tuned eye movements enhance visual acuity”  In _Nature Communications_ 11.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1038/s41467-020-14616-2
  * [30] Nathaniel J. Killian, Michael J. Jutras and Elizabeth A. Buffalo  “A map of visual space in the primate entorhinal cortex”  In _Nature_ 491.7426 Nature Research, 2012, pp. 761–764  DOI: 10.1038/nature11587
  * [31] Hee-kyoung Ko, Martina Poletti and Michele Rucci  “Microsaccades precisely relocate gaze in a high visual acuity task”  In _Nature Neuroscience_ 13.12 Springer ScienceBusiness Media LLC, 2010, pp. 1549–1553  DOI: 10.1038/nn.2663
  * [32] Hee-kyoung Ko, D. Snodderly and Martina Poletti  “Eye movements between saccades: Measuring ocular drift and tremor”  In _Vision Research_ 122 Elsevier BV, 2016, pp. 93–104  DOI: 10.1016/j.visres.2016.03.006
  * [33] Eileen Kowler  “Eye movements: The past 25years”  In _Vision Research_ 51.13 Elsevier BV, 2011, pp. 1457–1483  DOI: 10.1016/j.visres.2010.12.014
  * [34] Eileen Kowler and Robert M. Steinman  “Small saccades serve no useful purpose: Reply to a letter by R. W. Ditchburn”  In _Vision Research_ 20.3 Elsevier BV, 1980, pp. 273–276  DOI: 10.1016/0042-6989(80)90113-3
  * [35] John Kruschke  “Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan”  Academic Press, 2014
  * [36] Xutao Kuang, Martina Poletti, Jonathan D. Victor and Michele Rucci  “Temporal Encoding of Spatial Information during Active Visual Fixation”  In _Current Biology_ 22.6 Elsevier BV, 2012, pp. 510–514  DOI: 10.1016/j.cub.2012.01.050
  * [37] Eric Laloy and Jasper A. Vrugt  “High-dimensional posterior exploration of hydrologic models using multiple-try DREAM(ZS) and high-performance computing”  In _Water Resources Research_ 48.1 American Geophysical Union (AGU), 2012  DOI: 10.1029/2011wr010608
  * [38] Susana Martinez-Conde, Stephen L. Macknik and David H. Hubel  “The role of fixational eye movements in visual perception”  In _Nature Reviews Neuroscience_ 5.3 Springer ScienceBusiness Media LLC, 2004, pp. 229–240  DOI: 10.1038/nrn1348
  * [39] Susana Martinez-Conde, Stephen L. Macknik, Xoana G. Troncoso and Thomas A. Dyar  “Microsaccades Counteract Visual Fading during Fixation”  In _Neuron_ 49.2 Elsevier BV, 2006, pp. 297–305  DOI: 10.1016/j.neuron.2005.11.033
  * [40] M. McCamy, S.. Macknik and S. Martinez-Conde  “Different fixational eye movements mediate the prevention and the reversal of visual fading”  In _The Journal of Physiology_ 592.19 Wiley, 2014, pp. 4381–4394  DOI: 10.1113/jphysiol.2014.279059
  * [41] M. McCamy et al.  “Microsaccadic efficacy and contribution to foveal and peripheral vision”  In _Journal of Vision_ 12.9 Association for Research in VisionOphthalmology (ARVO), 2012, pp. 1015–1015  DOI: 10.1167/12.9.1015
  * [42] Ralf Metzler and Joseph Klafter  “The random walk’s guide to anomalous diffusion: a fractional dynamics approach”  In _Physics Reports_ 339.1 Elsevier BV, 2000, pp. 1–77  DOI: 10.1016/s0370-1573(00)00070-3
  * [43] Jacob Nachmias  “Two-Dimensional Motion of the Retinal Image during Monocular Fixation”  In _Journal of the Optical Society of America_ 49.9 The Optical Society, 1959, pp. 901  DOI: 10.1364/josa.49.000901
  * [44] Orin Packer and David R. Williams  “Blurring by fixational eye movements”  In _Vision Research_ 32.10 Elsevier BV, 1992, pp. 1931–1939  DOI: 10.1016/0042-6989(92)90052-k
  * [45] Xaq Pitkow, Haim Sompolinsky and Markus Meister  “A Neural Computation for Visual Acuity in the Presence of Eye Movements”  In _PLoS Biology_ 5.12 Public Library of Science (PLoS), 2007, pp. e331  DOI: 10.1371/journal.pbio.0050331
  * [46] Martina Poletti, Chiara Listorti and Michele Rucci  “Microscopic eye movements compensatefor nonhomogeneous vision within the fovea”  In _Current Biology_ 23.17 Elsevier BV, 2013, pp. 1691–1695  DOI: 10.1016/j.cub.2013.07.007
  * [47] Martina Poletti and Michele Rucci  “A compact field guide to the study of microsaccades: Challenges and functions”  In _Vision Research_ 118 Elsevier BV, 2016, pp. 83–97  DOI: 10.1016/j.visres.2015.01.018
  * [48] William Poynter, Megan Barber, Jason Inman and Coral Wiggins  “Individuals exhibit idiosyncratic eye-movement behavior profiles across tasks”  In _Vision Research_ 89 Elsevier BV, 2013, pp. 32–38  DOI: 10.1016/j.visres.2013.07.002
  * [49] James A. Roberts, Guy Wallis and Michael Breakspear  “Fixational eye movements during viewing of dynamic natural scenes”  In _Frontiers in Psychology_ 4 Frontiers Media SA, 2013  DOI: 10.3389/fpsyg.2013.00797
  * [50] Michele Rucci and Martina Poletti  “Control and Functions of Fixational Eye Movements”  In _Annual Review of Vision Science_ 1.1 Annual Reviews, 2015, pp. 499–518  DOI: 10.1146/annurev-vision-082114-035742
  * [51] Michele Rucci and Jonathan D. Victor  “The unsteady eye: an information-processing stage, not a bug”  In _Trends in Neurosciences_ 38.4 Elsevier BV, 2015, pp. 195–206  DOI: 10.1016/j.tins.2015.01.005
  * [52] Lynn Schmittwilken and Marianne Maertens  “Fixational eye movements enable robust edge detection”  In _Journal of Vision_ 22.8 Association for Research in VisionOphthalmology (ARVO), 2022, pp. 5  DOI: 10.1167/jov.22.8.5
  * [53] Heiko H. Schütt et al.  “Likelihood-based parameter estimation and comparison of dynamical cognitive models.”  In _Psychological Review_ 124.4 American Psychological Association (APA), 2017, pp. 505–524  DOI: 10.1037/rev0000068
  * [54] Lisa Schwetlick, Lars Oliver Martin Rothkegel, Hans Arne Trukenbrod and Ralf Engbert  “Modeling the effects of perisaccadic attention on gaze statistics during scene viewing”  In _Communications Biology_ 3.727 Nature Publishing Group, 2020, pp. 1–11  DOI: 10.1038/s42003-020-01429-8
  * [55] Natalya Shelchkova, Christie Tang and Martina Poletti  “Task-driven visual exploration at the foveal scale”  In _Proceedings of the National Academy of Sciences_ 116.12 Proceedings of the National Academy of Sciences, 2019, pp. 5811–5818  DOI: 10.1073/pnas.1812222116
  * [56] Erin M. Shockley, Jasper A. Vrugt and Carlos F. Lopez  “PyDREAM: high-dimensional parameter inference for biological models in python”  In _Bioinformatics_ 34.4, 2018, pp. 695–697  DOI: 10.1093/bioinformatics/btx626
  * [57] Petra Sinn and Ralf Engbert  “Small saccades versus microsaccades: Experimental distinction and model-based unification”  In _Vision Research_ 118 Elsevier BV, 2016, pp. 132–143  DOI: 10.1016/j.visres.2015.05.012
  * [58] David L Sparks  “The brainstem control of saccadic eye movements”  In _Nature Reviews Neuroscience_ 3.12 Nature Publishing Group UK London, 2002, pp. 952–964  DOI: ttps://doi.org/10.1038/nrn986
  * [59] Barbara J. Winterson and Han Collewijn  “Microsaccades during finely guided visuomotor tasks”  In _Vision Research_ 16.12 Elsevier BV, 1976, pp. 1387–1390  DOI: 10.1016/0042-6989(76)90156-5

##  Appendix A Appendix

###  A.1 Discretization

The model is defined on a $100\times 100$ lattice. In order to evaluate
experimental eye movement traces (which are typically given in degrees of
visual angle), we discretized the data. Each data point was multiplied by 350
and then applying the floor function. This discretization value for the entire
data set was chosen by visual inspection, as it allowed all eye movement
traces to stay within the confines of the grid, but also efficiently used the
space. Parameter values from the stepping distribution to the potential
critically depend on the specifics of this discretization.

###  A.2 Simulated data

Figure A1 shows some examples of experimental gaze trajectories which
illustrate that the model captures individual differences in the data that can
be validated by visual inspection.

Figure A1:  Some examples of fixational eye movement traces. The left column
shows the experimental data. The right column shows data simulated using the
individually fitted models. These examples qualitatively illustrate the
captured statistical properties of the drift movements.

###  A.3 Parameter recovery

Parameter recovery analyses are an important step in evaluating the stability
and reliability of a mathematical model. This is often done as a first step in
model building to ensure that the model is able to accurately capture the
underlying dynamics of the system being studied. Parameter recovery analyses
involves generating synthetic data with known (“true”) parameter values of the
model and then estimating the parameters from this simulated data. This
procedure permits the evaluation of the accuracy and precision of the
parameter estimates, as well as the sensitivity of the estimates to the choice
of estimation method and the quality and amount of data. It is an important
step to assess the overall performance of the model and identify potential
issues that may arise in the estimation process.

Figure A2:  Parameter Recovery Analysis. For each parameter we show that a
recovery analysis converges to the correct value. The top row shows the
posteriors relative to the priors and reveals strong convergence in all
parameters. The bottom row shows the caterpillar plots of the estimation
procedure, i.e., the parameter value that the three chains assumed in each
iteration. The horizontal line indicated the burn in period.

###  A.4 Priors

For the parameter estimation we chose truncated Gaussian priors. Table A1
details the numerical values that define the truncated Gaussian priors.

|  Potential Slope: $\lambda$ |  Relaxation: $\gamma$ |  Stepping Prior Radius i: $r_{i}$ |  Stepping Prior Radius j: $r_{j}$ |  Stepping Prior Slope: $\phi$
---|---|---|---|---|---
mean |  2 |  -2.5 |  5 |  5 |  1.5
sd |  2 |  0.3 |  5 |  5 |  0.5
lower |  0.3 |  -3.8 |  0.1 |  0.1 |  0.5
upper |  8 |  0 |  15 |  15 |  3
Table A1:  The parameters that define the prior distributions used during
parameter inference. For each parameter we report the mean, standard
deviation, as well as the upper and lower bounds of the truncated Gaussians.

###  A.5 Parameter estimation results

Table A2 gives the detailed point estimates for all estimated parameters for
all participants, as well as 98% confidence intervals. Note that the
participant IDs start at 20, since we estimated parameters for the final model
for participant IDs 20 to 39 only. The data for participant IDs 1 to 19 were
used for model building and are omitted here to prevent overfitting.

Subject |  Potential Slope: $\lambda$ |  Relaxation: $\gamma$ |  Stepping Prior Radius i: $r_{i}$ |  Stepping Prior Radius j: $r_{j}$ |  Stepping Prior Slope: $\phi$
---|---|---|---|---|---
|  mean | +/- |  mean | +/- |  mean | +/- |  mean | +/- |  mean | +/-
20 |  4.969 | 0.194 |  -3.585 | 0.209 |  3.304 | 0.124 |  2.664 | 0.083 |  1.080 | 0.021
21 |  4.497 | 0.162 |  -3.688 | 0.110 |  3.850 | 0.204 |  3.378 | 0.166 |  1.072 | 0.032
22 |  4.917 | 0.185 |  -3.727 | 0.072 |  3.153 | 0.106 |  2.333 | 0.081 |  1.062 | 0.028
23 |  4.533 | 0.045 |  -3.773 | 0.027 |  12.064 | 0.572 |  8.181 | 0.376 |  1.203 | 0.055
24 |  4.796 | 0.179 |  -3.693 | 0.104 |  3.402 | 0.152 |  2.644 | 0.091 |  1.077 | 0.030
25 |  4.563 | 0.167 |  -3.757 | 0.043 |  5.156 | 0.254 |  4.172 | 0.215 |  1.010 | 0.032
26 |  4.646 | 0.144 |  -3.767 | 0.033 |  7.301 | 0.266 |  4.562 | 0.195 |  1.199 | 0.026
27 |  5.100 | 0.163 |  -3.692 | 0.105 |  2.766 | 0.143 |  2.171 | 0.109 |  0.963 | 0.031
28 |  4.940 | 0.165 |  -3.700 | 0.098 |  3.558 | 0.171 |  2.152 | 0.071 |  1.052 | 0.035
29 |  5.211 | 0.165 |  -3.723 | 0.076 |  3.353 | 0.093 |  2.700 | 0.077 |  1.111 | 0.022
30 |  4.778 | 0.217 |  -3.750 | 0.050 |  3.700 | 0.133 |  2.193 | 0.104 |  0.972 | 0.028
31 |  4.621 | 0.115 |  -3.771 | 0.028 |  7.844 | 0.296 |  5.613 | 0.209 |  1.191 | 0.031
32 |  4.504 | 0.146 |  -3.759 | 0.041 |  11.304 | 0.393 |  6.701 | 0.223 |  1.217 | 0.039
33 |  4.565 | 0.223 |  -3.696 | 0.103 |  3.527 | 0.157 |  2.499 | 0.102 |  1.012 | 0.028
34 |  4.922 | 0.160 |  -3.720 | 0.080 |  4.331 | 0.207 |  3.204 | 0.165 |  1.085 | 0.037
35 |  4.693 | 0.209 |  -3.736 | 0.061 |  3.418 | 0.125 |  2.502 | 0.087 |  1.005 | 0.021
36 |  4.708 | 0.132 |  -3.765 | 0.035 |  7.780 | 0.336 |  5.351 | 0.220 |  1.220 | 0.037
37 |  4.223 | 0.104 |  -3.755 | 0.044 |  8.409 | 0.286 |  6.516 | 0.222 |  1.211 | 0.030
38 |  4.868 | 0.142 |  -3.732 | 0.067 |  3.914 | 0.116 |  2.751 | 0.091 |  1.092 | 0.025
39 |  5.061 | 0.091 |  -3.712 | 0.084 |  4.460 | 0.188 |  3.176 | 0.128 |  1.163 | 0.035
Table A2:  Point estimates and 98% confidence intervals of the 5 estimated
parameters for each subject.