SSVEP phase synchronies and propagation during repetitive visual stimulation at high frequencies

Steady-state visual evoked potentials (SSVEPs), the brain response to visual flicker stimulation, have proven beneficial in both research and clinical applications. Despite the practical advantages of stimulation at high frequencies in terms of visual comfort and safety, high frequency-SSVEPs have not received enough attention and little is known about the mechanisms behind their generation and propagation in time and space. In this study, we investigated the origin and propagation of SSVEPs in the gamma frequency band (40–60 Hz) by studying the dynamic properties of EEG in 32 subjects. Using low-resolution brain electromagnetic tomography (sLORETA) we identified the cortical sources involved in SSVEP generation in that frequency range to be in the primary visual cortex, Brodmann areas 17, 18 and 19 with minor contribution from sources in central and frontal sites. We investigated the SSVEP propagation as measured on the scalp in the framework of the existing theories regarding the neurophysiological mechanism through which the SSVEP spreads through the cortex. We found a progressive phase shift from posterior parieto-occipital sites over the cortex with a phase velocity of approx. 8–14 m/s and wavelength of about 21 and 24 cm. The SSVEP spatial properties appear sensitive to input frequency with higher stimulation frequencies showing a faster propagation speed.

Repetitive visual stimulation (RVS), or visual flicker, at a frequency in the range from 1 to 100 Hz elicits an oscillatory response in the electroencephalogram (EEG) at the same frequency and/or its harmonics known as steady-state visual evoked potentials (SSVEPs). This response is time and phase-locked to the driving stimulus. It appears most prominently over parietal and occipital areas as they are closer to the visual information processing centers of the brain cortex but can also be detected in central and frontal areas 1 .
Because of their high signal to noise ratio and robustness to artifacts, SSVEPs have proven beneficial in both research and clinical applications.
The wide spatial extent of the SSVEPs across the scalp can provide means for studying the collective behaviour of large groups of neurons 2 allowing the investigation of neural processes underlying brain rhythmic activity and functional brain connectivity. When RVS is superimposed on a cognitive task, especially around the alpha frequency band (8)(9)(10)(11)(12), specific topographic changes in amplitude and phase of the SSVEP appear analogous to the site-specific reductions in alpha EEG amplitude associated with cognitive and motor tasks 3 . This technique is known as steady-state topography and has been applied for studying working memory 3,4 , binocular rivalry 5 , selective attention 6,7 , vigilance 8 and even emotions 9 .
SSVEPs are also used as a diagnostic tool to study pathological brain dynamics and as an investigative probe of the cortical mechanisms of visual perception. The amplitude and phase of the SSVEP vary as a function of stimulus parameters such as temporal and spatial frequency, contrast, luminance, hue, visual angle and spatial attention 1,10-13 and can be a sensitive electrophysiological index of a variety of visual-perceptual and cognitive functions 14 . SSVEP-based techniques such as photic driving and sweep-VEP are used for evaluating visual acuity of infants and children 15 , studying infant vision 16 , assessment of neuro-visual function 17 , spatial neglect 18 and epilepsy 19 . SSVEPs are also applied for the investigation and treatment monitoring of age-related www.nature.com/scientificreports/ the spatiotemporal oscillatory properties of SSVEPs. Particularly, SSVEP responses at high frequencies, e.g. in the gamma band , which have not received enough attention and little is known about their generation and factors affecting the response. We set to investigate their site of origin and their cortical propagation properties in terms of speed and phase.

Results
SSVEP response. Time frequency decomposition analysis revealed a stable SSVEP response starting within 250 ms from stimulus onset at each fundamental and its harmonic frequencies. Figure 1 shows the percentage power change (event related synchronization, ERS) for the first second of stimulation for three of the conditions at each fundamental frequency. As expected the response is stronger over parieto-occipital sites. Lower frequency conditions are characterized by a broader response as compared to higher frequencies, showing some contribution from fronto-central sites. In all cases, after the transient response the first 250 ms, we observe a sustained steady state response for the whole duration of stimulation. These observations were consistent across all stimulation conditions. For a more comprehensive overview of the neural dynamics at the fundamental and the harmonic frequencies of this dataset as well as the interaction with ongoing oscillatory activity see Tsoneva et al. 52 .

Source analysis.
To identify the spread of cortical sources involved in SSVEP generation we performed source localization analysis using the sLORETA software. This analysis revealed strongest activation in response to frequencies between 40 and 56 Hz at coordinates: X = 0, Y = −85, Z = 30 , identified as the Cuneus in Brodmann area 19 of the Occipital cortex (see Fig. 2a,b). Other areas of the visual cortex such as V1 (Brodmann area 17) and V2 (Brodmann area 18) were also identified as contributors. Additionally, minor contribution from sources in central and frontal sites could be observed. These coordinates were consistent across all conditions under 56 Hz for both, the fundamental and subharmonic components (see Supplementary information Figure S1ab). We did not find a consistent activation for higher order harmonic frequencies.
The source distribution for stimulation at 58 Hz and 60 Hz was slightly different (see Fig. 2c). While activity in the visual cortex was present, the strongest activation was found in the inferior frontal gyrus (Brodmann area 47) for the fundamental frequencies and middle frontal gyrus (Brodmann area 11) for the sub-harmonic frequencies (see Supplementary information Figure S1c). Statistical testing, however, did not show a significantly higher current source density (CSD) for those coordinates when compared to the second before stimulation onset as baseline.
The average log-F-ratio statistics was quite consistent across all stimulation conditions highlighting cortical areas in the occipital, parietal, limbic and frontal lobe. The average log-F-ratio across all conditions is shown in Fig. 2d. The log-F-ratio thresholds for p < 0.01 and p < 0.05 for the fundamental component at each stimulation condition are listed in Table 1. The results for the sub-harmonic components of each stimulation condition can be found in the Supplementary information Figure S1d and Table S1. Higher order harmonics did not show any significant difference from baseline (the second before stimulation onset).

Source of origin.
While LORETA is very useful for localizing the activity, it cannot tell us about the chronology of events. Generally, the SSVEP response is time-and phase-locked to the driving stimulus. Figure 3a shows the instantaneous phase locking value (PLV) between midline electrodes and the stimulation for one of the stimulation conditions (44 Hz) averaged over all subjects. Phase locking ( PLV > 0.5 ) occurs already within the first 150 ms and saturates after approximately 250 ms, following the transient response. The same behaviour was observed for all other conditions. Because we are particularly interested in the phase properties www.nature.com/scientificreports/ of the response, we use finite impulse response (FIR) filters at the stimulation frequency, which have zero phase distortion after correcting for their linear delay. The inherently high order and the equiripple nature of those filters, however, can cause discontinuities at the head and tail of their impulse response. This is mostly visible at the start and the end of stimulation where the steady-state response is still weak. To deal with this issue we only consider PLVs for samples where the phase locking statistic (PLS) show significance of the phase covariance ( PLS < 0.05 ). Figure 3b shows the decrease of PLS shortly after the stimulus onset and its relatively low value for the whole stimulation period, especially for parieto-occipital channels.
To identify the initial source locking to the driving stimulus first, we look at the significant PLVs > 0.5 for all channels. The locking delays for different channels at stimulation frequency of 44 Hz are visualized in Fig. 3c. Note that for some channels PLS never reaches significance level. As expected, the channels with minimum locking delays are in the parietal and occipital areas for all frequencies, with channel Pz being the first to lock for 7 out of the 10 frequencies. For consistency, henceforth, we select channel Pz as source of wave propagation.
Phase velocity. The phase difference between each electrode and the source appears to be constant for the whole stimulation period (see Fig. 4b). The value, calculated over 5 cycles for each frequency and averaged over all trials and all subjects every 250 ms shows a progressive phase shift in the EEG, suggesting spherical propagation in all directions from the source.
To characterise the phase velocity of the propagation we estimate the spatial phase gradients by fitting a first order polynomial to the phase estimates for a set of electrodes covering the sites identified in the source analysis, namely O1, O2, Oz, PO3, PO4, P3, Pz, P4, CP1, CP2, Cz, FC1, FC2 and Fz as shown in Fig. 4a ordered in increasing distance from the source Pz.
The phase velocity is inversely proportional to the phase gradient and is calculated using Eq. (6) adjusting for the stimulation frequency for each condition. Figure 4c shows the first order polynomial fit and the phase velocity estimates over the area identified by the source analysis for all conditions 250 ms after stimulus onset when the steady-state response should be already stable. The correlation between phase difference and electrode distance (see Eq. 5) is significant for all conditions with mean correlation coefficient r = 0.85(±0.013) and    www.nature.com/scientificreports/ The phase velocity estimates are by definition sensitive to stimulation frequency (see Eq. 6). Figure 5 shows a linear relationship and a significant positive correlation between our velocity estimates and frequency of stimulation ( r = 0.99 , p = 3.852e − 08 ). We also estimate the linear regression model parameters that describe this relationship. An increase in the stimulation frequency of 1 Hz results in about 0.24 m/s increase in phase velocity. It is important to note that the validity of this model applies to the range of frequency tested as different mechanisms might apply to oscillations in different frequency sub-systems.

Discussion
We set to investigate the SSVEP responses to periodic visual stimulation at high frequencies within the gamma band (40-60 Hz) by studying the dynamic properties of EEG measured by sparse electrode montage in 32 subjects. Not surprisingly the strongest cortical sources involved in SSVEP generation, identified by the sLORETA analysis, were found in the primary visual cortex, Brodmann areas 17, 18 and 19. Minor contribution from sources in central and frontal sites could also be observed. Looking at the phase response we found phase locking occurring already within the first 150 ms with channel Pz being the one that synchronized to the driving stimuli first for the majority of the tested conditions. While it is difficult to distinguish between the neurophysiological mechanism through which the SSVEP spreads through the cortex, being a finite number of electrical dipoles that are activated sequentially in time, or through cortical travelling waves, the SSVEP response as measured by the surface EEG shows a progressive phase shift in all directions from the source. Our results suggest that highfrequency-evoked SSVEPs propagate over the cortex with a phase velocity of approx. 8-14 m/s and wavelength of about 21 and 24 cm.
It should be emphasized that estimates of phase velocity depend on a variety of different factors. Some authors highlight the fact that the choice of EEG recording reference might play an important role with bipolar recordings giving lower estimates as compared to reference methods 2,49 . Other estimates are based on reference-free methods such as surface Laplacian methods 38,49 . MEG and ECoG recordings produce lower velocity estimates than EEG measurements 40,46 . Apart from the above mentioned factors affecting the phase estimates, inter-electrode distance could also vary. Some authors choose to correct for cortical folding (fissures) 2,53 , while others do not 40,43 . Furthermore, using euclidian or great-circle distance between electrodes results in different velocity estimates. All these factors inherently affect the phase velocity estimates making it difficult to compare between studies.
There is also lack of literature reporting SSVEP propagation in the gamma band. However, frequency seems to be an important factor and a number of authors show that SSVEP properties are sensitive to flicker frequency 38,46 . Sleep slow waves (0.5-4 Hz) are reported to be travelling at 1-7 m/s, with most of the scalp velocities being in the 2-3 m/s range 43 . Alpha waves seem to be travelling at 3-8 m/s, with lower alpha (8 Hz) around 6.4 m/s and upper alpha (12 Hz) around 8.4 m/s 49 . One of the few studies looking at higher frequencies reports phase velocities ranging from 0.2 m/s at 0.5 Hz to 14 m/s at 32Hz for eye-movement evoked responses. This progressive increase of velocity with the frequency of oscillations (spontaneous or evoked) seems to be in line with our estimates and the correlation with frequency of stimulation we found (see Fig. 5). It is, however, important to note that SSVEPs www.nature.com/scientificreports/ are sensitive to cortical resonances in at least three frequency subsystems 10,38 and the mechanisms that apply at higher frequencies might not directly translate to oscillations at lower frequencies.
In their theoretical framework, considering the rough boundaries of propagation velocities within corticocortical fibers and accounting for experimental and theoretical uncertainties, some authors 53 expect the phase velocity of cortical waves traveling away from primary visual areas to not exceed 30 m/s (up to a factor of three from observed estimates between 6 and 9 m/s). Our estimates of SSVEP propagation speed using Laplacian derivation and great-circle distance are within these loose boundaries.
Scalp recorded patterns of cortical activity also reflect volume conduction across the scalp, skull and other tissues. Volume conduction is especially a problem of reference EEG recordings, which are known to overestimate propagation velocities. Alternatively, close bipolar electrodes generally do better, but may be biased toward shorter wavelengths 54 . In this study, the spatial distortion from conductivity variations was minimized by applying a surface Laplacian, which is reference independent linear data transformation that maintains the invariant aspects of the EEG signal and provides means for obtaining continuous estimates of radial current flow at the scalp for both low-and high-density EEG montages 55 . It enhances the spatial information of the EEG signal by emphasizing sources within 2-3 cm of the electrode and providing good representation of the underlying current dipoles 38 . In our initial analysis we compared two preprocessing approaches, namely the widely used common average reference (CAR) and surface Laplacian. Signals preprocessed using CAR seem to show progressive phase shifts even outside the stimulation period (350 ms before stimulation onset) showing clearly the effects of volume conduction (see Fig. 6). Applying surface Laplacian eliminated this effect and the phase shift was only observed during the stimulation period (see Fig. 4c). Additionally, the source localization analysis, revealing neuronal current generators underlying EEG scalp topographies, confirmed the spread of the observed propagation and was used as a basis for our further analysis.
The sources identified by sLORETA for the fundamental and sub-harmonic components seemed to agree with what has been previously reported 14,56 . We, however, did not find a significant activation for the first harmonic component. A reason for that might be that all our stimulation condition were in the high frequency region, resulting in first harmonic component in the range from 80 to 120 Hz. The decrease in SNR for higher order harmonics may distort the LORETA analysis. Wagner et al. suggest that simultaneously active sources can only be separated well by sLORETA analysis if their fields are distinct enough and of similar strength, and in the context of a strong or superficial source, weak or deep sources remain invisible 57 . In our previous work we did observe a linear decrease in SSVEP strength at the first harmonic especially at higher frequencies 52 . It was encouraging to see that the cortical sources identified by the sLORETA analysis for the first harmonic of our lowest stimulation condition (40 Hz) were consistent with those identified by Pastor et al. 56 for the same frequency (their highest frequency condition).
Regarding the differences observed for the highest stimulation conditions at 58 and 60 Hz, stimuli at those frequencies at average light level of about 130 cd/m 2 and assuming an average pupil diameter of about 4 mm, would fall very close to the visual perception threshold as defined by the temporal sensitivity curve at retinal illuminance of 1000 Td and above 58 . The responses in the high frequency range could possibly have different neural representations and the associated neural mechanisms could be different from those governing conscious visual perception at lower frequencies 11,59 .
Broadband pseudo random noise sequences produce a code-modulated visual evoked potential response (cVEPs) that shares many characteristics with SSVEPs, including topology. This approach has been successful in BCIs 60   www.nature.com/scientificreports/ While no conceptual framework could possibly account for the complex properties and interactions in the human brain, models considering both the spatial and temporal aspects of the cortical dynamics would provide better representation of the underlining mechanisms.
The broad range of SSVEP applications in research and clinical practice call for reduced complexity, which would positively impact the ease of use, patient throughput, comfort and safety. In this manuscript, we investigated the spatiotemporal oscillatory properties of SSVEPs using a sparse electrode montage and high stimulation frequencies. We show what are the cortical sources underlying HF-SSVEPs as well as their spread and propagation velocity. This work could help building more practical solutions, improve processing speed, reduce erroneous interpretation and enhance the patient experience. Furthermore, this can open new doors for studying the neural dynamics in the gamma band, which functional significance is only starting to be understood.

Methods
Participants. For this study we invited 32 healthy volunteers with normal or corrected to normal vision (14 male, 18 female, age range 20-30 years old). Because of the nature of the stimuli, the participants were carefully screened for a history of photosensitivity, seizures, migraine, or headaches and only healthy participants not reporting any of these conditions were enrolled in the study. They were informed about the goals of the study and asked to sign an informed consent before the start of the experiment. The study was approved by the Philips internal ethics committee on biomedical experiments and conducted in accordance with the declaration of Helsinki Guidelines for Ethics in Research.
Experimental paradigm. During the experimental session the participants were instructed to focus their attention on an LED panel with dimensions 44 × 36 cm, and to avoid eye blinks and movement as much as possible. The panel was placed centrally on a desk at a distance of about 70 cm and spanning a visual angle of about 35 • . The panel consisted of twenty 2 mm white LEDs (Philips Lumileds LXHL) evenly distributed in 4 rows and 5 columns (approx. 8 cm apart in each dimension) covered by a diffusing screen. It delivered repetitive visual stimulation in the form of square waves at 100% modulation depth and 50% duty cycle. The average light output of the stimuli was approx. 130 cd/m 2 . The experimental conditions comprised RVS at 10 different frequencies from 40 to 60 Hz with a step size of 2 Hz, excluding 50 Hz as it coincides with the power line frequency at the study site. Stimuli presentation was randomized and each condition was presented 10 times, which is expected to improve the SNR by a factor of √ (10) ≈ 3 61 . Each trial consisted of an inter-stimulus interval, where no light was present, of random duration between 3 and 6 s, followed by a 2 s long stimulation period (see Fig. 7). Trials were grouped in blocks of 10 separated by a break of 10 s during which the participants could rest, move or blink. The average experimental duration was approx. 13 min and ensured that effects of fatigue or decreased alertness, which might affect the response 62 , are minimized. Data acquisition. BioSemi ActiveTwo signal acquisition system (BioSemi B.V., Amsterdam, The Netherlands) was used to record EEG via 32 sintered Ag/AgCl active electrodes at a sampling frequency of 2048 Hz. The electrodes were mounted on an elastic cap according to the 32-channel extension of the international 10-20 positioning system (see Fig. 4a). The setup used an added common mode sense and driven right leg as "ground" electrodes. The signal from a photo diode was jointly recorded with the EEG to mark the presentation of the stimuli. Data analysis. EEG preprocessing. The recorded signals were analysed using custom MATLAB scripts and the open source EEGLAB toolbox 63 . First, a 50 Hz notch filter was applied to attenuate the power line interference. Eye-blinks and other ocular artifacts were removed using independent component analysis (ICA) 53 by excluding the component showing the highest average spatial difference between front and back electrodes. The signals were subsampled to 256 Hz to reduce the computation time, and surface Laplacian 64,65 was applied to reduce the effect of widely distributed background neuronal activity and to enhance the spatial localization of the signals.
The signals were then segmented into non-overlapping epochs time-locked to stimulus onset. The epochs started 1000 ms before stimulus onset and ended 3000 ms after, including 2000 ms of stimulation and a 1000 ms long period after stimulus offset.
Time frequency decomposition. All preprocessed epochs were subjected to a time-frequency decomposition using the fast Fourier transform (FFT) with a 126 samples (500 ms) wide Hanning window and 4.5 samples (17.59 ms) overlap (a total of 200 windows). The average power at each stimulation frequency was, then, calculated for five 250 ms-wide windows, starting with a baseline interval 250 ms before and continuing for 1000 ms after stimulation onset. Using this estimate we define the event related synchronization (ERS) as the fraction power change in each 250 ms window referenced to the 250 ms baseline interval. The resulting dynamics in www.nature.com/scientificreports/ each channel were then visualized on a topographic map. For this analysis the window length and baseline were optimized to improve the visual representation of the spatial maps.
Source localization. We estimated EEG cortical source current density with low-resolution brain electromagnetic tomography using the sLORETA software 66 . sLORETA calculates the standardized current source density (CSD) in the gray matter and the hippocampus at each of the 6239 voxels of the Montreal Neurological Institute (MNI) reference brain. CSD at the frequency of stimulation and its harmonics for each condition was estimated using a regularization factor of 10 (SNR = 10). For this analysis we computed the average cross spectra for each condition over 1 second long epochs starting 250 ms after stimulus onset. The results were then compared to the 1 second period before stimulation onset using sLORETA's built-in non-parametric statistical analyses (statistical non-parametric maps, SnPM) using log-F-ratio statistic for paired groups with a threshold of p < 0.01 and p < 0.05 and permutation test repeated 5000 times, correcting for multiple comparisons.
Phase synchrony analysis. Phase synchrony between two signals was estimated by looking for latencies at which the phase difference between signals varies little across trials 67 . For this analysis, all signals including the signal from the photo diode monitoring the light stimulation, were first filtered with a bandpass Parks-McClellan finite impulse response (FIR) filters centered around each stimulation frequency ( Fstop : F ± 3Hz, Fpass : F ± 1Hz, Dpass : 0.058, Dstop : 10 −4 , densityfactor = 20 ). The filter properties were chosen such that the filter is stable, while ensuring a narrow bandpass response. FIR filters were chosen because they can be designed to have an exact linear phase response. The constant group delay after filtering was estimated and signals were corrected for it. Next, the instantaneous phase of the signal y(t) was computed using the Hilbert transform 68 . The Hilbert transform y h (t) is the convolution of the signal y(t) with the function h(t) = 1/πt such that: where τ is the integration variable.
The Hilbert transform enables the analytical representation of the signal as follows: Comparison of the Hilbert transform and the FFT-based method revealed no substantial difference in phase estimates. The phase difference �φ 12 (t) between two time series y 1 (t) and y 2 (t) with instantaneous phases φ 1 (t) and φ 2 (t) , respectively, was calculated as follows: The phase locking value (PLV) for each condition, channel and subject at time t was then defined as the average value over trials ( n = 1, . . . , N): The phase locking statistics (PLS) 67 measures the significance of the phase covariance between the two signals. It is a statistical test based on a surrogate data aiming at differentiating significant PLVs against noise. PLS was calculated by shuffling the trials for another randomly selected stimulation condition and calculating the respective surrogate PLV. The PLS comprises the proportion of the surrogate values higher than the original PLV. We chose PLS ≤ 5% to ensure the probability of having false positive rate below 5%.
Distance estimates. The distance �δ 12 between electrode locations was estimated as the great-circle distance between two points on a sphere using the haversine formula (see Eq. (5)). We assumed a spherical head model with a radius of 9 cm.
where r is the radius of the sphere; φ 1 , φ 2 are the latitudes of the two points in radians; and 1 , 2 are the longitudes of the two points in radians.
Phase velocity. Phase velocity v is the speed at which the phase of a wave propagates over the scalp and is calculated as follows: where f is the stimulus frequency, and �φ 12 (t) is the phase difference as defined in Eq. (3) and �δ 12 is distance between electrodes as defined in Eq. (5).

Data availability
The data that support the findings of this manuscript are not publicly available, as it will potentially compromise research participant consent under the European GDRP law regarding reuse of data for secondary research and development purposes.