Modulation of Gamma Spectral Amplitude and Connectivity During Reaching Predicts Peak Velocity and Movement Duration

Modulation of gamma oscillations recorded from the human motor cortex and basal ganglia appears to play a key role in movement execution. However, there are still major questions to be answered about the specific role of cortical gamma activity in both the planning and execution of movement features such as the scaling of peak velocity and movement time. In this study, we characterized movement-related gamma oscillatory dynamics and its relationship with kinematic parameters based on 256-channels EEG recordings in 64 healthy subjects while performing fast and uncorrected reaching movements to targets located at three distances. In keeping with previous studies, we found that movement-related gamma synchronization occurred during movement execution. As a new finding, we showed that gamma synchronization occurred also before movement onset, with planning and execution phases involving different gamma peak frequencies and topographies. Importantly, the amplitude of gamma synchronization in both planning and execution increased with target distance and predicted peak velocity and movement time. Additional analysis of phase coherence revealed a gamma-coordinated long-range network involving occipital, frontal and central regions during movement execution that was positively related to kinematic features. This is the first evidence in humans supporting the notion that gamma synchronization amplitude and phase coherence pattern can reliably predict peak velocity amplitude and movement time. Therefore, these findings suggest that cortical gamma oscillations have a crucial role for the selection, implementation and control of the appropriate kinematic parameters of goal-directed reaching movements.


INTRODUCTION
Voluntary movements are accompanied by modulation of oscillatory activity of both the beta (13.5-25 Hz) and gamma frequency ranges (25.5-90 Hz) that have been consistently observed with EEG, MEG and ECoG over the sensorimotor cortex and with electrode recordings in basal ganglia structures (Kilavik et al., 2013;Nowak et al., 2018). Movement-related beta oscillatory event-related desynchronization (ERD) and synchronization (ERS) have been characterized to a great extent in both healthy and clinical populations (Kilavik et al., 2013;Little and Brown, 2014;Cao and Hu, 2016; Barone and Rossiter, 2021). On the other hand, the current knowledge on the functional significance of gamma ERS is more limited and mostly stems from a small number of studies, some with ECoG and intracranial recording in small samples of patients with various neurological problems (Crone et al., 1998;Pfurtscheller et al., 2003;Androulidakis et al., 2007;Miller et al., 2007;Ball et al., 2008;Brücke et al., 2008Brücke et al., , 2012Cheyne et al., 2008;Tan et al., 2013aTan et al., ,b, 2016Flint et al., 2014;Miocinovic et al., 2015;Rowland et al., 2015;Gunduz et al., 2016;Ryun et al., 2017;Lofredi et al., 2018;Jiang et al., 2020), and a very few others with EEG or MEG in healthy subjects (Muthukumaraswamy, 2010;Gaetz et al., 2013;Amo et al., 2016).
For instance, the results of some ECoG studies of gamma ERS in the human primary somatosensory cortex (Miller et al., 2007;Avanzini et al., 2016;Ryun et al., 2017) have suggested an association between gamma oscillations and the proprioceptive feedback ensuing movement and sensorimotor input. However, this hypothesis has been challenged by reports showing that passive movements do not elicit gamma synchronization (Muthukumaraswamy, 2010;Brücke et al., 2012) and that mirror illusion in absence of proprioceptive feedback prompts movement-related gamma ERS (Butorina et al., 2014). Together, these results led into considering movement-related gamma ERS as directly functional to the active control of the ongoing movement. Indeed, differently from the beta ERD/ERS dynamics, gamma activity seems to display a more direct association with some characteristics of the motor output. In fact, sensorimotor and subcortical gamma ERS is, in general, time-locked to movement onset (Crone et al., 1998;Ball et al., 2008;Cheyne et al., 2008;Muthukumaraswamy, 2010;Brücke et al., 2012;Joundi et al., 2012;Cheyne and Ferrari, 2013;Gunduz et al., 2016;Lofredi et al., 2018). Also, the amplitude of gamma ERS during movement is variably linked to some of the movement features. In particular, links between gamma ERS amplitude and the amount of force generated have been described by EcoG studies of the motor cortex of epileptic patients (Flint et al., 2014;Jiang et al., 2020), by recordings in the subthalamic nucleus of patients with Parkinson's disease (Tan et al., 2013a(Tan et al., , 2016 and by a MEG study over the contralateral sensorimotor region of eight normal subjects (Muthukumaraswamy, 2010). Also, some relationships between gamma power and movement amplitude and velocity have been found with recordings in the globus pallidus of dystonic patients (Brücke et al., 2012) and in the subthalamic nucleus of patients with Parkinson's disease (Lofredi et al., 2018), respectively.
The notion of a prokinetic nature of gamma oscillations is also supported by correlations of abnormal gamma oscillations with motor symptoms of Parkinson's disease, essential tremor and dystonia (Rowland et al., 2015;Swann et al., 2016;Nowak et al., 2018;Guerra et al., 2020Guerra et al., , 2021. Although gamma ERS is prominent during movement execution and linked to motor symptoms, it worth noting that a few studies reported increased gamma activity even before movement onset (Brücke et al., 2008;Gaetz et al., 2013;Gunduz et al., 2016;Jiang et al., 2020), thus suggesting that gamma oscillations might represent a signature of processes linked to goal-directed movement representation, planning and execution.
In recent studies, we have extensively used a reaching task that allows for the parameterization of peak velocity and acceleration. In fact, in that task, the fast and uncorrected movements to different target distances result from a more prominent scaling of peak velocity rather than of movement duration (Tatti et al., 2019(Tatti et al., , 2020(Tatti et al., , 2021. Using that task, we have demonstrated that the movement-related beta ERD-ERS dynamics does not depend on target distance, movement length or peak velocity (Tatti et al., 2019). In the present study, we parametrically explored the relationship between movement distance and peak velocity on movement-related gamma ERS in a large group of healthy subjects. Specifically, we wished to investigate whether the amplitude of movement-related gamma ERS scales with movement distance and peak velocity. We thus recorded highdensity EEG activity in 64 healthy young subjects performing reaching movements toward targets at three different distances and characterized changes of gamma oscillatory amplitude and phase-coherence activity. According to the scanty and fragmentary evidence on the pro-kinetic nature of movementrelated gamma ERS, we expected that, compared to shorter movements, longer movements would be associated not only with higher peak velocities but also with greater gamma ERS (Muthukumaraswamy, 2010;Brücke et al., 2012;Lofredi et al., 2018). We also verified whether the scaling of gamma ERS would occur to the same extent and topography during movement planning and execution. Further, we tested whether the amplitude and phase-synchronization of gamma oscillations could reliably predict such kinematic parameters.

Subjects and Experimental Design
We enrolled 64 right-handed healthy subjects (mean 23.9 ± 4.5 years, 39 women) with normal or corrected vision, and without known disorders of the nervous system. This investigation was approved by the CUNY University Integrated Institutional Review Board (UI-IRB) and performed in accordance with the ethical principles of the Declaration of Helsinki and its subsequent amendments. Each participant signed an IRB-approved informed consent form before completing the experiment.
High-density EEG was recorded with the 256-channel HydroCel Geodesic Sensor Net (Electrical Geodesics Inc., Eugene, OR) while the participants performed planar reaching movements toward target placed at three distances (mov test).
Participants performed 96 out-and-back reaching movements with their right dominant hand by moving a cursor on a digitizing tablet to targets appearing on a computer monitor. The targets were circles randomly presented every 3 s at three distances (4, 7, and 10 cm; radius: 0.5, 1, 1.25 cm, respectively) in eight directions (45 • separation) ( Figure 1A). The central starting point and the cursor position were always visible. Instructions were: to move the cursor as soon as possible without corrections, but only after the target presentation, with movements as fast and as accurate as possible, with overlapping strokes with fast reversals in the target circle.

Kinematic Data Recording and Analyses
The (x,y) coordinates of each trajectory were recorded with a custom-designed software and analyzed using an ad hoc MATLAB-based pipeline. First, we filtered the coordinates with a Butterworth filter and then computed the first, second, and third derivative of the trajectory to obtain velocity, acceleration, and jerk for all the movements.
As detailed in previous publications (Ghilardi et al., 2000;Perfetti et al., 2011;Nelson et al., 2017), we computed several measures for each movement. In this study, we focused on: reaction time (i.e., the time from target appearance to movement onset), movement time (i.e., the duration of the outgoing movement), total movement time (i.e., the duration of the out and back movement), amplitude of peak velocity of the out-going segment ( Figure 1A). Movements with any of such measures outside two standard deviations and those rejected from EEG preprocessing were excluded from EEG analyses.
Following Kolmogorov-Smirnov and Shapiro-Wilk normality tests on standardized residuals of each kinematic parameter, two participants were excluded from the analyses. Therefore, analyses of kinematic data were performed on the movements of 62 participants. Repeated measures ANOVA were run on reaction time, movement time, peak velocity amplitude and timing, and movement extent, with target distance (short, medium and long) as within-subjects factor. Violations of sphericity assumption were Greenhouse-Geisser-corrected and significant main effects (p < 0.05) were followed by Bonferroni-corrected pairwise comparisons. Further, we computed linear mixed-effect regression models using the MATLAB function fitlme to unveil the specific contribution of peak velocity and movement time on movement extent. Linear mixed-effects regression analysis is a versatile extension of simple linear regression models that allows the estimate of both fixed and random effects, thus controlling some expected variation on the independent variable (e.g., intersubject variability). Importantly, as mixed-effect models fit an intercept and/or a slope for each random-effect, they address the problem of the non-independency of the data (i.e., the inclusion of multiple trials per subject), while respecting betweensubject variability.
For each subject, we included all the available trials and set peak velocity and movement time as fixed-effect factors; the 62 participants were instead included as a random-effect factor. To best account for between-subjects variability, we tested the fit of a model including either an individual intercept or intercept and slope for each participant. To ascertain whether and which fitted model provided the best fit to the data, the following metrics were assessed: adjusted R 2 , Bayesian Information Criterion (BIC, an index used in Bayesian statistics to select among two or more models), and the Theoretical Likelihood Ratio Test (TLRT), which is commonly used to compare the goodness of fit of two statistical models. Specifically, the TLRT compared the model with random intercept and slope and the one with random intercept by computing their likelihood ratio test under the Chi-square distribution.
Importantly, even if visual inspection of residual plots did not reveal heteroscedasticity and deviations from normality, and Kolmogorov-Smirnov and Shapiro-Wilk tests on the residuals were not significant, all the included values were z-transformed to obtain standardized estimates from the regression models.

EEG Recording and Analyses
High-density (HD) EEG data were acquired using a 256-channel HydroCel Geodesic Sensor Net (Electrical Geodesic Inc.) with a Net Amp 300 amplifier (250 Hz sampling rate, online reference electrode: Cz) and Net Station software (version 5.0). Sampling frequency was 250 Hz and channel impedances were maintained below 50 k throughout the recording to preserve a good signalto-noise ratio.
All recorded data were preprocessed using the public Matlab toolbox EEGLAB version 13.6.5b (v.2016b) (Delorme and Makeig, 2004). The continuous signal was first filtered using a Finite Impulse Response Filter (FIR) between 1 and 80 Hz and Notch filtered at 60 Hz (59-61 Hz). Then, the signal was divided in 4-s epochs centered on target onset (-1 to 3 s) and visually inspected to remove sporadic artifacts and channels with poor signal quality.
Independent Component Analysis (ICA) with Principal Component Analysis (PCA)-based dimension reduction (max 108 components) was run to delete stereotypical artifacts, such as eye blinks, muscular activity and heartbeat. After a visual inspection of the power spectral density, topographical maps and time course of each estimated component, we retained an average of 16.05 ± 6.27 components per subject. Channels previously removed due to bad signal quality were reconstructed using spherical spline interpolation, whereas those located on the cheeks and neck were removed. Re-reference to overall signal average was finally applied on the resulting 180 channels.
All the subsequent analyses were carried out using custom data analysis scripts with the MATLAB-based Fieldtrip Toolbox .
Importantly, to avoid ambiguous effects from improperly executed movements, after the preprocessing, we discarded epochs representing movements whose kinematic parameters exceeded two SD.
After trial rejection, the average number of trials per subject was 76.2 ± 13.2 SD, with a similar number of trials for each target distance (short: 25.6 ± 4.6, medium: 26.2 ± 5.9, long: 24.4 ± 5.2). Data were then time-locked to movement onset (-1 to 2.5 s). Time-frequency representations (1-80 Hz) were computed by convolving the signal using Complex Morlet Wavelets at linearly spaced frequencies (1-80 Hz, 0.5 Hz bins steps) and a constant time-window (-1 to 2.5 s). The number of wavelet cycles and length were increased as a function of frequency (cycles 3-10, 3.14-0.11 s). We applied this approach to obtain a good balance between frequency and temporal resolution (i.e., lower wavelet width will increase temporal resolution and reduce frequency resolution and vice versa).
Each trial was baseline corrected by subtracting and dividing the average signal of the entire time-window of all trials. In keeping with our previous works (Tatti et al., 2019(Tatti et al., , 2020(Tatti et al., , 2021 on movement-related beta oscillations (13.5-25 Hz), here we defined gamma oscillatory activity as ranging from 25.5 to 80 Hz.

Movement-Related Gamma Oscillatory Dynamics
EEG spectral and time-frequency analyses were run using the non-parametric cluster-based permutation procedure implemented in Fieldtrip (Maris and Oostenveld, 2007). Briefly, t or F statistics was first computed for each data point using a critical alpha of 0.001 and a minimum number of four significant neighboring electrodes to form a cluster. Cluster-level statistics was then run using the sum of the t/F values within each cluster of electrodes; the largest statistic value from the cluster-level analysis was then compared with a distribution of maximum cluster values obtained with 10,000 permutations (Monte Carlo method, alpha = 0.0005).
In order to characterize gamma oscillatory activity (25.5-80 Hz) during the planning and execution phases of the reaching movements, we first ran non-parametric permutation statistics and compared the broad gamma band time-course activity with the average spectral power (paired t-statistics, time window: -500 to 2,400 ms in 24 ms-time bins). Significant consecutive timewindows were then used to run second-level analyses against the average spectral power to summarize the topography of different sub-bands of gamma (gamma: 25.5-80 Hz, low gamma: 25.5-40 Hz, medium gamma: 40.5-55 Hz, high gamma: 55.5-80 Hz). With the resulting clusters of electrodes (ROIs), timewindows and gamma frequency bands, we then explored whether target distance affected gamma amplitude. Thus, cluster-based permutation statistics were run with target-distance as withinsubjects factor (repeated-measure ANOVA, planning timewindow: -152 to -52 ms; movement execution time-window: 52-500 ms), along with the correspondent permutation-based post-hoc tests.

Movement-Related Gamma Peak Frequency Analysis
To investigate whether the different phases of the reaching movements would be also characterized by spectral differences within the gamma range, we extracted for each subject the peak frequency value during the two significant time windows (planning and execution, see above) and a post-movement control time-window (1,500-2,000 ms after movement onset), when movements were certainly completed. Due to the distribution of the independent variable, normality assumption could not be satisfied (Shapiro-Wilk test on both data and residuals resulted in p < 0.05). Therefore, non-parametric related-samples Friedman's test was run to check for peak frequency differences in the three time-windows. Post-hoc pairwise comparisons were obtained with Dunn's test and Bonferroni correction for multiple-comparisons (alpha = 0.05).

EEG Functional Connectivity Analysis
To characterize gamma functional connectivity during movement planning, execution and post-movement control window, we computed the squared weighted phase lag index (wPLI), a measure of phase-synchronization implemented in Fieldtrip as the debiased wPLI (Vinck et al., 2011).
The wPLI derives from Stam's Phase Lag Index (Stam et al., 2007), as it introduces a phase-difference weighting normalization with the imaginary component of the crossspectrum (Nolte et al., 2004), thus improving robustness to noise. Also, the wPLI has the advantage of not being spuriously affected by the volume conduction of independent sources to different sensors or by a common reference, and shows increased statistical power to highlight true changes in phase-synchronization. The debiased WPLI estimator is obtained by computing (1) the spectrum and cross-spectrum of a signal, (2) weighting the crossspectrum with the average imaginary component of the crossspectrum. Thus, the debiased WPLI can summarized with the following formula, as detailed in Vinck et al. (2011): where E denotes the expected value operator and I{X} the imaginary component of the cross-spectrum. Therefore, we first computed the power and cross-spectrum during the planning and execution time-windows using Complex Morlet wavelets (0.5 Hz bins, width = 7) for two gamma frequency-bands (all gamma: 25.5-80 Hz, high-gamma: 55.5-80 Hz) and then estimated the debiased wPLI across all channelpairs. The result of this process was a weighted network (for each time-window and each frequency band) represented as a 180 x 180 adjacency matrix C = [c ij ], where each node is represented by a given electrode and each edge as the node-wise functional connectivity estimated by the wPLI. The functional connectivity matrices (Supplementary Figure 1) were exported to be statistically analyzed with Network Based Statistics (Zalesky et al., 2010).

Movement-Related Gamma Functional Connectivity
The Network Based Statistics (NBS) method was applied to detect eventual differences among the movement, planning and postmovement windows in their wPLI-based functional connectivity subnetworks [i.e., subsets of nodes and their connecting edge of a network that are separable from the remaining nodes and edges (Schirmer et al., 2019)]. This approach permits multiple hypothesis testing at the level of interconnected subnetworks, while controlling the family-wise error when performing analyses associated with a particular effect or contrast of interest (Zalesky et al., 2010).
NBS performs a mass univariate testing in order to identify the connections exceeding a test statistical threshold belonging to a given connected component. Then, a corrected p-value is computed for each component using the null distribution of maximal connected component size, which is empirically derived via a nonparametric permutation method. Here, a test with 5,000 random permutations was performed to compute statistical significance for the identified network component. Finally, the hypothesis test is performed for the empirically determined components by comparing their extent with the proportion of permutations yielding a component with equal or greater size, correcting for the family-wise error rate at cluster level with p < 0.05. Importantly, here we defined the subnetworks in terms of connected graph components, including nodes and edges, associated with statistical effects above a predefined threshold t-score of 4.8 (Zalesky et al., 2010).
Connectivity matrices of all subjects for the planning, execution and post-movement control windows were entered into one-way repeated measures ANOVA within the NBS approach. Post-hoc paired t-tests were performed to assess between-group differences.
The results of the NBS procedure are presented as threedimensional graph visualizations, which represented p < 0.05 connection pairs surviving multiple comparison correction. For each subject, we also computed the average of weights, in terms of wPLI values, of the edges belonging to each subnetwork identified by NBS. The resulting functional connectivity values were then used to fit linear regression models to unveil a possible relation between gamma phase-synchronization in the identified subnetwork and kinematic performance (see below).

Linear Regression Model on EEG Data
Finally, in order to unveil a possible relationship between the observed changes in gamma oscillatory activity and the peak velocity amplitude and movement time, we fitted linear regression models to predict the kinematic parameters based on the gamma power during the planning and execution timewindows and the network connectivity values. Specifically, for the regression models with gamma amplitude as the predictor, we entered the average gamma power for the significant time-window (planning and execution), gamma sub-band and electrodes as independent variable and peak velocity and movement time as dependent variables. Furthermore, to establish whether the average gamma power during the significant movement could be a good predictor, we also extracted the amplitude of gamma oscillations at the timing of the peak velocity for both the out and back movements. Specifically, for each trial, we computed the timing of the peak velocity for both the outgoing and return movements and extracted the corresponding broad and high-gamma amplitude in the electrodes that showed significant gamma ERS in the movement window (52-500 ms). These values were averaged across trials and used as linear predictors of the peak velocity amplitude. For these analyses, we had to exclude additional 3 participants due to technical problems in their raw kinematic data files (N = 59).
For the connectivity regression analyses, we extracted the average debiased wPLI value for the significant subnetwork differences between planning and execution and used to predict peak velocity and movement time.
As reported for the kinematic analysis, in order to satisfy normality assumption, the data were z-transformed before model fitting. Further, two subjects were removed from all regression analyses as their kinematic average values were more than 2 SD lower than the mean performance.

Movement Characteristics Depend on Target Distance
Sixty-four subjects performed 96 fast and uncorrected out and back movements from a central starting point to one of 24 targets (three distances and eight directions, Figure 1A) appearing in a random order. In general, the resulting movements were straight and their temporal velocity profiles were on average bell-shaped and appropriately scaled to the target distance, as displayed in Figure 1A and in previous publications (Tatti et al., 2019(Tatti et al., , 2020(Tatti et al., , 2021. Repeated measure ANOVAs showed that movement extent, peak velocity amplitude and timing, as well as movement duration increased with target distance  Figure 1B). Post-hoc tests revealed significant differences between the three target extents for all these measures. Target extent had a significant effect also on reaction time (Supplementary Table 1), with shorter reaction time for more distant targets. However, post-hoc analyses showed significant differences only between short targets and the other two targets.
We then performed linear mixed-effect regression modeling to characterize the contribution of peak velocity amplitude and movement time to movement extent.
Thus, we first assessed: adjusted R 2 , Bayesian Information Criterion (BIC, an index used in Bayesian statistics to select among two or more models), and the Theoretical Likelihood Ratio Test (TLRT, commonly used to compare the goodness of fit of two statistical models) (see section "Materials and Methods"). Following the indication provided by these indices, the model with random intercept and slope showed the best fitting of the data (Supplementary Table 2). The results showed that both peak velocity and movement time were strong predictors of movement extent variability: the model with random intercept and slope for each participant explained approximately 95% of the movement extent variance (R 2 adjusted = 0.95, Supplementary Altogether, these findings show that movement extent resulted mainly from the scaling of the force to the appropriate target distance with a lesser contribution of movement duration.

Movement-Related Gamma Oscillatory Dynamics
We next investigated the progression and the topographical distribution of gamma oscillatory activity (25.5-80 Hz) during movement planning and execution. We thus compared gamma oscillatory activity of each 24 ms-time bins to the average spectral power of the entire epoch (from 500 ms before the movement onset to 2,400 ms after) with non-parametric Monte Carlo permutation testing. This analysis revealed significant electrode clusters in two time-windows: one from 152 to 52 ms before movement onset, that is, during the reaction time period (Cluster t = 907.27, p = 0.0003, CI: 0.0004) and the other from 52 to 500 ms after moment onset, that is, during the movement execution time (Cluster t = 13.120, CI:19.598, p = 0.0001). As shown in Figure 2, before movement onset, gamma oscillatory activity increased in a cluster of electrodes over the centro-parietal region; during movement execution, the increase in gamma activity at first, involved electrodes over the occipital region, then spread to most of the scalp electrodes and ultimately returned to a centro-parietal cluster.
We then averaged the topography of gamma activity in the time bins of the reaction time period (from 152 to 52 ms before movement) and that in the time bins from 52 to 500 ms after moment onset, during the movement execution time. Cluster-based permutation analyses on four gamma bands (broad gamma: 25.5-80 Hz, low gamma: 25.5-40 Hz, medium gamma: 40.5-55 Hz, high gamma: 55.5-80 Hz) revealed for both time windows significant cluster of electrodes in the broad, medium and high gamma bands (Figure 3). In all these gamma bands, the planning window displayed greater synchronization over a centro-parietal region, whereas movement execution showed a widespread gamma synchronicity over all scalp channels, with maximum amplitude increase over the occipital region. No significant clusters were found in the low gamma range.
These results prompt some considerations about the role of gamma activity in both movement planning and execution. In line with a few reports (Gaetz et al., 2013;Gunduz et al., 2016), the occurrence of gamma synchronization over a centroparietal region during movement planning suggests that gamma oscillations may reflect the generation of a motor output. In addition to the creation of a motor plan, gamma ERS might also be involved in the control of the actual movement, as suggested by the presence of gamma ERS increase during the motor act, thus indicating the possible engagement of online feedback control processes.

Movement-Related Gamma Event-Related Synchronization Shifts Toward Higher Frequencies During Movement Execution
We then determined whether the observed gamma ERS during planning and execution differed in their peak frequency ranges (see section "Materials and Methods"). We used non-parametric related-samples Friedman's ANOVA to test possible differences of gamma ERS frequency between the planning, execution and post-movement time windows. As displayed in Figure 4, we found that gamma ERS peaked at different frequency ranges [χ 2(2) = 93.23, p < 0.0001]. During the planning window, gamma ERS was maximally expressed in the medium gamma range (53.73 Hz ± 10.04, mean ± SD) whereas during movement execution the peak shifted toward higher frequencies (65.8 Hz ± 10.72). Importantly, after movement completion, cortical activity synchronized back to the beta/low gamma frequency range (28.06 Hz ± 11.47; Dunn-Bonferroni posthoc tests, Planning vs. Movement: z = -0.594, p = 0.002; Post-movement vs. Planning: z = 1.68 p < 0.0001; Postmovement vs. Movement: z = 1.086 p < 0.0001). These findings demonstrate that gamma ERS occurring during movement planning and execution is maximally expressed in two distinctive frequency ranges, thus suggesting possible different functional properties of gamma ERS.

Movement Execution and Post-movement Windows Are Characterized by Greater Long-Range Connectivity Than Planning
To further explore gamma activity differences between planning, execution and post-movement windows, we compared their degree of gamma phase coupling across scalp channels. Phasesynchronization in the broad gamma and high gamma frequency   range was measured with the weighted phase lag index (wPLI), a functional connectivity metric that minimizes the impact of volume conduction effects (see section "Materials and Methods"). The wPLI was computed across all channelpairs to provide a whole-brain mapping of functional EEG networks and to ascertain possible differences of connectivity between the planning, execution and post-movement windows. In the broad gamma frequency, rmANOVA revealed significant functional connectivity differences between planning, execution and post-movement windows (p < 0.001, corrected for multiple comparison) in almost the whole brain. Post-hoc network-based statistics (NBS) analyses revealed a subnetwork of increased functional connectivity during movement execution compared to movement planning in the broad gamma frequency. The subnetwork consisted of 244 edges connecting 92 different electrodes (p < 0.001, corrected for multiple comparison). Interestingly, apart from a few pathways linking frontal to occipital regions, these patterns of increased connectivity mainly involved a left fronto-temporal-parietal network ( Figure 5A). No significant subnetworks of increased connectivity during movement planning were detected compared to the execution window. Post-hoc analysis also revealed a subnetwork of increased functional connectivity during the post-movement window compared to movement planning in the broad gamma frequency band. The network component consisted of 232 edges connecting 88 different electrodes (p < 0.001, corrected for multiple comparison). These patterns of increased connectivity mainly involved bilateral fronto-temporal and fronto-central electrodes ( Figure 6A). No significant connectivity differences were found between movement execution and postmovement windows.
Since movement-related gamma ERS shifted toward a higher frequency range during movement execution, we thus looked for possible subnetwork differences in the high-gamma range . In this frequency range, rm-ANOVA revealed significant widespread functional connectivity differences between planning, execution, and post-movement windows (p < 0.001, corrected for multiple comparison). Post-hoc NBS analysis highlighted a subnetwork of increased connectivity during movement execution compared to planning (p < 0.001, corrected for multiple comparison) that consisted of 571 edges connecting 122 different electrodes. In addition to patterns of connectivity mainly involving the fronto-temporal-parietal areas bilaterally, this subnetwork also involved fronto-occipital connections (Figure 5B), in line with our finding of greater high-gamma power increase over the occipital region during movement execution. No significant subnetworks of increased connectivity during the planning window were detected compared to movement execution. Post-hoc analysis highlighted a subnetwork of increased functional connectivity during the post-movement window compared to movement planning also in the high-gamma range. The network component consisted of 122 edges connecting 50 different electrodes (p < 0.001, corrected for multiple comparison), that mainly involved the right fronto-temporal network ( Figure 6B). No significant connectivity differences were found between the movement execution and post-movement windows.

Movement-Related Gamma Event-Related Synchronization Is Modulated by Target Distance
We then investigated possible relationships between gamma activity indices and movement features. We first ascertained whether target distance would affect the amplitude of gamma ERS by comparing gamma oscillatory activity for short, medium and long target distances with cluster-based permutation analyses based on the significant two temporal windows, two electrode clusters and gamma bands identified in the previous analyses (Figure 3).  Both the omnibus repeated measure ANOVA and the subsequent post-hoc tests confirmed that, during both movement planning and execution, broad gamma oscillatory activity increased with target distance (Table 1A) with significant differences between short and medium and short and long target distance trials (Table 1B and Figure 4). Despite greater gamma power for long than medium distances in both planning (mean difference ± SD: 8.1 ± 18.4%) and movement execution (7.9 ± 13.0%), post-hoc tests did not reveal any significant differences between medium and long-distance targets. Similar results were obtained for medium gamma and high gamma frequency bands, except for medium and short movement trials during planning in the medium gamma frequency band (Table 1B).
Nonetheless, these results show that gamma ERS amplitude is modulated by target distance and further indicate that such modulation might reflect some kinematic features planned before movement onset.

Movement-Related Gamma Event-Related Synchronization and Phase Connectivity Predicts Peak Velocity Amplitude and Timing
As both kinematic parameters and gamma ERS showed strong dependence on target extent, we explored the relationship between the amplitude and phase-synchronization of gamma oscillations and kinematic features by fitting linear regression models.
To model the relationship between gamma ERS amplitude and kinematics, we first extracted the z-transformed power of all the gamma bands from the two significant clusters of electrodes and used it as a predictor of peak velocity amplitude and movement time. We found that gamma amplitude during the planning and movement windows were significant predictors of movement time and peak velocity amplitude for all the tested gamma frequency bands ( Table 2 and Supplementary Figure 2). Regression analyses on the gamma amplitude at the timing of the peak velocity for both the out-going and return movements confirmed that gamma power is a good predictor of peak velocity, albeit the model with average gamma amplitude provided a better fitting (see Table 2 and Supplementary Table 3 for a comparison). Interestingly, the regression models on the timing of peak velocity revealed that only the gamma ERS during the out-going movement, but not during the return movement, was a significant predictor of the peak velocity (Supplementary Table 3). As the instructions provided to the participants were focused on the out-going movement (i.e., reach the displayed target as fast as accurately as possible), it might be that the out-going movement was more accurately planned and executed than the return one.
We also modeled the amplitude of gamma phasesynchronization in the -subnetworks that were activated during the movement and post-movement windows. In this case, we extracted and z-transformed the wPLI values from the network edges that showed greatest broad gamma and high gamma functional connectivity during movement execution and the post-movement window compared to planning. Linear regression models showed that high gamma wPLI values in the subnetwork activated during the movement significantly predicted peak velocity amplitude (Table 3).
Interestingly, no significant correlations were observed for the whole gamma band (Table 3), thus suggesting that peak velocity might specifically depend on the amount of coherent high gamma activity in the highlighted occipital and contralateral parieto-frontal network. Both for gamma amplitude and wPLI, greater values were associated with higher peak velocity amplitude, in line with the evidence supporting the prokinetic role of gamma oscillations. Interestingly, no significant linear relationship was observed for movement time, thus indicating that gamma synchronicity might be linked to the specification of the scaling force necessary to reach the target at its location. Moreover, the right fronto-temproal subnetwork activated after movement termination was also not related to the kinematic characteristics of the movement (Table 3).

DISCUSSION
This study used a parametric approach to demonstrate in a large sample of healthy subjects that the amplitude of movementrelated gamma ERS scales with target distance. While confirming that gamma synchronization occurs both during movement planning and execution, the present findings provide the first evidence that different gamma peak frequency and topography characterize the planning and execution phases. This study also reveals the presence of a gamma-coordinated long-range network involving occipital, frontal and central regions during movement. Both during the planning and execution phases, the amplitude of the movement-related gamma ERS and its connectivity pattern parametrically increased with target distance. This is the first evidence in humans supporting the notion that both gamma synchronization amplitude and gamma phase coherence pattern are significant predictors of peak velocity amplitude and duration of goal-directed reaching movements.
It has been proposed that neocortical gamma oscillatory activity is central to sensorimotor and cognitive functions, specifically for information processing through the integration of neuronal assemblies' activity, enabling a transfer of information across connected regions. Movement-related changes in the amplitude of gamma oscillations have been observed with ECoG in a variety of motor tasks specifically linked to some kinematic features such as movement trajectory (Talakoub et al., 2017), velocity (Bundy et al., 2016;Wang et al., 2017a,b;Lofredi et al., 2018), and force (Flint et al., 2014;Jiang et al., 2020).
A notable finding of the present study is that cortical gamma synchronization occurs during both planning and execution of the reaching movements and is parametrically modulated by movement extent. Increased oscillatory gamma activity, mostly in the medium (40.5-55 Hz) and high (55.5-80 Hz) frequency ranges, was visible immediately before and after movement onset, albeit with different topographies that may reflect different functions carried out during movement planning and execution. Interestingly, during movement planning, the first burst of gamma ERS was maximally expressed over the centro-parietal region, whereas the two gamma ERS bursts corresponding to the out-and back-movements were mostly visible over the parietooccipital region.
The occurrence of gamma synchronization during the planning phase has been reported in previous studies (Brücke et al., 2008;Gaetz et al., 2013;Gunduz et al., 2016;Ryun et al., 2017) and supports the idea that gamma ERS might be associated with top-down mechanisms engaged for movement preparation and, thus, with the creation of an efferent copy of   the upcoming movement. Indeed, our task instructions (that is, to perform fast, uncorrected, out-and back-movements to targets appearing in an unpredictable order) mostly engaged feed-forward specification of a force-scaling factor according to movement extent (Ghez and Gordon, 1987;Gordon et al., 1994a,b). Planning of movement distance between the starting hand position and the target requires proper jerk, acceleration and velocity scaling prior to movement onset (Gielen et al., 1985), a process that engages several regions including the primary sensorimotor cortices, supplementary motor area and the basal ganglia (Bell et al., 1994;Turner et al., 1998;Ghilardi et al., 2000;Tankus et al., 2009). Accordingly, the topographical distribution of gamma ERS before movement onset and its correlation with movement features suggest that the first burst of gamma ERS might reflect the engagement of sensorimotor neuronal assemblies appropriate to the characteristics of the upcoming movement. Indeed, our regression analyses confirmed that the amplitude of the pre-movement ERS positively predicts the amplitude of the peak velocity and, conversely, negatively predicts the movement duration.
Although the observed scaling of the pre-movement gamma ERS supports our interpretation of its role in movement planning, we cannot completely exclude that gamma ERS prior to movement onset might arise from subtle muscle contractions occurring before the motor act that can be captured only by comprehensive electromyography (EMG) recordings (Ryun et al., 2017).
Therefore, future analyses should combine EEG, EMG and kinematic recordings in order to clarify the relationship between pre-movement cortical gamma oscillations and muscle activation.
After the initial gamma ERS, gamma oscillatory activity visibly increased after movement onset and strikingly tracked the temporal profile of movements' velocity. Compared to the initial burst that peaked around 54 Hz, this movementrelated synchronous activity was most prominent in the high-gamma range, in agreement with different studies (Pfurtscheller et al., 2003;Muthukumaraswamy, 2010;Joundi et al., 2012;Cheyne and Ferrari, 2013). Although gamma ERS was visible on the majority of the scalp sensors, our finding of greater synchronization over the posterior regions might reflect the engagement of the associative posterior parietal area, which is thought to integrate efferent motor signals with visual and proprioceptive information (Hyävrien and Poranen, 1974;Mountcastle et al., 1975;Andersen et al., 1997;Wolpert et al., 1998;Galletti et al., 2003;Bernier and Grafton, 2010;Shadmehr et al., 2010;Elliott et al., 2017), as well as to encode and integrate visual and proprioceptive information during movement execution (Buneo and Andersen, 2012). This hypothesis is in line with the results of our phase synchronization-based connectivity analyses: compared to the planning window, during movement execution, we observed greater wPLI connectivity in the high gamma band between the electrodes located over the right parieto-occipital region and those over the bilateral fronto-central region. The same connectivity analysis on the broad gamma band displayed a more lateralized connectivity pattern between the left centro-parietal and fronto-central regions with a minor contribution of the parieto-occipital areas. Importantly, our regression analyses further showed that the network phase-synchronization index was a significant predictor of peak velocity in high gamma, but not in the broad gamma range.
Altogether, we may speculate that movement-related ERS in the high gamma range during movement execution could reflect the specific activity of the dorso-medial stream, integrating online afferent sensorimotor and visual input and sending realtime feedback to premotor and motor regions to control the ongoing movement. Moreover, the observed positive linear relation between the functional connectivity and ERS amplitude during the movement and peak velocity might reflect either the amount of neural assemblies recruited to produce more force to increase speed (Muthukumaraswamy, 2010) or the amount of coordinative control needed to sustain the incoming afferent information about the ongoing movement.
Interestingly, we also found increased connectivity (mostly in the high gamma range) in a subnetwork that mainly involved the right fronto-temporal areas after movement termination. Differently from gamma connectivity during the movement, this subnetwork was not linked to movement kinematics, thus it might reflect a different function. A possible interpretation of this result is that it might reflect working-memory processes occurring after movement termination and the update of the activated motor representations. Alternatively, this finding may represent re-orienting attentional processes to unattended locations that precede the appearance of the next target. Indeed, re-orienting attentional processes engage the ventral attentional network that is lateralized to the right hemisphere and includes the temporo-parietal junction as well as the frontal cortex at different levels (Corbetta et al., 2008). Although plausible, these conclusions remain for the moment only hypotheses to be tested in future experiments.

CONCLUSION
In conclusion, this study provides an extensive characterization of scalp recorded movement-related gamma synchronization during planar reaching movements in a large sample of healthy participants. Gamma ERS amplitude during movement planning and execution, as well as increased parieto-occipital and fronto-central connectivity pattern reliably predict the specific movement features.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by CUNY University Integrated Institutional Review Board (UI-IRB). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
MFG, ET, and AQ: study conception and design. ET, MFG, FF, AC, and CC: data collection and analyses. ET, MFG, AC, AQ, and FF: interpretation of the results and manuscript drafting. ET, FF, MFG, and AQ: manuscript revising. All authors agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.