**Micromovements of the eye during visual fixations provide clues about how our visual system acquires information. The analysis of fixational eye movements can thus serve as a noninvasive means to detect age-related or pathological changes in visual processing, which can in turn reflect associated cognitive or neurological disorders. However, the utility of such diagnostic approaches relies on the quality and usability of detection methods applied for the eye movement analysis. Here, we propose a novel method for (micro)saccade detection that is resistant to high-frequency recording noise, a frequent problem in video-based eye tracking in either aged subjects or subjects suffering from a vision-related pathology. The method is fast, it does not require manual noise removal, and it can work with position, velocity, or acceleration features, or a combination thereof. The detection accuracy of the proposed method is assessed on a new dataset of manually labeled recordings acquired from 14 subjects of advanced age (69–81 years old), performing an ocular fixation task. It is demonstrated that the detection accuracy of the new method compares favorably to that of two frequently used reference methods and that it is comparable to the best of the two algorithms when tested on an existing low-noise eye-tracking dataset.**

*δ*or lower than −4.00

*δ*, astigmatism larger than 1

*δ*, visual acuity lower than 8/10 (0.1 logMAR), or dissociate phoria larger than 6

*δ*were excluded. The clinical screening was ethically approved by the “CPP Ile de France,” and it was carried out under medical supervision at the Clinical Investigation Center of the Quinze-Vingts Ophtalmologic Hospital, Paris.

*δ*was added to the refractive error. Five sessions were performed per subject, 10 trials per session. Each trial lasted 32 s, resulting in 26.66 min total recording time per subject.

*x*and

*y*traces. These jumps were removed by aligning data traces independently for the horizontal and vertical components. No other data preprocessing was performed.

*v*, with Δ

_{y}*t*= 1 ms (Engbert & Kliegl, 2003). This is equivalent to smoothing the position by a triangular normalized Bartlett window (of size 12) and then differentiating (Otero-Millan et al., 2014). Velocity smoothing using this method was very close to smoothing using Savitzky-Golay differentiation filters of second order (filter size 11; Nyström & Holmqvist, 2010).

*η*and

_{x}*η*were computed for the horizontal and vertical components using a multiple (

_{y}*λ*) of median

*SD*estimator:

*λ*controls the initial separation between the potential microsaccades from noise, and was shown to produce the best results with empirically defined values set to 5 or 6 (Engbert & Kliegl, 2003; Engbert & Mergenthaler, 2006). In our algorithm we used a permissive value

*λ*= 5 to detect most of the velocity peaks, under the constraint that consecutive peaks are separated by at least 30 ms. In addition, peaks with unnaturally high velocity (1000°/s) were discarded.

*λ*is considered as a free parameter in the Engbert and Kliegl's method, it is fixed to the relatively low value in our algorithm, since it is only used to implement initial selection of candidate events. A higher value of this parameter will result in a lower sensitivity of the algorithm (i.e., its ability to detect a true microsaccade), since it limits the number of candidate events even before noise removal. A lower value will result in a large number of candidate events without further increasing the sensitivity, since a vast majority of the microsaccades is detected at this value at low noise. Separation of putative microsaccades from noise is performed as a next step of the algorithm.

*p*= 50). The feature vectors can be easily extended to also include position and acceleration data around the peak (see Results) without any changes to the algorithm. The unsupervised clustering is performed in the high-dimensional (i.e., of dimension 2Δ

*p*+ 1 in the case of velocity-based features) space using a recently proposed fast density-based clustering method (Rodriguez & Laio, 2014). This method estimates, for each data point

*i*, its local density

*ρ*and the nearest neighbor distance

_{i}*δ*between this point and its closest point having a higher density. Cluster centers are then determined as widely separated points of highest density—that is, points with large values of

_{i}*ρ*and

_{i}*δ*.

_{i}*d*between feature vectors. Our method uses a correlation-based distance metric:

_{ij}*r*is the Pearson correlation between the two vectors. This distance corresponds to the cosine of the angle between the two vectors (in the case of centered data) and hence does not depend on the vector lengths. Since microsaccades are events with a stereotyped shape but generally different amplitudes, the correlation distance is well suited to detect them. A high correlation between two data points (e.g., between two microsaccades) corresponds to the correlation distances close to zero, whereas a low correlation (e.g., between a microsaccade and a noise event) gives

_{ij}*d*≈ 1. The space spanned by all pairwise distances between feature vectors is limited to the volume of maximal size 2, in which all microsaccades are expected to form a tight cluster, while noise events are expected to lie far from the cluster center (Figure 1A). The number of clusters may be higher than 1 depending on the chosen features (see Results).

_{ij}*d*the distance between data points

_{ij}*i*and

*j*, and

*d*is a parameter determining a characteristic distance between feature vectors in the data set (see the following section for the proposed procedure of estimation of this parameter based on the data). The nearest neighbor distance is computed as (Rodriguez & Laio, 2014):

_{c}*δ*as a function of

_{i}*ρ*for all data points for one sample subject shows an example of the decision graph, with a clearly separated cluster center (Figure 1B). The cluster center (i.e., a prototypical microsaccade) is determined automatically as the data point with the largest product

_{i}*ρ*. In the case of multiple clusters, a corresponding number of points with the largest product

_{i}δ_{i}*ρδ*is determined, and each point in the dataset is assigned to the same cluster as its closest neighbor with a higher local density (Rodriguez & Laio, 2014). This clustering approach is not iterative; it is thus fast compared to classical algorithms such as

*k*-means, and it does not suffer from related overfitting problems.

*d*from its nearest neighbor with a higher density. The value of

_{c}*d*is determined based on the density of the microsaccade cluster as follows. First, the density profile of the data is estimated using a sorted

_{c}*k*-dist graph (Ester, Kriegel, Sander, & Xu, 1996). Namely, for each data point, we calculate the distance (Equation 3) to its

*k*th nearest neighbor with

*k*= 4 (we checked that values

*k*= 1, …, 5 produce similar results). These distances are then sorted in the descending order, giving rise to the curve that reflects the density distribution of the data points (Figure 1C). The data points with largest indices (corresponding to the lowest

*k*-dist values) are those located near the cluster centers. As the data point index decreases, there is a sharp increase in the associated

*k*-dist value, indicating a fast decrease in the density. We set the characteristic distance

*d*to the

_{c}*k*-dist value corresponding to the point of the steepest decrease in density (i.e., the fastest change in the

*k*-dist curve), which is in turn calculated as the maximum of the derivative of the

*k*-dist function. The derivative is estimated using Savitzky-Golay differentiation filters of fourth order (filter size was set to 10% of the number of data points to smooth out noise). Depending on the size of the dataset, other methods can be used to detect the maximum of the derivative, including semi-automatic approaches. Thus, distances between neighboring data points are (a) smaller than

*d*for data points close to the cluster center; (b) approximately equal to

_{c}*d*near the border of the cluster (i.e., when its density drops fast); and (c) larger than

_{c}*d*far from the clusters, since noise events have lower local density (i.e., noise events are less correlated between each other than microsaccade events).

_{c}*k*-means (Otero-Millan et al., 2014), no separate cluster is assigned to noise. Rather, noise events are those that are far from the microsaccade center (although not necessarily close to each other; Dave, 1991). This essentially means that no specific assumptions are made about noise, apart from it being different from microsaccades.

*λ*parameter varied from 3 to 15, with minimal duration of microsaccade fixed to 6. The OM algorithm (Version 1.1, downloaded from http://smc.neuralcorrelate.com/sw/microsaccade-detection/) does not have free parameters, since it automatically detects the best separation between the microsaccade cluster and noise cluster in a chosen low-dimensional feature space. Therefore, when showing ROC curves, we indicate the performance of this method by a point corresponding to the optimal performance, instead of a curve.

*adjusted false positive rate*in the figures. While in the classical ROC analysis, the false positive rate equals to 0.5 when the number of false positives is equal to the number of true negatives, the same value of the adjusted rate corresponds to the case when the number of false positives equals the number of true microsaccades executed by the subject.

*p*< 10

^{–6}), and over the OM method in terms of sensitivity (i.e., it has a higher true positive rate, Wilcoxon

*p*< 10

^{–5}). Although the EK method had a higher mean sensitivity, compared to the new algorithm, it also had a much larger number of false positives. The number of false positives did not differ between the new method and the OM method.

*γ*of the ratio of velocity distributions predicted on successive time steps. If the assumptions of the model are not violated and the two noise parameters are well chosen, the application of the method results in a strong decrease of noise, after which saccades can be detected by thresholding. Because of the stochastic nature of Bayesian prediction, several passes of the algorithm are required, after which candidate microsaccades are grouped together, or discarded, based on three heuristic grouping parameters. In total, at least six parameters have to be suitably chosen to ensure correct detection. Figure 6A shows a

*γ*-plot corresponding to the velocity trace from Figure 2B for noise parameters adjusted by hand. Even though the noise is strongly reduced, the issue of the threshold assignment and the choice of grouping parameters remains. While it is possible that some parameter values could result in a high-detection performance, the choice of these parameters on a per-subject basis is a problem that makes practical application of this method and its comparison with our algorithm difficult.

*r*= 0.94,

*p*< 10

^{–6}; saccades smaller than 1°:

*r*= 0.9,

*p*< 10

^{–6}; Figure 7A and D). Mean (±

*SEM*) amplitude (Figure 7B and E) and velocity (Figure 7C and F) of fixational saccades were equal to 1.13° ± 0.01° and 85.11°/s ± 0.39°/s, respectively, with a mean frequency of 1.69 ± 0.15 Hz. Microsaccade frequency was below 1 Hz on average (Figure 7G). A vast majority of fixational saccades had a larger horizontal component (Figure 7H), independently of their amplitude. This analysis can be compared to a previous study that addressed saccades and microsaccades in subjects with ages between 60 and 70 years old (as a part of a larger cohort of younger subjects; Port et al., 2016). In their visual search task (“Where's Waldo” puzzles) Port et al. (2016) observed similar saccade frequencies, but higher amplitudes and mean velocities, as well as a large subset of saccades with a vertical bias. These differences can be attributed to the differences in the task requirements, since stimulus characteristics during visual search influence the spatial distribution of visual fixations (Najemnik & Geisler, 2008).

*x*) component of the velocity traces. Since the vast majority of microsaccades in the fixation task have a nonzero component along the horizontal direction (Figure 7D), such a reduction of the feature space is expected to increase the detection rate (since an approximately constant nonzero

*y*amplitude during a horizontal microsaccade leads to a decreased signal-to-noise ratio of the microsaccade velocity signal). Reducing the input data to only include the horizontal velocity component indeed resulted in a better detection performance (Figure 8A and B). To test whether adding position and acceleration data can further improve the detection performance, and to illustrate the ability of the algorithm to work with multiple clusters, the feature space was extended to include full position, velocity, and acceleration data along the horizontal component. Since the shape of a microsaccade in the positive

*x*direction is inverted with respect to the negative one, the clustering step results in two clusters in the feature space (Figure 8C and D). Such an extension of the feature space did not result in any further improvement in accuracy (in fact it slightly decreased the sensitivity of the algorithm, see Figure 8A), thus indicating that velocity profile contains most of the information useful for microsaccade detection. Note that for all the tested feature combinations the performance of the algorithm was better than that of the two other methods (Figure 8B).

*d*between microsaccade candidates in the correlation feature space. For each subject, the ROC-optimal value of the parameter corresponds to the point with the maximal Youden's index, given by the difference between the true positive rate and the (adjusted) false positive rate (Youden, 1950). We thus estimated the drop in the Youden's index due to an incorrectly estimated

_{c}*d*value (Figure 9A and B). An underestimation of

_{c}*d*by about 0.05 or an overestimation by about 0.1 (in units of correlation distance, 0 <

_{c}*d*< 2) resulted in the drop in performance by about 10% on average, even though it is different from subject to subject (Figure 9B). To assess the efficiency of our proposed estimation procedure for the value of the free parameter, we quantified how close the estimated value was to the ROC-optimal value (Figure 9C). Thus, our procedure for the choice of

_{c}*d*based on cluster density overestimates the optimal value by about 0.04, leading to the corresponding performance drop of only about 2.1% on average.

_{c}*k*-dist curve, so that the drop in density can be detected (once per subject). Alternatively, the experimenter can manually choose the

*d*value that corresponds to the fastest density drop by visual inspection of the

_{c}*k*-dist graph. As our analysis shows, an overestimation of this value by as much as 0.1 will not result in a substantial drop in performance.

*(pp. 213–229). New York: Springer.*

*Levels of perception**, 11 (11): 473, doi:10.1167/11.11.473. [Abstract]*

*Journal of Vision**, 8 (14): 20, 1–21, https://doi.org/10.1167/8.14.20. [PubMed] [Article]*

*Journal of Vision**, 46 (11), 987–993.*

*Journal of the Optical Society of America**, 6, 312, https://doi.org/10.3389/fnagi.2014.00312.*

*Frontiers in Aging Neuroscience**, 78 (5), 328–331, https://doi.org/10.5935/0004-2749.20150087.*

*Arquivos Brasileiros de Oftalmologia**, 12 (11), 657–664, https://doi.org/10.1016/0167-8655(91)90002-4.*

*Pattern Recognition Letters**, 235, 15–168, https://doi.org/10.1016/j.jneumeth.2014.06.020.*

*Journal of Neuroscience Methods**, 145 (1), 98–107.*

*The Journal of Physiology**, 154, 177–192, https://doi.org/10.1016/S0079-6123(06)54009-9.*

*Progress in Brain Research**, 43 (9), 1035–1045, https://doi.org/10.1016/S0042-6989(03)00084-1.*

*Vision Research**, 103 (18), 7192–7197, https://doi.org/10.1073/pnas.0509557103.*

*Proceedings of the National Academy of Sciences, USA**Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining*, 226–231.

*, 9, 167, https://doi.org/10.3389/fnsys.2015.00167.*

*Frontiers in Systems Neuroscience**, 42 (22), 2533–2545, https://doi.org/10.1016/S0042-6989(02)00263-8.*

*Vision Research**. New York: Oxford University Press.*

*Eye tracking: A comprehensive guide to methods and measures**, 13 (1), 41–48, https://doi.org/10.1016/S0966-6362(00)00089-8.*

*Gait & Posture**, 47 (6), 2478–2484, https://doi.org/10.1167/iovs.05-1311.*

*Investigative Ophthalmology and Visual Science**, 125 (9), 2005–2011, https://doi.org/10.1093/brain/awf204.*

*Brain**, 8 (14): 19, 1–25, https://doi.org/10.1167/8.14.19. [PubMed] [Article]*

*Journal of Vision**, 36 (2), 535–543, https://doi.org/10.1007/s11357-013-9582-3.*

*Age (Dordrecht, Netherlands)**, 13 (12), 1549–1553, https://doi.org/10.1038/nn.2663.*

*Nature Neuroscience**, 22 (6), 510–514, https://doi.org/10.1016/j.cub.2012.01.050.*

*Current Biology**. New York: Oxford University Press.*

*The neurology of eye movements*(5th ed.)*, 52 (3), 1275, https://doi.org/10.1167/iovs.09-4334.*

*Investigative Opthalmology & Visual Science**, 3 (3), 251–258, https://doi.org/10.1038/72961.*

*Nature Neuroscience**, 49 (2), 297–305, https://doi.org/10.1016/j.neuron.2005.11.033.*

*Neuron**, 17 (1): 13, 1–23, https://doi.org/10.1167/17.1.13. [PubMed] [Article]*

*Journal of Vision**, 76 (1), 38–42, https://doi.org/10.1034/j.1600-0420.1998.760107.x.*

*Acta Ophthalmologica Scandinavica**, 8 (3)4, 1–14, https://doi.org/10.1167/8.3.4. [PubMed] [Article]*

*Journal of Vision**, 45 (1), 272–288, https://doi.org/10.3758/s13428-012-0247-4.*

*Behavior Research Methods**, 42 (1), 188–204, https://doi.org/10.3758/BRM.42.1.188.*

*Behavior Research Methods**, 14 (2): 18, 1–17, https://doi.org/10.1167/14.2.18. [PubMed] [Article]*

*Journal of Vision**, 31 (12), 4379–4387, https://doi.org/10.1523/JNEUROSCI.2600-10.2011.*

*The Journal of Neuroscience**, 55 (1), 90–96, https://doi.org/10.3368/aoj.55.1.90.*

*American Orthoptic Journal**, 10 (3): 6, 1–18, https://doi.org/10.1167/10.3.6. [PubMed] [Article]*

*Journal of Vision**, 20 (10), 1413–1417, https://doi.org/10.1038/nn.4622.*

*Nature Neuroscience**, 118, 144–157, https://doi.org/10.1016/j.visres.2015.05.013.*

*Vision Research**, 344 (6191), 1492–1496, https://doi.org/10.1126/science.1242072.*

*Science**, 4 (8): 655, https://doi.org/10.1167/4.8.655. [Abstract]*

*Journal of Vision**, 181 (4102), 810–819, https://doi.org/10.1126/science.181.4102.810.*

*Science**, 10, 23, https://doi.org/10.3389/fnsys.2016.00023.*

*Frontiers in Systems Neuroscience**, 3 (1), 32–35, https://doi.org/10.1002/1097-0142(1950)3:1.*

*Cancer**, 34 (41), 13693–13700, https://doi.org/10.1523/JNEUROSCI.0582-14.2014.*

*Journal of Neuroscience*