The participant was a 25-year-old, right-handed male, who has been playing the NBA2K basketball games for 5+ years, totaling over 4,000 h. The participant was seated in front of a 24'' monitor at roughly 15" distance from a Tobii 5L Eye Tracker (Tobii AB, Sweden). The monitor was connected to a PlayStation 4 (Sony Interactive Entertainment, San Mateo, CA) console. The video was captured through an Elgato HD60 S+ (Elgato Systems GmbH, Germany) capture card at native 60fps. Using corneal reflection (by illuminating the eye with infrared light and measuring the reflection from the cornea and the pupil), pupil diameter data was collected through Tobii 5L Eye Tracker sampling at 120 hz (For information on the accuracy and precision of Tobii 5L, please, see). EEG data were collected using a Neurosity Crown (Neurosity Inc., San Francisco, CA) EEG device (8 channels, 256 hz). The game was NBA2K 22 Basketball for the PlayStation 4, and the primary game mode was the Recreational Center, which is a 5-on-5 competitive game mode but is not ranked. The participant played 11 times in 1 day with breaks of 10-60 min in between. The recordings took place on 22 April 2022. To investigate the dynamic interactions over time, we chose three recordings (first, sixth, and eleventh), namely, Sessions 1, 2, and 3. Session 1 started at 14:51 (Pacific Time), Session 2 started at 19:12, and Session 3 started at 23:58. They all lasted about 30 min. Lighting was the same throughout and no stimulants were used.
EEG signals were recorded from four brain locations: P03, P04, CP3, and CP4 via the Neurosity Crown. That EEG device was connected via Wi-Fi and using the LabStreamingLayer (LSL) functionality provided by Neurosity. The LSL stream was connected to CuriaRecorder (Curia LLC, Dublin, CA), which is a fork of the open-source Labrecorder app provided by LSL. The CuriaRecorder is designed to subscribe to and record LSL streams, keeping timestamping and jitter information. The raw data was stored and preprocessed into native XDF format, after which it was run through basic de-noising (i.e., data >4SD of the dataset is removed). Finally, the raw and preprocessed data was made available and a final processing pass was performed. The recording was cut into 15-s intervals. Calculations were performed on these 15-s intervals as single data points. This technique allows long recordings to show the trends over time while being broken into discrete 'bins' that calculate the data independently.
For example, in order to collect data on eye movements, the utilized eye tracking device recorded at a sampling rate of 120 Hz. This resulted in a total of 1800 samples being collected within a 15-s processing window, representing the left pupil size in millimeters. The mean of these 1800 values was calculated and stored as the average for that 15-s window. Over the course of an 8-h data collection period, this design resulted in the storage of approximately 3.5 million values for the left pupil size alone. Given the large number of data points collected, and the potential for device noise and cumulative error over the prolonged data collection period, subsampling of the data was employed to "sparsify" the data and facilitate the identification of long-term trends while minimizing short-term variation in any given sample or bin. The subsampled data was extensively tested against the full sample data for both accuracy and various parameters, ranging from 5 to 60-s "chunks". It was determined that a 15-s processing window provided the desired benefits in processing without altering the shape of the data or the trends that emerged when compared to full sample recordings.
The Tobii 5L Eye Tracker was connected via USB3 to the collection system and was using Curia-written LSL broadcasters that take the Tobii datastream from the API and also push it to LSL as a datastream. Using the same acquisition method (i.e., the CuriaRecorder), the Tobii LSL streams were recorded with accurate timestamps and jitter information. They were also stored in LSL native XDF format. The data were preprocessed for outliers, stored as raw and preprocessed data, and finally cut into the same 15-s subintervals for comparison between time points.
Concerning extraction of outliers, the process for cleaning EEG data involves performing a 4-SD removal pass, followed by a 2 SD removal pass. The first step in the process is to calculate a rolling standard deviation of the data with a 1-s window. A threshold is then calculated by taking the 0.3 quantile of these values and multiplying it by 4, resulting in a threshold of 691 muV in this example. Any values (absolute values) exceeding this threshold, plus or minus 691 muV, are rejected and set to 0 for subsequent calculations. For an example channel with 380,000 samples, 866 samples were rejected for exceeding the threshold. The second pass is a standard z-score and 2 SD removal for each channel once the first pass is complete. For pupil data, an exponentially weighted smoothing algorithm is applied to the pupil and gaze values to compensate for the potential introduction of noise by the 120 Hz eye tracking device. This algorithm uses the pandas exponentially weighted mean (EWM) functionality, which is similar to a median filter. The smoothing factor used in this implementation is 0.3 (equivalent to moving the smoothing slider in "Tensorflow" for Python). The relevant code section is provided as Supplementary Material.
Of note, in the present study, we employed two distinct methods to address noise (i.e., any undesired or extraneous signal that is not generated by the brain) and outlier (i.e., data points that deviate significantly from the rest of the data and may indicate errors in the recording or data processing) removal from EEG data. The fundamental difference between the two methods is that the 4-SD removal pass was solely directed towards noise reduction in the raw data, while the combined usage of 4-SD and 2-SD removal passes were primarily employed for removing outliers from preprocessed data.
As stated above, the sampling frequency of the Tobii 5L eye tracker and Neurosity Crown EEG were 120 Hz and 256 Hz, respectively. In order to facilitate a comparison between the two datasets, it was necessary to establish a common sampling frequency. Utilizing the "scipy" resample function within the Python programming language, the eye tracking data, stored in a Pandas dataframe, was upsampled to 256 Hz to align with the EEG data. The decision to upsample the eye tracking data rather than downsample the EEG data was made based on the tradeoffs between frequency resolution and the potential loss of information. Therefore, the sampling frequency for both Neurosity Crown and Tobii 5L Eye Tracker was 256 Hz.
 Section 2.2-Section 2.5 summarize the main steps for the spectral decomposition and cross-correlations analysis (for a detailed explanation of the entire procedure, please see).
We first segmented the previously pre-processed signals (see Section 2.1. for details) from each brain location and eye into 2-s time windows with a 1-s overlap across each session. Within each 2-s time window, we extracted the spectral power S(f) from each EEG/EYE signal using the 'pwelch' function in Matlab, based on the discrete Fourier transform (DFT) and the Welch's overlapped segment averaging estimator. For each time window, we obtained a spectral power value in bins of 0.5 Hz (fbin) for the range from near 0 Hz to the Nyquist frequency, fN, of 128 Hz (i.e., the half of the sampling frequency of 256 Hz), which renders the total number of spectral estimates, N, for each window of 2 s N = fN * (1/fbin) = 256; that means that N = 256 is the number of spectral power data points for each window of 2 s. To probe specific contributions from different frequency bands fi to the spectral power within each 2-s time window of the EEG/Eye signal, we considered seven frequency bands corresponding to the commonly utilized cortical rhythms: delta (0.5-3.5 Hz), theta (4-7.5 Hz), alpha (8-11.5 Hz), sigma (12-15.5 Hz), beta (16-19.5 Hz), gamma1 (20-33.5 Hz), and gamma2 (34-98.5 Hz). We calculated the average power  across all frequency bins within each frequency: , where  are all frequencies considered in each frequency band. Thus, we obtained seven time series of EEG/Eye band power  with 1-s resolution for each Brain Location/Eye.
Based on the exploratory nature of this research (e.g., investigating the linear relationship between two time series data sets), for each protocol session and each pair of Brain location-Eye (a total of eight different Brain Location/Eye combinations), we calculated the bivariate equal-time Pearson's cross-correlation for all pairs of time series representing EEG/Eye spectral power in the frequency bands fi, where i = 1, ... ,7. This led to 49 (7 x 7) cross-correlation values Ci,j for each pair of Brain Location-Eye, as shown in the cortico-ocular cross-correlation matrices (Figure 1). The cross-correlation values range from Ci,j = -1 (fully anti-correlated) to Ci,j = 1 (fully correlated), with Ci,j = 0 indicating the absence of a linear relationship between the power spectra of two EEG/EYE frequency bands.
Cortico-ocular cross-correlation matrices represent pairwise coupling strength between the seven frequency bands fi of one brain location with the same bands derived from one Eye (i.e., eight distinct pairs: CP3-Eye Left, CP3-Eye Right, P03-Eye Left, P03-Eye Right, P04-Eye Left, P04-Eye Right, CP4-Eye Left, and CP4-Eye Right) during a given session (Figure 1A). Matrix elements indicate the coupling strength between the band-wise power of EEG signal in a given Brain Location and the pupil diameter dynamics in one Eye.
To visualize the information provided by the cortico-ocular cross-correlation matrices, we mapped the matrices in Figure 1 into different networks for Session 1, Session 2, and Session 3. Each brain location and the Eye is represented by a semicircle/circle, where color nodes represent distinct frequency bands fi. Network links correspond to the values of cross-correlation matrix elements Ci,j in Figure 1. Link strength is marked by line color and width. To illustrate the differences in network organization, we used the following link strength classification: positive links [Ci,j > 0.17; dark red); weak positive links (0.12 < Ci,j < 0.17; light red); negative links (Ci,j < -0.17; dark blue)] and weak negative links (-0.17 < Ci,j < -0.12; light blue). Links corresponding to cross-correlation values -0.1< Cij < 0.1 (Figure 2) are not shown in the network maps. Note that this set of thresholds was arbitrarily selected to provide the best visualization for the main results of the study shown in the cross-correlation matrices (Figure 1). To visualize the hierarchical organization of cortico-ocular network interactions among Brain Locations and Eyes, all multiplex networks (Figure 2) are presented separately as sub-network maps for all pairs of Brain Location-Eye (Figures 3, 4). Each sub-network map corresponds to a pair for a given Brain Location-Eye and follows the same color code as in the original network. Due to space limitations, we only present the sub-networks for the right hemisphere: P04-Eye Left, P04-Eye Right, CP4-Eye Left, and CP4-Eye Right (Figures 3, 4). The results for the left hemisphere (P03-Eye Left, and P03-Eye Right, CP3-Eye Left, CP3-Eye Right) are shown in Supplementary Figures S1, S2.
To examine the changes in the relationships between the frequency bins of each subnetwork pair across time, we used a multi-step process. It should first be noted that the correlation between sub-network pairs at each time point was based on 900 observations (i.e., one observation every 2 seconds for each 30-min session), which rendered statistical significance even for very small correlations. Thus, we examined the changes in standardized estimates of association (i.e., Cohen's q) across time. For all sub-networks examined, the correlation between frequency bin 1 for one node and frequency bin 7 for the second node of each pair changed considerably over time (see Supplementary Figure S3). For example, the correlation in the P04-Left Eye sub-network changed from 0.12 in Session 1 to -0.51 in Session 3.