Oculomotor randomness is higher in autistic children and increases with the severity of symptoms

A variety of studies have suggested that at least some children with autism spectrum disorder (ASD) view the world differently. Differences in gaze patterns as measured by eye tracking have been demonstrated during visual exploration of images and natural viewing of movies with social content. Here we analyzed the temporal randomness of saccades and blinks during natural viewing of movies, inspired by a recent measure of “randomness” applied to micro‐movements of the hand and head in ASD (Torres et al., 2013; Torres & Denisova, 2016). We analyzed a large eye‐tracking dataset of 189 ASD and 41 typically developing (TD) children (1–11 years old) who watched three movie clips with social content, each repeated twice. We found that oculomotor measures of randomness, obtained from gamma parameters of inter‐saccade intervals (ISI) and blink duration distributions, were significantly higher in the ASD group compared with the TD group and were correlated with the ADOS comparison score, reflecting increased “randomness” in more severe cases. Moreover, these measures of randomness decreased with age, as well as with higher cognitive scores in both groups and were consistent across repeated viewing of each movie clip. Highly “random” eye movements in ASD children could be associated with high “neural variability” or noise, poor sensory‐motor control, or weak engagement with the movies. These findings could contribute to the future development of oculomotor biomarkers as part of an integrative diagnostic tool for ASD.


INTRODUCTION
Eye movements provide a rich source of information on sensory, perceptual, cognitive, and motor brain processes that govern a variety of voluntary and involuntary behaviors.As such, they have been studied extensively in people with autism spectrum disorder (ASD).Early eyetracking studies of ASD focused on gaze patterns related to social interaction (Klin et al., 2003), motivated by the clinical observations of reduced eye contact and joint attention in ASD (Senju & Johnson, 2009).Abnormal eye contact, as part of non-verbal communication, is included as an indicator of ASD in DSM-5 (American Psychiatric Association, 2013) and the ADOS (Lord et al., 1989;Lord et al., 2000).More recent studies recorded eye movements while the participants watched movies or pictures.The general finding could be described as reduced social interest or avoiding eye gaze with social content.This includes gaze preference in children with ASD for dynamic geometric shapes over social videos (Pierce et al., 2016), fewer fixations towards faces (Chawarska et al., 2012;Frazier et al., 2017;Shic et al., 2011) and specifically towards the eyes (Frazier et al., 2017;Jones et al., 2008;Jones & Klin, 2013), as well as a lack of preference for biological motion found already in 2-year-old infants (Klin et al., 2009).
Other studies have measured the divergence of gaze patterns in ASD individuals from those of TD individuals when naturally watching movie clips (Avni et al., 2020;Ramot et al., 2020).These studies mostly used movies with social content, hypothesizing that differences in social preferences and visual exploration would generate different gaze patterns in ASD versus TD children.The result of these studies revealed that the gaze patterns of TD individuals are very similar to their group average, whereas gaze patterns of ASD individuals are more idiosyncratic and differ substantially from their group average as well as from the TD average (Avni et al., 2020;Ramot et al., 2020).Such gaze differences could be due to differences in the way ASD individuals interpret the movies, to altered preferences and interests including an attraction to low-level visual saliency (Wang et al., 2018;Wang et al., 2015), or avoidance of social content (e.g., faces and eyes) (Baron-Cohen et al., 1997;Chawarska et al., 2012;Chita-Tegmark, 2016;Constantino et al., 2017;Jones et al., 2008;Jones & Klin, 2013;Klin et al., 2002;Papagiannopoulou et al., 2014;Riby & Hancock, 2009;Rice et al., 2012;Q. Wang et al., 2018).Alternatively, gaze differences could be due, perhaps partially, to basic differences in oculomotor control.Although some studies have reported abnormal saccade kinematics (including latency, duration, and accuracy) in ASD (Caldani et al., 2023;Goldberg et al., 2002;Luna et al., 2007;Luna et al., 2008;Schmitt et al., 2014;Wilkes et al., 2015), others have not (Avni et al., 2021;Johnson et al., 2016).Yet, another alternative is that the eye gaze patterns were determined not only by the dynamic changes in the stimulus but also due to excessive sensorimotor variability, found in autism in other domains such as body movements (Torres et al., 2013;Torres & Denisova, 2016); which was not captured by the analyses of basic saccade and fixation properties.
Previous studies have suggested that ASD may be characterized by excessive levels of neural variability (Dinstein et al., 2015).In the motor domain, Torres (2013) assessed the randomness of hand movement during simple pointing tasks, by examining the distribution of its peak velocity across trials.The distribution was fit with a gamma function and groups were clustered in the Gamma plane using the gamma parameters shape and scale.The results showed higher randomness in the ASD group compared with controls, that is, the gamma parameters showed lower gamma shape values approaching a memoryless exponential distribution with higher scale values, suggesting an elevated noise-to-signal ratio (NSR) in the ASD group.Moreover, the effect was stronger in the more severe cases, including low IQ, and nonverbal individuals (Torres et al., 2013).In later studies, Torres and colleagues extended their findings to micro head movements during MRI scanning in a large sample (N > 1000) (Torres & Denisova, 2016), and to an alternative method for assessing randomness using the "R" parameter (Wu et al., 2018).Zhao, Xing, et al. (2022) assessed the randomness of eye movements in ASD and found excessive and more random eye movements in ASD children compared with TD during face-to-face conversations as quantified with an entropy measure (Shannon, 1948;Zhao, Xing, et al., 2022).
In the current study, we applied the "gamma fit" analysis of Torres et al. (2013) to a large dataset of eyetracking recordings from 230 children naturally watching short video clips.A subset of this data has already been used by studies reporting idiosyncratic gaze patterns (Avni et al., 2020) and typical basic oculomotor function (Avni et al., 2021) in young children with ASD.Here we extended the previous analysis to examine the randomness of saccade timing (inter-saccadic intervals) and eye blink properties.Our results show a striking similarity to those obtained by Torres et al. (2013) and Torres and Denisova (2016) when examining micro-movements of the hand and head, including a correlation with severity, and group differences that allow sensitive classification.

Participants
A total of 230 children, mean age 4.50 years (SD = 2.00): N = 189 children diagnosed with ASD according to DSM-5 criteria, and N = 41 TD.All of the children underwent a routine eye examination, and none of the subjects were visually impaired.All children were tested at the Azrieli National Center for Autism Research (Meiri et al., 2017), with analyses of a subset of the current data previously reported (Avni et al., 2020(Avni et al., , 2021;;Meiri et al., 2017).The study was approved by the Soroka Medical Center Helsinki committee and the Ben Gurion University internal review board committee.Written informed consent was obtained from all parents.For more information, see Avni et al. (2020Avni et al. ( , 2021)).

Clinical evaluation
All ASD children were diagnosed by a developmental psychologist and a child psychiatrist or pediatric neurologist, according to DSM-5 criteria (American Psychiatric Association, 2013).Of the 189 ASD children, 162 completed the Autism Diagnostic Observation Schedule 2 (ADOS-2) (Lord, 2012).Raw ADOS scores were converted into calibrated severity scores (CSS), also known as the ADOS-2 comparison score, to enable a fair comparison of core ASD symptoms severity across multiple ages and language abilities (Hus Bal & Lord, 2015;Wiggins et al., 2019).The ASD mean severity across ASD children was 6.69 (SD = 2.25), ranging from 1 to 10, where 10 represents the highest ASD severity (Hus et al., 2014).Cognitive evaluations were also performed in 153 of the ASD children using one or more of the following tests: the Bayley Scales of Infant and Toddler Development (Bayley, 2006), the Wechsler Preschool and Primary Scale of Intelligence (Wechsler, 2002), or the Mullen Scales of Early Learning (Mullen, 1995).In cases where children were invited to follow-up visits, we used the average score per test.To ensure that TD children did not exhibit significant ASD symptoms, all parents completed the social responsiveness scale (SRS) (Constantino & Gruber, 2012).The mean SRS score of the TD group was 34.38 (SD = 14.25), ranging from 5 to 66 (all scores were under 75).

Stimuli and procedure
Each session consisted of 3 movies with social content: (a) Naturalistic (Nat): an unedited home video that shows the daily interaction between two sisters in a domestic room; the soundtrack is in Hebrew; (b) "Jungle-Book" (JB): a scene from the original Walt Disney cartoon movie with the song "I wanna be like you"; the background showed many details, with monkeys dancing and singing, with a soundtrack in English; (c) "Jack-Jack Attack" (JJ): a scene from the animated Pixar movie in which a babysitter interacts with a baby with supernatural powers, with a soundtrack in English.Each movie was $1.5 min long and was presented twice, 6 movies in total in a fixed order, $10 min per session including calibration tasks (with 9-10 target locations) between movies (for more details, see Avni et al. (2020Avni et al. ( , 2021)).

Apparatus
The children sat on a chair adjusted to their age; young children were seated in a safety seat with straps.The screen (with a size of 32 Â 27 cm and a resolution of 1280 Â 1024 pixels at 60 Hz) was attached to the wall by an adjustable arm so that the children could rest their heads on the chair to reduce head movements; the sitting distance was 60 cm.Eye movements were recorded using the EyeLink 1000+ head-free eye-tracking system (SR Research, Ontario, Canada), with a sampling rate of 500 Hz, with the sticker mode used for head tracking.
The EyeLink camera was located below the display screen.All recordings were done from the left eye.A fivepoint calibration was used before each session.

Saccade detection
Saccades were detected using the algorithm of (Engbert & Kliegl, 2003).To enable a good detection of saccades, data were smoothed with a local linear regression fitting (the LOWESS method), with a size of 25 ms before saccade extraction.The included saccade amplitude range was 0.08 -40 (to ensure that saccades will not exceed the screen size), a maximum velocity below 1000 /sec, and a minimum duration of 9 ms constraints were applied.After saccades were detected, we computed the inter-saccadic intervals (ISI) from the saccade onset to the following saccade onset.This onsetto-onset definition of ISI is not standard, but it is more appropriate here since we are specifically interested in the saccadic events and their temporal intervals.We removed ISIs shorter than 50 ms and fit the remaining ISIs with the gamma function.In addition, we computed the basic saccade properties for each movie clip, that is, the saccade rate was calculated as the number of saccades per movie divided by the movie's recorded length (calculated from valid samples including the detected blink periods), the mean saccade amplitude was calculated as the angular displacement of the saccade in each movie clip, the mean saccade duration, and the mean saccade peak velocity.

Blink detection
Blinks were detected as periods when the pupil size was marked as zero.In addition, we used the EyeLink blink detection events to enhance the reliability of identifying periods of missing data as actual blinks.To spot the onset and offset of each blink, the vertical traces were examined in data windows of 100 ms before and 150 ms after each missing data period (Yablonski et al., 2017).The blink onset and offset were detected when the vertical trace passed the 50th percentile of the selected data window around each blink.In addition, to avoid false detection of blinks in periods of fragmented recording, the blinks whose end time of the previous blink and the start of the next blink are shorter than 50 ms were combined into a single period that could be too long for a blink.Blinks outside the range of 250-700 ms were rejected as possible measurement noise (Yablonski et al., 2017).We fit the duration of the blinks, extracted according to the criteria described above, with the gamma function.In addition, we calculated the basic eye blink properties for each movie clip: the blink rate as the number of blinks divided by the length of each movie, the mean blink duration, and the mean inter-blink interval (IBI) calculated from the blink onset to the following blink onset.

Recording eye-tracking screening quality
We used two measures for testing the data for low quality: (1) the percentage of valid samples per movie, and (2) the percentage of valid samples between saccades.For (1), we rejected a movie clip data set when the percentage of valid pupil samples relative to the total samples, excluding the blink periods (see Section 2 of blink detection) was less than 50%; we used this method to reject data in all analyses reported in the main text.For (2), we rejected inter-saccadic intervals with less than 80% valid samples; the results are shown in Figure S1 and S3.

Gamma fit analysis
We adapted the Torres et al. (2013) gamma analysis originally used for kinesthetic data to measure the extent of oculomotor randomness.The gamma distribution is characterized by two parameters: shape (a) and scale (b).The Gamma distribution mean is represented by μ ¼ a Â b, and the variance by σ 2 ¼ a Â b 2 .The skewness of the gamma distribution is controlled by the shape parameter (2= ffiffi ffi a p ), whereas the scale parameter (b) represents the dispersion according to the Fano Factor; a higher scale value corresponds to a higher noise-tosignal-ratio (NSR) and higher dispersion (Torres et al., 2013;Torres & Denisova, 2016).The gamma distribution can vary from an exponential distribution when the shape parameter is equal to one (a = 1), and it approaches a normal distribution with large shape values.Here, we fit a gamma distribution, using a built-in MATLAB function (The MathWorks, 2020), to the: (1) inter-saccadic intervals (ISI), and (2) eye-blink duration (BD) and (3) the basic oculomotor properties including saccade peak velocity, amplitude, duration, and inter-blink intervals.Gamma function parameters: shape (a) and scale (b), were calculated per recording of one movie clip and were averaged across movies.We included gamma function values within a 95% confidence interval (CI) and excluded gamma function values with CI larger than 2 standard deviations of the CI mean to ensure the fit quality.These conditions left us with the minimum value of 46 saccades or 2 blinks per gamma fitting.In addition, we did not include gamma fitting measures for data sets with less than 50% valid samples (see "Eye-tracking recording quality screening" above, and Figure S2).In our results, a gamma distribution with a lower shape value and a higher scale value corresponds to a more exponential-like distribution that is indicative of a more random process.For statistical convenience, we created a single parameter representing the gamma parameter ratio (GR).The GR was computed by dividing the gamma scale by the gamma shape and transforming it to log 10 units such that a higher GR corresponded to a more random process.

Assessment of the statistical significance
We used Type III Analysis of Variance (ANOVA) with Satterthwaite's method to compute the denominator degrees of freedom and p values for the F-test (Kuznetsova et al., 2017).The linear mixed effect model (LMER) was fitted to the data set with the average result per movie for a subject, the response variable was set to the gamma parameter ratio (e.g., GR-ISI and GR-BD) or the basic saccade and blink properties.The predictors were set as the main effects and the interaction of the movie type and the group type, and the subjects were set as the random effect.In addition, we used a linear model (LM) on the data set containing the average result across movies per subject, the response variable was set as the GR-ISI or GR-BD, and the predictor was set as the ADOS comparison scores (i.e., the severity).We also used LM with GR-ISI as the response variable, and the predictors were set as the main effects and the interaction of the rounded age groups (to the nearest whole number) and the group type.After the models were made, we used ANOVA to generate the F-test results.Pearson's correlation was tested between GR and the cognitive scores (i.e., Bayley, Mullen, and WPPSI) and severity.Since this study, despite a specific original goal, also measured multiple effects, we applied the FDR method (Benjamini & Hochberg, 1995) to assess the effect of multiple comparisons on the statistical significance.We took a conservative approach and applied the FDR to all measures, even though some involved assessing the data quality or replicating previous findings.We therefore reported the original p-values and the FDR-corrected threshold for significance when relevant.All statistics were performed using R (Bates et al., 2015;Kuznetsova et al., 2017;R Core Team, 2021).

RESULTS
We analyzed the temporal aspects of eye movements in children while they watched three types of short video clips.We started by investigating basic oculomotor parameters, followed by Gamma distribution analyses of (1) the time intervals between saccades (inter-saccadic intervals), and (2) the duration of blinks.

Basic saccade and eye blink characteristics
The basic saccade properties, including the mean saccade peak velocity, amplitude, duration, rate, and ISI were calculated and are summarized in Table 1 (see also Figure S4a-e).Similar analyses were previously done on a subset of the data (Avni et al., 2020) and were repeated here for comparison with the previous findings, as well as with the gamma distribution analyses.As shown, the basic saccade properties did not show significant group effects, except for a weak group effect of the mean inter-saccadic interval (ISI) (F(1, 227) = 4.48, p = 0.035, see Table 1, Figure S4a), which was found to be non-significant after FDR correction.However, the saccade amplitude, duration, and peak velocity showed significant movie effects, with a significant movie-group interaction for the mean saccade duration (Table 1).The basic blink properties are summarized in Table 1 (see also Figure S4g-i), showing a significant group effect for all the basic blink properties including the mean inter-blink intervals (IBI), blink duration, and the blink rate; however, following FDR showing the results of the group (ASD and DT) and movie (JB, JJ, and Nat) main effects and their interaction (see the Methods).The pairwise comparisons in the right columns were averaged over the group level, the degrees of freedom were estimated using the Kenward-Roger method, and p-values adjusted using the Tukey method for comparing a family of 3 estimates.The saccade amplitude and duration showed a significant interaction, and their post hoc test is shown at the bottom of each oculomotor measure.The main effects were further tested for multiple comparisons using FDR (Benjamini & Hochberg, 1995); p-values below p < 0.035 were found to be significant after FDR and were marked in bold.ISI and IBI stand for the inter-saccadic and inter-blink intervals, respectively.correction, the IBI group effect was found to be nonsignificant.

Gamma distribution of inter-saccadic intervals
The inter-saccadic intervals (ISI) distribution was fit with the gamma function, and two measures were obtained: gamma shape and gamma scale (see Section 2).The Gamma fit measures for both groups, ASD and TD, are shown in Figure 1; the statistics are summarized in Table 2. Figure 1a shows a scatter plot of the two gamma parameters, shape and scale, in log units, with one point per observer, averaged across the 3 movie clips.As shown, the ASD children had a higher gamma scale and a lower gamma shape in comparison with the TD children, the ASD group results correspond to a more exponential-like distribution, implying a more random process (see examples of the ISI distribution with a gamma fit in Figure S1).The group averages of the gamma parameters across subjects and movies are shown in Figure 1b for the Jungle Book (JB), Jack-Jack (JJ), and Naturalistic (Nat) clips, all containing social interactions.As shown, there was a clear difference between groups and different randomness levels per movie clip.
The highest randomness was found for the Jungle Book movie clip in both groups.For statistical assessment, we computed a single combined measure for randomness, GR, as the log gamma parameters ratio (the gamma scale divided by shape in log units), with a higher ratio representing higher randomness.The GR for the intersaccadic intervals (GR-ISI) is plotted in Figure 1c per group and in Figure 1b per movie clip and group.The GR-ISI showed a significant movie effect (F(2, 427.23) = 10.22,p < 0.001), and a significant group effect (F (1, 217.49) = 33.93,p < 0.001); however, no significant interaction was found (F(2, 427.23) = 1.09, p = 0.34).A Tukey post hoc test showed significant pairwise differences between JB-JJ (t(430) = 4.44, p < 0.0001) and JB-Nat (t(432) = 2.96, p < 0.01); however, JJ-Nat (t (432) = À1.44,p = 0.32) did not show significant differences (see Table 2).In addition, when computing the Note: The GRs were calculated for the inter-saccadic intervals (GR-ISI), saccade amplitude (GR-SA), saccade durations (GR-SD), saccade peak velocity (GR-SPV), inter-blink intervals (GR-IBI) and blink duration (GR-BD).The percentage of good recording shows the results of the mean recording percentage between saccades (intervals) and in one movie clip (Movie).The additional section refers to the statistics under the condition of at least 80% recording in the intervals between the saccades, calculated for the intervals that fulfill the condition (%Rec.Intervals) and GR-ISI.The main effects and the pairwise comparisons, including the FDR analysis for multiple comparisons, were calculated as in Table 1.
Area-Under-the-Curve (AUC) of logistic regression (Figure 1e) using the GR-ISI as the predictor, there was a good discriminability of AUC = 0.79 for the observers' mean across movies, and a lower AUC, yet an acceptable value of 0.75 for the Jangle-book, an AUC of 0.74 for Jack-Jack, and an AUC of 0.69 for the Naturalistic movie clip.The saccade rate statistics were examined as a possible contributor to the GR-ISI significant results, yet no significant group and movie effects were found (see Table 1).In summary, the results showed higher randomness in the GR-ISI for the ASD group, with significantly higher randomness in the JB movie clip compared to the other clips.

Test-retest reliability of the inter-saccadic randomness measure
Pearson correlation coefficients were computed between the first and the last run GR-ISIs in the three movie clips: JB, JJ, and Nat, separated by group: ASD (a), and TD (b) (Figure 2).The correlations in the ASD group (Figure 2a) show; r(178) = 0.48 for JB, r(173) = 0.35 for Nat, and r(183) = 0.43 for JJ.The correlations in the TD group (Figure 2b) show; r(39) = 0.63 for JB, r(38) = 0.60 for Nat, and r(39) = 0.63 for JJ.All correlations show a significant positive correlation (p < 0.001), with stronger correlations in the TD group (Figure 2b).The significant correlations between the first and last runs may indicate an individual "randomness" measure attributed to each child.The test-retest results of the GR-BD showed lower correlations than the GR-ISI, yet significant in the ASD group in all three movie clips (JB: r(178) = 0.29, p = 0.001; JJ: r(183) = 0.26, p = 0.002; Nat: r(173) = 0.26, p = 0.003).However, the TD group did not show significant correlations (JB: r(39) = 0.29, p = 0.074; JJ: r(39) = 0.29, p = 0.074; Nat: r(38) = 0.32, p = 0.057) between the first and last run of the GR-BD results, presumably because of the smaller sample size.There was also a significant correlation between the GR-ISI measures across movie clips (each clip averaged across viewings) for the ASD children: r = 0.42, p < 0.001 for JB and JJ, r = 0.44, p < 0.001 for JB and Nat, and r = 0.47, p < 0.001 for JJ and Nat.The correlations in the TD group were lower and did not reach significance for JJ and JB.

Gamma distribution of blink duration
The blink duration distributions were fit with the gamma function; the results are shown in Figure 3. ASD children showed a more random blink duration distribution, with higher scale and lower shape values.The gamma ratio of the blink duration (GR-BD) is shown in Figure 3, with subplot c for the group means and subplot d for the group and movie means; the statistics are summarized in Table 2.We found a significant group effect (F (1, 213.33) = 43.52,p < 0.001) and a movie effect (F (2, 418.37) = 3.03, p = 0.049); however, the movie effect was weak and was found to be non-significant after FDR correction.In addition, the movie group interaction effect F I G U R E 2 Consistency of the GR-ISI measure across runs for the three clips.Each data point represents one observer's GR-ISI.Each plot shows the correlation between the first and last runs of each movie clip: Jungle Book, Jack-Jack Attack, and Naturalistic, in the ASD group (a), and the TD group (b).All correlations were significant (p < 0.0001), suggesting that the GR-ISI is an individual trait and a reproducible measure per movie.In addition, note that there is less variability in the TD group.
was found to be non-significant (F(2, 418.37) = 0.17, p = 0.85).The post hoc test did not show significant differences between movies (see Table 2).In addition, the AUC of the logistic regression using the GR-BD as a predictor showed an acceptable discrimination of 0.76 for the observers' mean across movies.All the basic blink properties showed significant, yet lower, group effects (Table 1), including the average blink duration (F (1, 222.76) = 19.43,p < 0.001), which can explain part of the gamma results obtained for the blink duration.In addition, unlike for GR-ISI, we did not find significant group, movie, or interaction effects for the GR calculated for the intervals between the blinks (onset to onset times), see the statistical results in Table 2.Note that the results of this analysis are based on a careful detection of blinks (see Section 2), and are unlikely to be based on arbitrary gaps in the recording of eye-tracking data.

Gamma distribution of basic saccade properties
We conducted gamma analyses of additional saccade parameters.First, we analyzed the saccade peak velocity gamma distribution, since it is the most natural adaptation of Torres et al. (2013) micro-movement analysis to eye movements.Contrary to expectations, the GR of the saccade peak velocity did not show a significant group effect (F(1, 222.45) = 2.34, p = 0.13, AUC = 0.57) (see also Table 2 and Figure S5a-c).However, gamma analysis (GR) of a related measure of saccade amplitude (known to be proportional to the peak velocity by the saccades main sequence) showed a small group effect (F (1, 218.24) = 4.36, p = 0.038, AUC = 0.58), which was found to be non-significant after FDR correction (see Figure S5d-f).In contrast, the gamma analysis (GR) of saccade duration did not show a significant difference between the groups (F(1, 220.36) = 3.20, p = 0.08, AUC = 0.57).All the GR measures of the basic saccade properties, that is, peak velocity, amplitude, and duration, showed significant movie effects, and the pairwise comparison showed significant results between JB-Nat and JJ-Nat.The results may indicate that the basic saccade characteristics of randomness were affected similarly by the movie content.

Gamma parameters and ASD severity
The gamma parameters of the inter-saccadic intervals and the blink durations were compared to the ADOS-2 comparison scores; the results are shown in Figure 4.The ADOS comparison scores show severity scores ranging from 1 to 10, where 10 represents the highest severity.All TD children were denoted as severity 0. The results for GR-ISI are shown in Figure 4a; the LM showed significant results with F(10, 190) = 7.67, p < 0.001, R 2 = 0.29.The results for GR-BD are shown in Figure 4b; the LM showed significant results with F(10, 189) = 4.54, p < 0.001, R 2 = 0.19.The statistics included comparison scores from 1 to 3, although some regard them as "nonspectrum" cases (Hus et al., 2014).When omitting 1-3 severity cases, including TD with low comparison scores, the LM results for GR-ISI showed F(7,177) = 10.59,p < 0.001, R 2 = 0.30, and the GR-BD showed F(7,176) = 5.88, p < 0.001, R 2 = 0.20.In summary, the gamma parameters' ratios for both inter-saccadic intervals and blink durations showed a significant positive correlation (denoted on each plot), and a significant LM with the ADOS comparison score, with higher randomness for higher severity.

Gamma parameters and age
We assessed whether age was significantly associated with the gamma parameters (Figure 5).Age was rounded to the nearest whole number, creating comparable age groups ranging from 1 to 7 years in both groups.Analysis revealed distinct age clusters within the ASD group (1-3 and 4-7 years), evident in the scale and shape plot (Figure 5a) and confirmed by K-means clustering, forming four clusters.In contrast, the TD group exhibited less pronounced age clusters and did not align with the same age division (1-3 clusters also contained age 6).In addition, LM was used to determine whether the age (1-7 years) and group could predict the GR-ISI (Figure 5b).The LM showed significant results with R 2 = 0.24, F (13,200) = 4.74, p < 0.001, along with a significant age effect (F (6,200) = 3.21, p = 0.005) and group effect (F (1,200) = 40.18,p < 0.001), without a significant interaction between group and age (F (6,200) = 0.37, p = 0.89).The post-hoc tests did not show significant pairwise differences between the age groups.In addition, statistics were checked after dividing the data into two age groups (age 1-3 years, and 4-7 years), derived from the observed decrease in GR-ISI after age 3, shown in Figure 5b, and the clustering data obtained from the ASD age groups in Figure 5a.The results yielded a higher age effect (F(1,210) = 20.72,p < 0.001), and a significant group effect (F (1,210) = 39.92,p < 0.001).The interaction between the group and the two age clusters was non-significant (F (1,210) = 0.69, p = 0.41).In conclusion, both age and group effects were found to be significant in predicting GR-ISI; the younger age was associated with higher GR values or increased "randomness" in both groups.

The effect of recording gaps on the gamma distribution analyses
In every head-free eye-tracking experiment, tracking could sometimes be lost due to movement that exceeds the tracking limits, producing missing data or gaps, with a percentage that could vary between individuals, groups, and conditions.More gaps could be expected in ASD children who often show excessive motor noise (Torres & Denisova, 2016).We first examined the quality of the recordings using two methods: (a) the percentage of recording per inter-saccadic interval in movies with more than 50% recording, and (b) the percentage of recording in one movie.The results showed that the recording quality in the ASD group was lower than in the TD group; both measures showed significant group effects of (a) F(1, 213.26) = 21.19,p < 0.001, and (b) F(1, 228) = 16.49,p < 0.001.These results, along with additional statistical information, are summarized in Table 2 (see Figure S2).We conducted an additional analysis to test the potential role of the excessive tracking gaps in ASD in producing the gamma parameter differences we found for the inter-saccadic intervals (GR-ISI).These intervals could have been affected by the gaps, that is, longer intervals produced by a lack of recording periods and missing saccades.For that purpose, we selected only intersaccadic intervals with more than 80% valid samples and a total recording above 50%.Using these thresholds, we obtained non-significant differences in the ISI recording percentage between groups (F(1, 222.80) = 2.96, p = 0.087); Thus, using a threshold of 80% recording of ISI equalized the average ISI gaps between the groups.
We then conducted the same analysis shown in Figure 1 only for ISIs with more than 80% recording between saccades (see Figure S3).The GR-ISI results under these conditions showed a lower AUC than the full data set (AUC = 0.65 and AUC = 0.79 respectively); however, the group effect (F(1, 213.19) = 9.00, p = 0.003) and movie effect (F(2, 420.47) = 9.29, p = 0.0001) showed significant results, with similar movie order mean and significant pairwise differences between JB-JJ and JB-Nat movie clips (see Table 2, Figures 1d, and S3d).These results imply that the GR-ISI effect is not due to excessive gaps in the ASD groups, since it persisted even when these gaps were equalized between the groups.

DISCUSSION
We assessed the "randomness" of oculomotor parameters in ASD using the Gamma fit parameters of their distribution.This method was inspired by the work of Torres et al. (2013), who used Gamma parameterization of body movement to obtain a "personalized signature of motor variability" shown for micro-movements such as the hands and head (Torres, 2011(Torres, , 2013;;Torres et al., 2013Torres et al., , 2016;;Torres & Denisova, 2016).In our main analysis, we assessed the distribution of the inter-saccadic intervals and blink durations; although we studied the time series distribution, and not the velocity distributions, we predicted similar results, where random oculomotor behavior is expected to generate nearly an exponential distribution, whereas a controlled process is expected to produce a close-to-normal distribution.Here we showed that ASD children have more temporal randomness in their eye movements while watching short video clips compared with TD children.Higher randomness corresponds to a higher gamma scale and a lower shape, producing a higher gamma ratio (log 10 of gamma scale divided by the gamma shape) in ASD children for both the inter-saccadic intervals and blink duration.This "temporal randomness" adds to the idiosyncratic viewing patterns in ASD (Avni et al., 2020) shown for the same movie clips, and it adds to the randomness of pointing kinesthetic (Torres et al., 2013) and micro-movements of the head (Torres & Denisova, 2016).

Basic saccades and blink properties
Overall, the basic saccade properties of rate, amplitude, duration, and peak velocity did not show a significant group effect, as previously found by Avni et al. (2020) in a subset of the same data set.On the other hand, all of the blink basic measures, including IBI, duration, and rate showed a significant group effect.Previous studies have shown that blinking and especially blink inhibition reflect visual processing and attentional engagement (Bonneh et al., 2016;Nakano et al., 2009;Ranti et al., 2020;Shin et al., 2015).The current finding of a higher blink rate suggests a reduced engagement with the video clips in the ASD group (Ranti et al., 2020).

Oculomotor randomness
There are many ways to assess oculomotor "randomness" while watching video clips, although it is not clear to what extent these measures are related to the content of the clip or represent a more fundamental property of gaze patterns.Our initial goal was to assess the oculomotor "randomness" in ASD, using the gamma parameters.We analyzed the inter-saccadic intervals as the time from one saccade onset to another (Figure 1) and the blink duration (Figure 3).The gamma analysis revealed a pattern similar to that of Torres et al. (2013), who tested the peak velocity of hand pointing and head movements (Torres et al., 2013;Torres & Denisova, 2016), showing a higher scale and a lower shape in the ASD group, suggesting a more random and noisier process in the ASD group.In our oculomotor results, we did not find significant differences between the groups using the GR of the saccades' peak velocity; however, the saccade peak velocity and the saccade amplitude are related by the "main sequence" (Bahill et al., 1975).Moreover, the GR of the saccades' amplitude showed a significant group effect (see Figure S5); however, it was a weak group effect, which was found to be non-significant following FDR correction.The different GR and gamma parameters' results of the saccade amplitude and the peak velocity could stem from the different units of measure (degrees vs. degrees per second), in addition to the variations that exist in the peak velocity of saccades with a similar amplitude (Leigh & Zee, 1999).
The recording quality in ASD and TD We considered the possibility that the results of higher randomness could be due to the children's movement during the experiment, which results in periods of a lack of good or valid eye-tracking recording.To determine whether the increased randomness in the saccade intervals was due to gaps in the recording, we checked the percentage of time of good recording (i.e., good or valid data) in both groups using two measures: (1) the percentage of good recording in the intervals between the saccades, and (2) the percentage of good recording in a full movie clip including permitted blinks.The results clearly showed that the recording quality in ASD was lower than in the TD group (Figure S2).However, when intervals with less than 80% recording were excluded, the quality of the remaining intervals was equalized between groups, but the significant difference in "randomness", as measured by the gamma parameters persisted (Figure S3), despite a reduction in its discriminatory power.This result rules out the possibility that the gamma parameter difference is merely a reflection of recording gaps.Note that eye blinks are a source of gaps in the recording and that their detection and dissociation from other recording gaps are not fully certain despite the restricted method for detecting the blinks (see Section 2).We found a significantly higher blink rate in the ASD group (Figure S5i, and statistical results in Table 1), even when excluding long periods with consecutive or long blinks, which could increase the blink rate.However, the reduced recording quality we found in ASD could not be attributed to the higher blink rate, since we included the blink periods when calculating the recording percentage.

Different movies show different oculomotor randomness
When comparing different video clips, the Jungle Book was found to have the highest ISI randomness and the Jack-Jack had the lowest in both groups (Figure 1, sub-plots b and d).Moreover, the ISI showed a similar randomness order comparing the three movie clips in both groups, suggesting that although all the movies contained social content, the oculomotor "randomness" was affected by the visual input (e.g., the number of figures in each frame and elements in the background) and the auditory input (e.g., musical, loud effects, or normal domestic soundtrack), with apparently additional internal randomness that is unique for each child (Figure 2).The results of a higher gamma ratio in JB could reflect, for example, the difficulty of ASD children to synchronize their rhythm with music (Dvir et al., 2020).

Gamma parameters and ASD severity
The GR, for both inter-saccadic intervals (ISI) and blink duration (BD), was found to be correlated with the ADOS-2 comparison score (Figure 4a,b), showing higher randomness for higher severity.This result indicates that more severe cases of autism are associated with higher randomness in their saccade timing and blink duration.Yet, despite the high significance of the correlation, there was much overlap in the GR values between severity levels.The different levels of randomness within the same comparison score group could stem from the known symptom variability within the same comparison score, including different levels of social-communication impairments and repetitive behaviors (Hus et al., 2014).Torres and Denisova (2016) showed similar results with micro-movements obtained from head movement during MRI recordings.They found a higher NSR (comparable to a higher gamma scale), and a more skewed distribution, characterized by low gamma shape values, with higher ASD severity.These results are comparable to the higher GR results in more severe ASD cases (Figure 4a,b).In addition, Torres and Denisova (2016) showed that higher social deficits measured by the ADOS Social Affect (SA) scores were associated with a higher NSR in comparison to the Restricted, Repetitive Behavior (RRB) domain scores (Gotham et al., 2009;Torres & Denisova, 2016).Here we found higher randomness with higher ASD severity scores; the comparison with other ADOS domains was left for future work.

Gamma parameters and the ASD cognitive scores
We found a negative correlation between the GR-ISI of children with ASD and their cognitive scores on the Mullen and WPPSI assessments (see Figure 4c-e).This negative correlation suggests that higher levels of randomness were associated with lower IQ scores.However, no significant correlation was found between GR-ISI and the Bayley cognitive scores.Our results are similar and consistent with those of Torres and Denisova (2016), who found that ASD individuals with higher verbal and performance IQ scores exhibited lower motor NSR (Noise-to-Signal Ratio) values and higher shape values, indicative of distributions closer to normal.In contrast, those with lower IQ scores displayed a noisier, memoryless exponential distribution.The consistency between the results of eye movements and micro-movements, as well as their relationship to IQ scores, suggests that oculomotor randomness, particularly in the form of ISI, can serve as a valuable measure for quantifying cognitive abilities.

Gamma parameters and age
The age-related findings unveiled a distinct pattern where children aged 1-3 exhibited higher values for GR-ISI and scale but lower values for shape when compared to children aged 4-7 in both groups (refer to Figure 5a,b).
Whereas the age groups within the ASD group showed evidence of clustering based on age, the TD group displayed a less distinct division into clusters.The variations observed in the gamma distribution outcomes across different age groups can be attributed to the developmental progress of oculomotor and visual functions in infants, as noted by Braddick and Atkinson (2011).Specifically, the reduction in saccade latencies with age (Bucci & Seassau, 2012;Ross et al., 1994;Yang et al., 2002) was taken into account in the GR-ISI analysis.However, despite the evolving oculomotor development in children, Torres et al. (2013) obtained analogous results, demonstrating higher NSR in younger child age groups when utilizing gamma distribution analysis to assess micromovements during hand pointing.In summary, the more distinct categorization of age groups within the ASD group, as opposed to the TD group, suggests that, in addition to the anticipated oculomotor development associated with age, the TD group gains better control over eye movements as they age, and their distribution approaches normality.In contrast, the ASD group's change in GR was less pronounced, and their distribution indicated more random and noisier eye movements.It is worth noting that additional factors such as gender did not yield significant differences between males and females (see Figure S6).

Oculomotor properties as a predictor of ASD
Previous studies have shown the feasibility of using classification models for the early detection of ASD by utilizing gaze features such as entropy, areas of interest, and fixation analysis (Gaspar et al., 2022;Kang et al., 2020;Mazumdar et al., 2021;Startsev & Dorr, 2019;Zhao et al., 2021;Zhao, Wei, et al., 2022).We have obtained preliminary results for predicting ASD in our sample using logistic regression along with the oculomotor parameters developed in this study.The parameters were derived from the gamma analysis of ISI and blink duration distributions (i.e., the scale, shape, and the ratio between them), as well as parameters that quantify the goodness of recording, and the basic characteristics of blinks.Owing to our unequal sample size between ASD and TD ($190 vs. $40), we provide the analysis here as a proof of concept.For this purpose, we performed bootstrapping by generating an additional 150 samples from the TD group data, resulting in a total sample size of n = 380.Using the leave-one-out methodology, the classification model applied to the enlarged sample showed an AUC value of 0.87, an accuracy of 85.5%, a false positive rate of 7.9%, and a true positive rate of 80.1%.These results suggest that mining eye movement data recorded during ordinary movie viewing, including parameters that quantify "oculomotor randomness", holds promise for the future development of an objective diagnostic approach for ASD.However, further research with larger and more balanced sample sizes is needed to validate and refine these results.

CONCLUSION
Inspired by studies of motor behavior in ASD, which found more "randomness" in simple pointing (Torres et al., 2013), we showed here that like with body movement, we can measure the level of control versus randomness merely by analyzing the time stamps of oculomotor events such as saccades and blinks.Our results show, using a large sample of children (n = 240), higher randomness in the ASD group, which allows good discrimination (AUC = 0.79 in ISI and AUC = 0.76 in BD) and shows significant correlations with the ASD severity.Importantly, the oculomotor randomness measure, unlike the gaze pattern measures, is not based on similarity to the typical average; however, it provides an independent assessment of a general individual property, which might not be specific to ASD, but it distinguishes ASD from TD.These findings could contribute to the future development of oculomotor biomarkers as part of an integrative diagnostic tool for ASD.
was obtained from all parents.For more information, see Avni et al. (2020Avni et al. ( , 2021)).

F
I G U R E 1 The inter-saccadic interval (ISI) gamma distribution parameters.(a) A scatter plot of the observers' mean across movies plotting the log units of the gamma scale and gamma shape, one point per observer.Note that the ASD group has a higher scale and lower shape values; the inner plot shows the groups' means.(b) The gamma parameters were averaged across movies and groups.Note that the randomness measurement was the highest in the Jungle Book movie clip in both groups.(c) Average GR-ISI (log gamma ratio) per group.(d) GR-ISI per movie, showing a significant group effect (p < 0.001) and movie effect ( p < 0.001).The post hoc test showed a significant difference between JB and Nat ( p < 0.01), and JB and JJ ( p < 0.0001), group-movie interaction was not significant.(e) The Area-Under-the-Curve (AUC) of the ROC discrimination function shows the best discriminability for the observers' mean across movies (AUC = 0.79) using the gamma parameters ratio (GR) of the inter-saccadic intervals as in (c,d).Error bars denote 1 SE of the mean across subjects.

F
I G U R E 3 The blink duration gamma distribution parameters.(a) A scatter plot of the observers' mean across movies, plotting the gamma scale versus the gamma shape (log units), one point per observer.The inner plot shows the groups' mean.Note the higher scale and lower shape values for the ASD group as found for the ISI.(b) Average GR-BD (log gamma ratio) per group.(c) GR-ISI per movie.Results show a significant group effect ( p < 0.001); however, the results did not show a significant movie effect or interaction.Error bars denote 1 SE of the mean across observers.

F
I G U R E 4 The gamma parameters and behavioral test scores; (a,b) ADOS-2 (c) Bayley, (d) Mullen, and (e) WPPSI.The ADOS-2 comparison score (a,b) ranges between 1 and 10, where 10 represents the highest severity.TD children were assigned a score of 0. As shown, there was a positive correlation between the ADOS comparison score and the gamma ratio (GR) of the (a) inter-saccadic intervals, and (b) the blink duration, showing a higher GR value in more severe cases.The average GR was calculated across movies with 1 point per observer; error bars denote 1 SE of the mean across observers.In (c-e), the results show significant negative correlations between GR-ISI and the cognitive scores of the Mullen and WPPSI tests, showing lower GR values in higher cognitive scores.The correlation results were denoted on each plot.Each dot in plot (c-e) represents the cognitive score of one ASD subject, averaged across repeated tests.

F
I G U R E 5 The gamma parameters and age.(a) Log gamma scale and the log gamma shape of the ISI according to age groups, one symbol per rounded age (year).(b) The gamma ratio of the inter-saccadic intervals (GR-ISI) in age groups from 1 to 7 years.(c) GR-ISI in 2 age groups: 1-3 and 4-7.In (a), the age clusters computed via k-means are marked in distinct colors on the markers' perimeter.The results show pronounced clustering in the ASD group.As shown in (b), the GR-ISI decreased with age and was significantly higher in the ASD group ( p < 0.001).Both (b) and (c) showed significant group effects (p < 0.001) and age effects ( p < 0.01 and p < 0.001, respectively).Error bars denote 1 SE of the mean across observers.
Basic oculomotor statistics for saccades and blinks.
T A B L E 1Note: The table shows the statistical results of Type III ANOVA calculated with Satterthwaite's method after fitting each oculomotor measure with a Linear mixed model, Gamma ratio (GR) statistics and the percent of good recording.