Impaired saccadic eye movements in multiple sclerosis are related to altered functional connectivity of the oculomotor brain network

Highlights • Impaired eye movements in multiple sclerosis (MS) and functional connectivity (FC)• Eye movements related to altered FC of the oculomotor brain network.• Lower (beta band) and higher (theta/delta band) FC related to abnormal eye movements.• Regional changes were more informative than whole-network measures.• Eye movement parameters also related to disability and cognitive dysfunction.


Introduction
In many neurological diseases, the shift in focus from studying focal pathology to functional brain networks has increased our understanding of the etiology of clinical dysfunction, which is especially relevant in multiple sclerosis (MS). (Mollison et al., 2017;Schoonheim et al., 2014;Schoonheim et al., 2015) Functional brain network inefficiency resulting from accumulating structural damage varies between individuals, Abbreviations: MS, multiple sclerosis; FC, functional connectivity; MEG, magnetoencephalography; INO, internuclear ophthalmoplegia; EDSS, Expanded Disability Status Scale; Pv/Am, peak velocity divided by amplitude; BNA, Brainnetome Network Atlas; AECc, corrected Amplitude Envelope Correlation. but eventually leads to clinical and cognitive decline in MS.  Network dysfunction has mainly been studied with functional (f)MRI in MS, a measure of brain function based on levels of oxygenated blood. Magneto-encephalography (MEG) is a neurophysiological modality that directly measures neuronal activity with high temporal resolution. Recently, studies have shown that MEG is able to detect clinically relevant disruptions in functional networks in MS, potentially with higher sensitivity than fMRI. (Nauta et al., 2020;Tewarie et al., 2015;Sjøgård et al., 2021) These insights have underlined the importance of brain network changes in MS, but network evaluation in routine clinical practice remains impossible. Therefore, the search for easily measured and objective clinical outcome measures that can reflect abnormalities of the functional network is remains ongoing. (van Munster and Uitdehaag, 2017) Eye movement measurement is a strong candidate for such a marker of network dysfunction.
Eye movements are the fastest movements of the human body with highly consistent patterns of movement. (Leigh and Zee, 2015) Eye movement abnormalities are common is MS and clinically relevant due their disabling nature in daily life. Jasse et al., 2013) Recently, it has been shown that they relate to general disability, cognitive function and neurodegeneration in MS. Nij Bijvank et al., 2020;Fielding et al., 2015;Kincses et al., 2019) Recent methodological advances now enable a precise and non-invasive way of detecting (subclinical) eye movement deficits using highfrequency infrared oculography. (Leigh and Zee, 2015;Nij Bijvank et al., 2018;Frohman et al., 2003) With this technique the well-known pathology of internuclear ophthalmoplegia (INO) can be readily diagnosed. Frohman et al., 2003) In addition, precise quantification of more subtle abnormalities has become feasible. (Nij Bijvank et al., 2020;Sheehy et al., 2020; Importantly, eye movements are generated by a wide-spread and coordinated functional network. This network integrates sensory, motivational, executive and motor information, which is also extensively involved in cognitive function. (Leigh and Zee, 2015;Fielding et al., 2015;Coiner et al., 2019) Despite these theoretical arguments, the oculomotor network and the associations with eye movement abnormalities have not yet been studied in MS. We hypothesize that eye movements are more strongly related to functional connectivity of this oculomotor network than to whole brain functional connectivity. We aim to identify which eye movement measures are most indicative of this network and evaluate if the network as a whole or regional changes are most relevant. In addition, we verified whether and which eye movement measures also reflect clinical and cognitive function.

Study design and patient population
For this observational cross-sectional study, MS patients and healthy controls were included from the Amsterdam MS cohort, as previously described.  This study was approved by the Medical Ethical Committee on Human research of the Amsterdam UMC and followed the tenets of the Declaration of Helsinki. Written informed consent was obtained from all participants before study inclusion. Included patients were diagnosed with clinically definite MS according to the revised McDonald criteria, (Polman et al., 2011) and subdivided into relapsing-remitting, secondary progressive and primary progressive MS. (Lublin and Reingold, 1996) The Expanded Disability Status Scale (EDSS) score was used to assess the level of disability of MS patients. (Kurtzke, 1983) All assessments and data collection (clinical, infrared oculography, magnetic resonance imaging (MRI), magnetoencephalography (MEG) and neuropsychological evaluation), were performed on the same day and in the same order, as previously described.

Infrared oculography
Eye movements were measured and analysed using our standardised infrared video-oculography protocol, the DEMoNS protocol. (Nij Bijvank et al., 2018) In brief, eye movements were measured binocularly with the Eyelink 1000 Plus eye tracker at 1000 Hz during a pro-saccadic and an anti-saccadic task as described previously. (Nij Bijvank et al., 2020;Nij Bijvank et al., 2018) An off-line analysis was performed in Matlab (Mathworks, Inc., Natick, MA; version 2017b) to automatically calculate movement parameters, (Nij Bijvank et al., 2018) including latency, peak velocity, and gain of the centrifugal saccade. (Nij Bijvank et al., 2018;Antoniades et al., 2013) In the pro-saccadic task only, a main sequence relationship was calculated by dividing the peak velocity by the saccadic amplitude (Pv/Am). (Leigh and Zee, 2015;Nij Bijvank et al., 2020;Nij Bijvank et al., 2018) For the anti-saccadic task a set of three additional parameters was calculated, namely: (1) the proportion of errors (first saccade in the same direction as the target), (2) the latency of a correction saccade after an incorrect response and (3) the spatial error of the final eye position after the anti-saccade (before reappearance of the target). (Nij Bijvank et al., 2020;Nij Bijvank et al., 2018) Parameters were averaged over the left and right eye. Presence of internuclear ophthalmoplegia (INO) was determined using Versional Dysconjugacy Index based thresholds. . A higher value of latencies and errors, and a lower value of peak velocity and gain were interpreted as worse eye movement performance. (Nij Bijvank et al., 2020)

Magnetoencephalography (MEG)
Resting-state MEG measurements were pre-processed using a standardized procedure. (Nauta et al., 2020; In short, MEG data were acquired using a 306-channel whole head MEG system (Elekta Neuromag Oy, Helsinki, Finland), situated in a magnetically shielded room. Eyes-closed resting state measurements were performed (5 min) at 1250 Hz. MEG data were visually inspected to discard malfunctioning channels and the temporal extension of Signal Space Separation (tSSS) was used to remove artefacts. (Taulu and Simola, 2006) Source-localized MEG data were constructed for 224 cortical and thalamic regions (listed in Supplementary Table 1) of the Brainnetome Network Atlas (BNA) (Fan et al., 2016) using an atlasbased beamformer. (Hillebrand et al., 2012) An advantage of the BNA atlas is that it is created using a connectivity-based parcellation applying multimodal connectivity information. (Fan et al., 2016) For each subject the first 18 epochs of 16,384 samples (13.11 s) were used. The (pairwise) corrected Amplitude Envelope Correlation (AECc) (Brookes et al., 2011;Bruns et al., 2000) was used as a measure of functional connectivity (FC) in the delta (0.5-4 Hz), theta (4-8 Hz), alpha1 (8-10 Hz), alpha2 (10-13 Hz), beta (13-30 Hz) and gamma (13-48 Hz) frequency bands (see Supplementary File 1). The AECc was computed per epoch for each pair of BNA regions, and then averaged over the epochs. Finally, AECc values were averaged globally (whole brain connectivity) and also averaged over regions in the oculomotor network.

Definition of saccades and the oculomotor network
Saccades are the fast eye movements that change our line of sight and can broadly be divided in reflexive and volitional saccades. Pro-saccades are made from a fixation point towards, and in response to, the appearance of a visual target, also called reflexive saccades. In contrast, the anti-saccade is one type of volitional saccade, (Leigh and Zee, 2015;Antoniades et al., 2013) evoked by generating a saccade in the direction away from the target, while suppressing reflexive pro-saccades and correctly estimating the mirror location. (Leigh and Zee, 2015;Nij Bijvank et al., 2018) Previous work has shown that especially anti-saccades involve widespread cortical processing, which has led to its use in cognitive, neurological and psychiatric research domains. (Leigh and Zee, 2015;Anderson and MacAskill, 2013) Fig. 1 provides a schematic overview and description of the oculomotor network involved in the control of eye movements. This network was defined using the FOCuS atlas, (Coiner et al., 2019) and divided in three main sub-regions (see Fig. 2 and Supplementary Table 1).

Neuropsychological evaluation
All subjects underwent an extensive neuropsychological evaluation as previously described, (Eijlers et al., 2017) using an extended Rao's Brief Repeatable Battery of Neuropsychological tests (BRB-N) (Rao, 1990), including: Selective Reminding Test (verbal memory), Symbol Digit Modalities Test (information processing speed), Word List Generation Test (verbal fluency), 10/36 Spatial Recall Test (visuospatial memory), Concept Shifting Test (executive functioning), Stroop Colour Word test (attention), and Memory Comparison Test (working memory). Raw test scores were corrected for effects of age, sex, and level of education based on healthy controls using linear regression. Z-scores based on the distribution of healthy control scores were calculated for each individual domain and averaged across domains. (Eijlers et al., 2017;Amato et al., 2006)
Next, associations between network and eye movement data were explored using linear regression models. A stepwise approach was used to evaluate if and to what extent narrowing down from global connectivity to regional connectivity within the oculomotor network regions was meaningful and to select the most relevant variables for our final multivariate model (i.e. a feasible number of variables for our sample size), in summary (for more details, see Supplementary File 2): Fig. 1. Schematic overview of the oculomotor network, in which the regions are highly interconnected to integrate sensory, motivational executive and motor information which are required for eye movement execution. Three main sub-regions (occipito-temporal, parietal and (pre)frontal) were defined and the corresponding subareas are indicated with the arrows. Visual stimuli travel from the retina, via the lateral geniculate nucleus, to the primary visual cortex located in the occipital cortex. Preliminary processing takes places in, amongst others, the middle temporal cortex, which is involved in the perception of motion. These regions are connected with the posterior parietal cortex (parietal eye field and precuneus), which are responsible for constructing a spatial representation of the environment and directing spatial attention. From the posterior parietal cortex, direct connections with the superior colliculus exist, which can generate reflexive eye movements. For more volitional eye movements and broad-sale cognitive processing, information passes from the parietal to the (pre)frontal cortex. Regions in this area are thought to be involved in, amongst others, decisional and predictive processes, performance monitoring and motivation and modulation of motor commands. The frontal regions are directly connected to the superior colliculus and indirectly via the basal ganglia, the latter connections are related to evaluating the significance of an action. A premotor circuit in the brainstem, including the superior colliculus, is responsible for initiation and direct control of eye movement. Cerebellar regions are crucial for fine-tuning of the eye movement, and project to the brainstem. The cerebellum, along with premotor regions in the brainstem, also projects to the thalamus. The thalamus is a widely connected region which is responsible for directing visual attention and relaying the various information from other sources necessary for oculomotor control.
1. Global assessment, to identify associations between pro-and antisaccadic performance and FC values in the different bands (whole brain, ocular motor network and its three main sub-regions). Variables showing significant relations were entered into step 2. 2. Regional associations, in which selected associations were further explored by zooming into individual areas within the oculomotor network 3. Multivariate regression, to evaluate the effect of adjustment for possible confounders. Next, we investigated which saccadic parameters were most strongly related to the specific FC value, by combining these parameters in one model. In the main text adjusted standardized beta's are reported. Associations with a p-value lower than 0.01 were indicated in the result table (see Supplementary File 2 for more information).

Fig. 2.
Oculomotor network definition. Top: functional connectivity matrix of the healthy control group (left) and the MS group (right). The colours in the matrices represent the values of the AECc between two ROIs, as indicated in the colour bar. The larger red square indicates the ROIs that are included in the oculomotor network. The smaller squares represent the three main sub-regions of the oculomotor network, that are indicated with the same colours in the brain plots (bottom): (pre)frontal (blue), occipito-temporal areas (green), and parietal areas (darker red). The thalamus is not shown in the brain plots, but is included in the (pre)frontal areas (see also Supplementary Table 1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 4. Clinical relevance, as post-hoc analysis, we investigated if and which saccadic parameters resulting from step 3 were related to measures of clinical and cognitive functioning with linear and logistic regression models. Associations that survived the Holm-Bonferroni correction for multiple comparisons (Holm, 1979) were indicated in the result table. Furthermore, MS patients were divided in subgroups to evaluate the direction of results in comparison to healthy controls.

Results
In total 176 MS patients and 33 healthy controls were included in this study (see Table 1). Based on previous quality control Nij Bijvank et al., 2018) we excluded the pro-saccadic task of 8 subjects (7 MS, 1 control) and the anti-saccadic task of 4 MS patients. Neuropsychological data were available for 145 MS patients and 32 healthy controls.
Patients had a mean disease duration of 20.6 (±8.4) years, a median EDSS score of 3.5 (IQR 2.5) with a relapsing-remitting disease course for 65%. MS patients showed a significant lower average cognition than healthy controls (Z-score of − 1.07 ± 0.82 and 0.03 ± 0.57 respectively, p < .001. Altered saccadic performance has been published before, (Nij Bijvank et al., 2020) showing delayed latencies of both pro-and antisaccades in MS compared to controls (latency of 15 degrees prosaccades 207 ms and 225 ms respectively, p = .028), reduced gain of pro-saccades (0.95 and 0.92 respectively (non-INO only), p = .003) and an increased proportion of errors (0.50 and 0.37 respectively, p = .004). As previously published, abnormalities of latencies and proportion of errors were more pronounced in progressive MS patients compared to relapsing-remitting MS patients. (Nij Bijvank et al., 2020)

Step 1 and 2: Variables of interest
The screening step aimed to identify associations between saccadic parameters and FC as listed in Supplementary Table 2. The final associations and effects of adjustment for confounders are presented in Table 2, and the directions of the main results in Fig. 3A. For these associations, zooming in on the oculomotor network resulted in, on average, 13% higher effect sizes compared to whole brain associations. Moreover, zooming in further on main sub-regions resulted in, on average, 78% higher effect sizes compared to associations for the oculomotor network as a whole. For all but one of these associations, higher effect sizes were found for individual areas of the oculomotor network than for the larger region, with an average increase in effect size of 37%. Adjustment for confounders (age, sex and disease type) did not substantially change these effect sizes ( Table 2). The results are summarized per type of saccadic parameter in the following section.

Latencies
A higher latency of pro-saccades was associated with lower beta band FC for parietal and occipito-temporal areas, with the largest effect size for the precuneus (β = -0.21, p = .009, see Fig. 3B). In the antisaccadic task (correct responses), higher latency associated with lower FC for the beta band in the precuneus (β = -0.17, p = .030). The latency of a correction of an incorrect response related to higher delta band FC, with the largest effect in the primary visual cortex (β = 0.19, p = .012). The latency of incorrect responses was not associated with FC.

Errors
The proportion of errors in the anti-saccadic task was associated with higher FC (theta and delta band) of the parietal eye field, middle temporal cortex and the (pre)frontal sub-region (Table 2) and the largest effect was observed for theta band FC of the parietal eye field (β = 0.17, p = .028).

Peak velocity
Lower peak velocity and Pv/Am were associated with higher theta band FC, as well as lower and upper alpha band FC ( Table 2). The largest effect was found for theta band FC of the precuneus (adjusted model: β = -0.25, p = .007).

Gain and final eye position
Lower gain (hypometric saccades) of pro-saccades associated with higher theta and delta band FC of distributed areas ( Table 2). The largest effect was found for the association between the gain of 15 degrees prosaccades and theta band FC within the parietal eye field (adjusted model: β = -0.24, p = .004), see Fig. 3C. Gain and error of the final eye position of anti-saccades were not related to FC.

Sex effects
For latency of 8 degrees pro-saccades and proportion of errors in the anti-saccadic task, relevant effect modification by sex was found, explored further in Supplementary Table 3. In summary, latency of 8 degrees pro-saccades was negatively related to lower alpha band FC of parietal areas in female patients (adjusted model parietal eye field: β = -0.25, p = 0.010). In male patients higher gamma band FC of the inferior and superior frontal eye field and thalamus associated with a higher proportion of errors (adjusted model inferior frontal eye field: β = 0.29, p = 0.042).

Combining saccadic parameters in one multivariate model
To assess which (combination of) eye movement parameter(s) related strongest to the functional network, FC values (specific frequency bandregion combination) that were related to more than one saccadic parameter (Table 2) were investigated in an additional multivariate regression model. All relevant combined models involved theta band FC: Theta band FC of the middle temporal cortex was most strongly related to a combination of the proportion of errors in the anti-saccadic task (β = 0.18, p = .020) and the peak velocity of 15 degrees prosaccades (β = -0.20, p = .026). Theta band FC for the parietal eye field was related to a combination of peak velocity (β = -0.23, p = .012) and gain (β = -0.22, p = .009) of 15 degrees. Similarly, theta band FC of the precuneus and inferior frontal eye field was related to a combination of peak velocity and gain of 15 degrees pro-saccades (precuneus: β = -0.22, p = .014 and β = -0.19, p = .025, respectively; inferior frontal eye field: β = -0.21, p = .021 and β = -0.20, p = .015, respectively).

Step 4: Disability and cognition
Relations between clinical measures and saccadic parameters are shown in Table 3. In summary, in MS, a higher latency (of pro-saccades and correction in the anti-saccadic task) related to longer disease durations and higher EDSS. Furthermore, higher latency of 15 degrees saccades related to worse executive functioning (β = -0.49, p < .001), information processing speed (β = -0.37, p < .001) and working memory (β = -0.36, p < .001). Latency of 8 degrees saccades and latency of correct anti-saccades showed similar relations. Higher latency of a correction in the anti-saccadic task related to worse attention (β = -0.22, p = .015). More errors associated with a lower score in five domains, most strongly information processing speed (β = -0.36, p < .001) and attention (β = -0.31, p = .001). Peak velocity, Pv/Am and gain of prosaccades were not strongly related to clinical and cognitive function. Fig. 4 shows examples of abovementioned relations.

Subgroup analysis
Finally, subgroups were formed using controls as a reference, to study effects in patients with saccadic parameters that were abnormal (Table 2). This was only done for variables showing the strongest associations with FC, i.e. for latencies of pro-saccades (beta band) and gain of pro-saccades (theta/delta bands). The MS patients were therefore divided in 1) a normal and a high latency subgroup, 2) a normal and a low gain subgroup. In Fig. 5, FC values (beta and theta band) for MS subgroups and healthy controls are visualized for areas of the oculomotor network. Patients who displayed high latency and low gain generally showed a larger FC deviation from controls (i.e. lower beta band FC, higher theta band FC), than did the patients with normal saccadic parameters. The same pattern was observed for low gain and higher delta band FC. Taken into account the regions of Table 2, the FC was significantly different between the abnormal saccadic group and healthy control group for the parietal eye field (latency/beta band and gain /theta band), inferior frontal eye field (gain/theta band) and the middle temporal cortex (gain/delta band).

Discussion
This study showed that altered saccadic parameters in MS are related to FC of the oculomotor network. Regionally, latency was most strongly related to FC of parietal areas, errors, gain and peak velocity additionally to frontal FC. Saccadic parameters were also related to clinical outcomes, especially cognitive function.
Our study identified that abnormal eye movements in MS are mostly related to higher delta and theta, but lower beta FC, which was not studied before. Most previous studies either used other modalities, other measures of network function (power, coherence or different FC measures), or investigate task-related responses, which makes direct comparison of our results not feasible. Signals in the delta band are often attributed to artefacts, typically in relation to blinks and eye movement (Bodala et al., 2016), or are linked to sleep and deep levels of relaxation. (Mandal et al., 2018) However, increased delta and theta band power and FC can also reflect pathologic conditions, amongst others shown in  Table 2. sFEF: superior frontal eye field; iFEF: inferior frontal eye field; PEF: parietal eye field; MT: middle temporal cortex; PCUN: precuneus; PV: primary visual cortex. B) and C) Exemplar scatterplots of associations between a saccadic parameter and FC (AECc) . The linear fit of the unadjusted association shown (solid line), with the corresponding unstandardized regression coefficient and 95% confidence interval (dashed lines). B) Association between latency of 15 degrees pro-saccades and beta band FC of the precuneus. C) Association between the gain of 15 degrees pro-saccades and theta band FC of the parietal eye field. malignancies of the brain and MS. (DEJONGH et al., 2003;Tewarie et al., 2014) Beta waves are dominant in normal conscious states and during attentive cognitive tasks, and theta band rhythms are associated with drowsiness or sleep, as well as memory and learning. (Mandal et al., 2018) In general, the pattern of lower FC in the beta band and higher FC in theta/delta band could indicate generalized slowing of the network, which can be considered as an indicator of neurodegenerative processes, such as in Lewy Bodies dementia, but has also been described in MS. (Tewarie et al., 2014;Dauwan et al., 2018) This generalized slowing has also been related to cognitive dysfunction in MS, (Schoonhoven et al., 2018) to which our results also add relevance for eye movement abnormalities.
Regarding the specific localization of abnormalities within the network, our results showed that a higher latency of pro-saccades (a longer reaction time of reflexive movements) related to a lower beta band FC of the precuneus and parietal eye field (PEF). The PEF is located in the posterior intraparietal sulcus (IPS). Previous work has shown that the PEF plays a role in attention and visuospatial integration in conjunction with other regions in the IPS and is crucial for generation of reflexive saccades. (Coiner et al., 2019) The precuneus has also been implicated in eye movement control. (Leigh and Zee, 2015;Coiner et al., 2019) Although this region is not specific for oculomotor function, it has an important role in overall network functioning as part of the defaultmode network, as well as for attention and visuospatial processing. (Coiner et al., 2019) A higher latency of pro-saccades can potentially result from a delay anywhere in the afferent or efferent visual system, for which we now shown involvement of these two parietal areas in MS.
This study supports the hypothesis that abnormal saccadic parameters can reflect clinically relevant dysfunction in the functional brain network. Especially, our results confirmed that an altered latency of saccades might reflect a broad range of cognitive functions, including processing of visual information, task planning, attention and selection of relevant stimuli. (Leigh and Zee, 2015) This might indicate that latency partly reflects the reaction time / speed component that contributes to the performance on different cognitive tests as well. This could indicate involvement of the same underpinning network, which should be studied further. Peak velocity and gain of pro-saccades showed strong relations with theta band FC in the regression models, but were less consistently related to clinical variables and cognitive function. This suggests that these parameters reflect network dysfunction that is not directly related to other (clinical) outcomes.
The anti-saccadic task of our study represents volitional eye movements and theoretically requires more broad scale cognitive processing and involvement of frontal regions. Although associations of antisaccadic parameters with FC were weaker than for pro-saccadic parameters, our results do suggest that the proportion of errors in this task was related to a more widespread and less region-specific increase in theta band FC of MS patients. We did not find a significant relation Bold p-values (p) represents values lower than 0.01. Unstandardized regression coefficients (B) are all multiplied by a factor of 10,000 and are presented per 10 ms for latency, per 10 degrees/second for peak velocity, per 1 degree/second/degree for Pv/Am per 1 degree for the error of the final eye position and per 0.1 for gain and proportion of errors. Model 1: raw association; Model 2: association adjusted for age and sex; Model 3: association adjusted for age, sex and disease type. In model 3 the 95% confidence interval (CI) of B and the standardized regression coefficient (β) are additionally listed. In all models, parameters that are directly influenced by internuclear ophthalmoplegia (peak velocity, Pv/Am and gain) are additionally adjusted for the presence of unilateral or bilateral internuclear ophthalmoplegia. FC: functional connectivity; PS: pro-saccades; AS: anti-saccades; 15: saccades made in response to target amplitude of 15 degrees of visual angle; 8: saccades made in response to target amplitude of 8 degrees of visual angle; Pv/Am: peak velocity divided by amplitude; CI: confidence interval.
specifically with FC of the dorsolateral prefrontal cortex though, a region traditionally associated with the number of errors in the antisaccadic task. (Leigh and Zee, 2015) The proportion of errors did relate strongly to cognitive function in multiple domains. The relation with network dysfunction was potentially not fully captured by the connectivity measure chosen in this study, which is a conservative measure that avoids spurious connections that are due to the effects of volume conduction/field spread at the cost of ignoring true zero-lag functional interactions. We also observed sex effects (Supplementary Table 3), which might even have been underestimated due to the stepwise approach in the analyses; we only tested effect modification in the final model. Nevertheless, our findings underline that sex differences are important to take into account in MS studies, especially in network studies, because there is evidence that both the alterations to functional wiring of the brain and the severity of MS pathology are especially severe in the male MS brain. (Schoonheim et al., 2012) Some potential limitations of our study have to be considered. MEG generally has a lower spatial resolution than functional MRI (fMRI), although its spatial resolution has recently improved in newer models. (Hillebrand et al., 2012;Gross, 2019) It is also a less commonly used modality in MS, which makes it more difficult to compare our results with previous studies. Although recent studies have demonstrated that the sensor-space data of MEG can be accurately projected to deeper brain regions,  whether the same also holds for the brainstem or cerebellum remains unclear. These structures are relevant for accurate eye movements, therefore associations with pathology in these structures would be relevant to study in more detail, e.g. relations between FC, lesion volume and measures of atrophy The higher temporal resolution of MEG is an advantage compared to fMRI, as it allows for extension of our approach to more sophisticated parameters, such as high-resolution dynamic FC, (Tewarie et al., 2019) in future studies. Another potential limitation is the definition of the oculomotor network. Exact locations of areas involved in eye movement control have not always been clearly delineated, (Leigh and Zee, 2015) can vary between individuals, or are not accurately represented in different atlases. Nevertheless, we did find associations for the defined areas, most frequently stronger than for the whole brain or the oculomotor network as a whole. However, we did not investigate specific regions outside this oculomotor network. A next step could be to use more advanced statistical approaches, for example deep-learning algorithms, which could potentially select the most relevant regions and interactions between regions that are involved in oculomotor dysfunction, in a nonhypothesis driven manner. Another point to consider is that we included MS patients with a relatively long disease duration and as part of an ongoing cohort. We cannot draw conclusions on the occurrence of eye movement abnormalities and the relation to network dysfunction in patients during early stages of the disease. A methodological challenge for the study was the number of potential statistical comparisons, which were reduced by our stepwise approach in the statistical analysis. This was specifically used to show that zooming in towards regional measures yields stronger clinical relations, which should be validated in another patient sample. Finally, the cross-sectional design of our study limits the ability to draw conclusions about the order in which brain and eye movement abnormalities occur, and their causal relationships. Mediation analyses in a longitudinal dataset could help to crystallize the expected causal chain of (1) structural damage, (2) brain network dysfunction and (3) eye movement abnormalities (as well as other clinical outcomes).
To conclude, this study indicates that eye movement disorders are related to functional network changes of the oculomotor network and also strongly relate to cognitive impairment in MS. Latency, gain and peak velocity of pro-saccades were most strongly related to FC and regional oculomotor network changes were more relevant than global network changes in this MEG study. To further elucidate the relation with network dysfunction, future work should focus on more complex network measures, longitudinal associations and also include FC of the brainstem and cerebellum.  Geurts is an editor of MS journal and serves on the editorial boards of Neurology and Frontiers of Neurology and is president of the Netherlands organization for health research and innovation andhas served as a consultant for Merck-Serono, Biogen, Novartis, Genzyme and Teva Pharmaceuticals. B.M.J. Uitdehaag has received consultancy fees from Biogen Idec, Genzyme, Merck Serono, Novartis, Roche and Teva. L.J. van Rijn reports no disclosures. A. Petzold reports personal fees from Novartis, Heidelberg Engineering, Zeiss, grants from Novartis,outside the submitted work; and is part of the steering committee of the OCTiMS study which issponsored by Novartis and the Angio-OCT steering committee which is sponsored by Zeiss. He doesnot receive compensation for these activities -M.M. Schoonheim serves on the editorial board of Frontiers of Neurology, receives research supportfrom the Dutch MS Research Foundation, and has received compensation for consulting services orspeaker honoraria from ExceMed, Genzyme and Biogen. Fig. 4. Exemplar scatterplots of associations between saccadic parameters and cognitive domain Z-scores of MS patients. The linear fit of the unadjusted association is shown, with the corresponding unstandardized regression coefficients and 95% confidence interval (dashed lines). A) Association between latency of 8 degrees prosaccades and the executive functioning Z-score. B) Association between the latency of correct responses in the anti-saccadic task and the attention Z-score. C) Association between the proportion of errors in the anti-saccadic task and the attention Z-score. D) Association between the proportion of errors in the anti-saccadic task and the information processing Z-score.