Response to photic stimulation as a measure of cortical excitability in epilepsy patients

Studying states and state transitions in the brain is challenging due to nonlinear, complex dynamics. In this research, we analyze the brain's response to non-invasive perturbations. Perturbation techniques offer a powerful method for studying complex dynamics, though their translation to human brain data is under-explored. This method involves applying small inputs, in this case via photic stimulation, to a system and measuring its response. Sensitivity to perturbations can forewarn a state transition. Therefore, biomarkers of the brain's perturbation response or “cortical excitability” could be used to indicate seizure transitions. However, perturbing the brain often involves invasive intracranial surgeries or expensive equipment such as transcranial magnetic stimulation (TMS) which is only accessible to a minority of patient groups, or animal model studies. Photic stimulation is a widely used diagnostic technique in epilepsy that can be used as a non-invasive perturbation paradigm to probe brain dynamics during routine electroencephalography (EEG) studies in humans. This involves changing the frequency of strobing light, sometimes triggering a photo-paroxysmal response (PPR), which is an electrographic event that can be studied as a state transition to a seizure state. We investigate alterations in the response to these perturbations in patients with genetic generalized epilepsy (GGE), with (n = 10) and without (n = 10) PPR, and patients with psychogenic non-epileptic seizures (PNES; n = 10), compared to resting controls (n = 10). Metrics of EEG time-series data were evaluated as biomarkers of the perturbation response including variance, autocorrelation, and phase-based synchrony measures. We observed considerable differences in all group biomarker distributions during stimulation compared to controls. In particular, variance and autocorrelation demonstrated greater changes in epochs close to PPR transitions compared to earlier stimulation epochs. Comparison of PPR and spontaneous seizure morphology found them indistinguishable, suggesting PPR is a valid proxy for seizure dynamics. Also, as expected, posterior channels demonstrated the greatest change in synchrony measures, possibly reflecting underlying PPR pathophysiologic mechanisms. We clearly demonstrate observable changes at a group level in cortical excitability in epilepsy patients as a response to perturbation in EEG data. Our work re-frames photic stimulation as a non-invasive perturbation paradigm capable of inducing measurable changes to brain dynamics.


Introduction
Understanding brain function is challenging due to its' nonlinear, complex dynamics.Methods from statistical physics and dynamical systems theory have been successful in characterizing other complex, nonlinear systems, including ecological population changes (Sugihara et al., 2012), weather (Peters and Neelin, 2006), and economics (Lux and Marchesi, 1999).These theories have also been successfully applied to the study of the brain, usually in animal models or using invasive procedures, including the study of brain state transitions in epilepsy patients (Kalitzin et al., 2002;Da Silva et al., 2003;da Silva and Harding, 2011;Maturana et al., 2020).
Epilepsy is a brain disease characterized by abnormal electrical brain activity and unprovoked recurrent seizures, that occurs in approximately 1% of the population (Fiest et al., 2017).These can be detected in electroencephalography (EEG), and the transition to seizure can be conceptualized as a critical transition (Da Silva et al., 2003;Kramer et al., 2012;Jirsa et al., 2014).Metrics that identify changes in system dynamics before these transitions, called "early warning signs, " have been observed in the EEG of epilepsy patients before seizures (Meisel et al., 2015;Maturana et al., 2020).These include increased variance, autocorrelation, and sensitivity to perturbation (a small external input) (Scheffer et al., 2009).
As a system approaches a state transition, it becomes more sensitive to perturbations, such that it is slower to return to equilibrium (Scheffer et al., 2009).Since increased sensitivity to perturbation can forewarn a state transition, biomarkers of the brain's perturbation response or "cortical excitability" (Freestone et al., 2011) could be used to track brain dynamics toward transitions to different states, including seizure states or other altered states of consciousness.Studying epilepsy in humans often involves trying to capture seizures during brain recordings, which can never be guaranteed and often requires long-term recordings.Perturbation methods have the benefit of reliably inducing responses of interest in real time.However, delivering perturbations to the brain is generally invasive in practice, requiring implantation of intracranial electrodes or other invasive devices.Animal studies in mice have successfully used electrical stimulation as a perturbation to track cortical excitability (Lamers et al., 2022).In humans, perturbation studies have used a subset of patients with focal epilepsy eligible for resection surgeries who are undergoing monitoring to confirm the seizure onset zone, where electrodes are placed strategically for that patient's unique seizure onset zone (Kalitzin et al., 2005;Cook et al., 2013;Bergey et al., 2015;Wendling et al., 2016;Oderiz et al., 2019).Noninvasive perturbation is also possible through technologies such as transcranial magnetic stimulation (TMS) (Manganotti et al., 2013;Ozdemir et al., 2021;Perellón-Alfonso et al., 2021), although this technology is relatively expensive and the perturbation itself is arguably diffuse and spatially non-specific.
Photic stimulation is a widely used diagnostic technique (Fisher et al., 2022) that could be used as a non-invasive perturbation paradigm to probe brain dynamics during routine EEG studies in humans.It involves changing the frequency of strobing light, sometimes triggering a photo-paroxysmal response (PPR), which is an electroencephalographic event visible on EEG.Generally, its primary clinical relevance is to diagnose photosensitive epilepsy, although we argue that it offers a safe, inexpensive research method for brain perturbation studies that is accessible to a much broader range of people.Photosensitivity is prevalent in the general population but is more common in those with epilepsy, with photosensitivity being significantly associated with genetic generalized epilepsy (GGE) (De Kovel et al., 2010).Photosensitivity has been shown to have a strong genetic basis, with family studies suggesting evidence of common gene loci in some subsets of families, though these vary geographically and depending on clinical phenotypes (De Kovel et al., 2010;Fisher et al., 2022).Sex has also been shown to be a risk factor for photosensitivity (Cerulli Irelli et al., 2023), with females exhibiting almost double the prevalence compared to males, and most prominently presenting during adolescence (Fisher et al., 2022).This highlights the importance of investigating sub-groups of epilepsy patients where possible, in investigations of photic stimulation.
Unlike stimulating the brain with electrodes or with TMS, photic stimulation itself is input via the visual system of the brain, engaging with the brain's functionality in a more "natural, " endogenous way.In previous studies, markers of critical transitions have previously only been observed in patients with focal epilepsy; we aim to observe them during the generalized PPR response.Also, unlike previous studies, we aim to compare different epileptic populations in order to compare how photic stimulation effects these different groups, including patients with genetic generalized epilepsy (GGE) who are photosensitive and GGE patients who are not photosensitive, as well as patients with Psychogenic Non-Epileptic Seizures (PNES) and healthy controls.Kalitzin et al. (2002) and Parra et al. (2003) found increased excitability in photosensitive individuals during stimulation using a measure of phase clustering called relative phase clustering index (rPCI).An extended study used rPCI to indicate greater probability of epileptic seizures (Kalitzin et al., 2005).In another study, higher visually evoked potential (VEP) amplitudes indicated greater responses to photic stimulation in individuals with photosensitive epilepsy compared to controls, especially in occipital areas (Wilkins et al., 2004).In fact, individuals with photosensitive epilepsy exhibited stronger suppression of visual perception after TMS over the visual cortex, and had lower phosphene thresholds compared to healthy controls, indicative of increased excitability of the visual cortex.
Studies have used power spectrum measures to indicate changes in excitability.For instance, reduced inhibition of alpha rhythm generating networks was associated with photosensitive epilepsy (Vaudano et al., 2017), and stimulation via flashed chromatic gratings induced alpha desynchronisation in posterior electrodes (Haigh et al., 2018).Changes in power spectral density (PSD) in lower frequencies have been demonstrated to predict seizures from EEG recordings (Chu et al., 2017), which is consistent with expectations of increased autocorrelation and has a theoretical relationship with the power spectrum (Chatfield and Xing, 2019), such that we may expect increased low frequency power and longrange temporal associations.Therefore, increased spatial similarity could also be a candidate EWS, that could be measurable via synchrony measures like phase clustering and phase synchrony.This is supported by findings of increased phase synchronization in the gamma frequency range in EEG, being associated with increased neuronal excitability in epilepsy patients (Meisel et al., 2015).We extend these analyses by using the PPR and photic stimulation from the perspective of dynamical systems theory.We study the transition to PPR as a critical transition, and expect to observe early warning signs prior to that transition.As mentioned, we expect increased response to perturbation, variance, and autocorrelation in the brain.The goal of this research is to measure cortical excitability in epilepsy patients via a perturbation response.In this case, the perturbation response is non-invasive and in the form of photic stimulation, which is a commonly used investigation in routine EEG studies.We hypothesize that we can use photic stimulation as a perturbation method for observing increased cortical excitability before PPR.We expect to observe early warning signs of a PPR transition that are identifiable via changes to a particular metric's distribution over time.We analyse distributions of statistical moments of the EEG timeseries as a datadriven method for inferring changes to system dynamics, because the underlying dynamics of brain activity are not observable.These include changes in variance and autocorrelation function widths of brain activity, as well as in phase-based synchrony measures.In particular, we expect the greatest increases to be demonstrated in the photosensitive GGE (PPR) group, just prior to a PPR, where they likely experience heightened cortical excitability.

Methods . Participants
The study was conducted using data obtained retrospectively from databases of the Department of Neurosciences, St. Vincent's Hospital, Melbourne, over the period May 2012-August 2022.This included selected EEG recordings with photic stimulation from 10 patients with genetic generalized epilepsy (GGE) and photoparoxysmal response (PPR) used in a previous study (Seneviratne et al., 2016), 10 GGE patients without PPR, and 10 patients with psychogenic non-epileptic seizures (PNES).Diagnosis of GGE was established using International League Against Epilepsy (ILAE) criteria (Scheffer et al., 2017).The diagnosis of PNES was confirmed with the consensus of at least two epilepsy specialists following video-EEG monitoring using criteria from a previous study (Seneviratne et al., 2010).Clinical data for all patients is summarized in Table 1, with a patient-specific breakdown in Supplementary Table 1.The use of this data for study was approved by the human research ethics committee of St Vincent's Hospital Melbourne, where the data was obtained.

. . Healthy controls
Healthy control EEG data was obtained from an open source repository, originally described and analyzed in a study by Torkamani-Azar et al. (2020).Control data included ten (six female) healthy adults aged 22-45 (M = 30.25,SD = 6.95).For each control subject, we analyzed 2.5 min of eyes-closed, resting-state EEG.The eyes-closed condition was chosen for consistency since our analysis focuses on stimulation epochs in which participants were instructed to close their eyes.More information regarding recording method is provided in the original study.Electrodes were selected to match our data, and sampling rate was conserved.The same EEG pre-processing and measure calculations were used for this control dataset.

. Electroencephalography recording
Electroencephalography (EEG) data were acquired with a Compumedics Grael 4K-EEG system (Compumedics Ltd., Melbourne, Australia) or Siesta 802 in the case of ambulatory recordings.Scalp electrodes were placed according to the 10-20 international system.All EEGs were recorded with a bandwidth of 0.15-120 Hz, with a sampling rate of 256 Hz, including provocation techniques such as hyperventilation and intermittent photic stimulation (IPS) (as illustrated in Figure 1).Photo-paroxysmal response (PPR) was identified through visual inspection of the recordings, defined as epileptiform discharges triggered by IPS visible on the EEG.For clinical reasons, a number of participants, EEG recordings were also continued for a full day recording using the ambulatory recording system, sometimes capturing spontaneous seizures.
For GGE (PPR) participants, the ordering of stimulation frequencies often differed after their first PPR was observed; stimulation frequencies were sometimes repeated in an attempt to elicit PPR again, and sometimes alternated to high frequencies to determine an upper sensitivity range, as is convention (Kasteleijn-Nolst Trenité et al., 2012).It is not clear whether stimulation order effects exist whether rest periods of approximately 30 seconds truly account for refractory effects of PPR (Brausch and Ferguson, 1965).Therefore, in an attempt to control for these confounding variables, we only analyzed data before the first PPR.For GGE (no PPR) and PNES participants, stimulation frequency ordering was the same.In order to account for confounding effects of frequency-dependent responses, we only included stimulation epochs up to 15Hz in GGE (no PPR) and PNES groups in our analyses, as we observed PPR at or before 15Hz in all GGE (PPR) participants.

. Pre-processing
Analyses were performed by importing raw data files from Profusion EEG 5 into MATLAB R2021a.
For variance and autocorrelation measures, we applied bandpass filters between 0.5 and 70 Hz to reduce high frequency artifact and remove DC drift.For phase-based measures, this filtering step was ommitted, as more narrow bandpass filtering was later applied to calculate the measure.
To normalize the data, we focused on making the majority of the signal comparable between different channels with varying ranges, allowing tail values but without being overly skewed by the presence of artifacts distorting the re-normalized range.Therefore, we applied z-score normalization using the native MATLAB normalize function, using the medianiqr argument to center the data at 0 and ensure an inter-quartile range of 1, with a window length of two seconds.In GGE (PPR) participants, since only data before their first PPR was analyzed, normalization was not affected by the high amplitude spiking during PPR.

. Measures
To capture the expected rapid dynamics, measures generally used small windows with maximal overlap, and are described in detail below.

. . Variance and autocorrelation function width
Variance and autocorrelation function width (ACFW) were calculated for each non-reference EEG channel (n = 19).The variance of the signal was calculated for each channel with a moving window of two seconds (512 samples), with overlapping steps of 1 sample using the native MATLAB movvar function.
To measure the autocorrelation function width (ACFW), using the same aforementioned window and overlap, the autocorrelation function was first calculated using MATLAB's xcorr function, and the width in samples was measured at the half prominence of the function peak, centered at it's maximum.

. . Phase-based synchrony measures
To calculate phase-based synchrony measures, we obtain the phase-angle time series φ(t) (the instantaneous phase of each sample in the time-series) by constructing an analytic signal where the real component is the original time-series x(t), and the complex component is the Hilbert transform of the real signal H[x(t)].In other words, the analytic signal is the product of the envelope and instantaneous phase of the original timeseries, (t) . (1) Note that this only holds if the original signal satisfies Bedrosian's Product Thoerem, which requires that the signal is narrow-banded using a band-pass filter (Ktonas and Papp, 1980).To do this, we filtered the normalized data using an elliptic filter, which is a zero-phase filter that circumvents phase shifting (Gustafsson, 1996).Specifically, a filter was constructed using MATLAB's designfilt function, using the bandpassiir filter, with a filter order of 10, pass-band frequency limits according to the upper and lower limits of standard EEG frequency bands (delta: 1.5-4 Hz; theta: 4-8 Hz; alpha: 8-12 Hz; beta: 13-30 Hz; gamma: 30-70 Hz), a pass-band ripple of 3, and stop-band attenuation of 40.

. . . Relative phase
Cosine of the Relative Phase (CRP): Honari et al. ( 2021) is calculated by taking the cosine of the instantaneous phase difference between two channels: cos( φ(t)), where φ x (t) − φ y (t) obtained from the original channel timeseries x(t) and y(t) with the method previously explained for a given frequency band.This produces a measure ranging from −1 to 1, where magnitudes close to 1 denote high synchrony.Values close to −1 also arguably demonstrate a highly synchronous event, since this is produced when there is a phase difference of exactly π , known as an antiphase relationship (Honari et al., 2021).Measures were considered for each combination of non-reference channels (n = 19), resulting in n(n − 1)/2 = 171 unique combinations of channels.

. . . Global and posterior synchrony
We also calculated the average synchrony between sets of channels using the average phase vector for a given sample.For N oscillators, with instantaneous phases θ j for j = 1, 2, ...N, their position can be represented on the unit circle by e iθ j .Then, the complex number z is the average of their positions,

FIGURE
Schematic of the photic stimulation procedure.Stimulation epochs last for s with a s rest in between.The frequency of subsequent stimulations increase until a Photo-Paroxysmal Response (PPR) is triggered, and this activating frequency is designated the lower bound.Then, the technician starts from a high frequency and decreases the frequency of subsequent stimulations until another PPR is triggered, and this frequency is designated the upper bound.These lower and upper bounds are indicative of the range of frequencies that trigger a PPR response specific to that patient.

FIGURE
Each of channels are shown on the unit circle as a blue dot, representing the instantaneous phase for a given time sample.The average phase is indicated by the direction of the red vector on the unit circle, and it's length corresponds to their similarity, where a perfect alignment of phase would produce a length of , while a uniform distribution around the circle would be . (2) We found the magnitude of z using MATLAB's native abs function; this produces a value from 0 to 1, where 0 represents no synchrony, or a uniform distribution, and 1 represents all phases being equal over all channels for a given time sample (see Figure 2 for an illustration).We calculated the magnitude of the average phase vector over all 19 channels and designate this global synchrony, as well as over six posterior channels only (P3, P4, T5, T6, O1, O2) and designate this posterior synchrony.
The magnitude of the mean phase vector can also be calculated over an epoch of time.For N oscillators over an epoch containing S samples, the instantaneous phases θ j,k for channels j = 1, 2, ...N and time samples k = 1, 2, ...S can be represented on the unit circle by e iθ j,k .Then, the complex number ẑ is the average of their positions, Again, the magnitude of this vector produces a value from 0 to 1, where 0 represents no synchrony, or a uniform distribution, and 1 results if all phases were equal between all channels over all time samples in the epoch.

. Analysis strategies and statistics . . Clinical characteristics
Differences in ages between the four groups were evaluated using a one-way ANOVA (using the anova1 function from the MATLAB Statistics and Machine Learning Toolbox: SML).If significance was observed, follow-up multiple comparisons were evaluated using Tukey's honestly significant difference (HSD) procedure (by passing the stats struct returned from the ANOVA into multcompare from the SML Toolbox).
Categorical clinical characteristics such as sex were examined between all four groups using Pearson's chi-squared test of independence (using the crosstab function from SML).Differences in ASM, syndrome, and seizure type(s) were also evaluated using this test between the two GGE groups (with and without PPR) only, since PNES participants were not taking ASMs and were currently exhibiting non-epileptic seizures.Categories of ASM and seizure type(s) were considered distinct for each unique combination, as grouped in Table 1, since this allows for distinction of combinations that may be clinically significant.
The aforementioned tests were deemed significant with a pvalue less than 0.05.Note that Tukey's HSD procedure adjusts p-values for pair-wise multiple comparisons automatically.

. . The photo-paroxysmal response
We compared time-frequency spectrograms of PPR events and spontaneous seizures for the four participants (P1, P2, P5, P6) for which both events were available.To generate the spectra, unfiltered epochs were fed to MATLAB's spectrogram with window size of 256 samples, overlap of 128, and transform (nfft) length of 512 samples.Other PPR characteristics including PPR count, duration, and lower and upper frequency bounds are also measured for each participant.

. . Measures
To investigate the hypothesis that photic stimulation perturbs brain dynamics toward a more excitable state that is measurable via changes to timeseries measures, we compared distributions of a number of measures during photic stimulation epochs compared to a control distribution.Substantial deviations from the control distribution suggests a shift in dynamics of the underlying system.
Photic stimulation distributions contained epochs before the first PPR in the GGE (PPR) group, or before the 15 Hz stimulation in the GGE (no PPR) and PNES groups.We also examined an additional pre-PPR epoch distribution in the GGE (PPR) group only, containing the closest full stimulation epoch before (i.e., not overlapping) the first PPR event.
Since sufficient quality resting state data from the photic stimulation recordings were not available, we compared our data to distributions of healthy controls during eyes-closed, resting EEG.To create the control distributions, we randomly sampled 10s epochs from the full measure timeseries, with equal representation from the individuals of the cohort.This constructed distribution emulated the process of pooling 10s stimulation epochs in the epileptic cohort groups.
Differences between pairs of group distributions were evaluated using permutation-based hypothesis tests, which are non-parametric tests that do not make assumptions about the distribution of the data (Ernst, 2004).This process involves creating a joint distribution of two considered groups, and performing the test statistic of interest on random permutations of the joint distribution (where each permutation is equally as likely as another) to construct a distribution of test statistics (n = 10, 000) that represents the null result distribution.The test statistic between the original separate distributions is then compared to the null distribution, where the p-value is the fraction of permutation tests that are at least as extreme as that test statistic.A pseudocount (addition of one to the fraction numerator) is included to avoid p-values of zero (where the test statistic is never surpassed by permutation values).This is possible because we used a subset of the possible permutation set (though subsets of sufficient size have been shown to be accurate) (Knijnenburg et al., 2009).In our case, this single test statistic value was the average value of a test distribution constructed from permutation subsets of the two considered distributions.From the central-limit theorem, this test statistic distribution will be normal with the mean at its center.We created confidence intervals around this statistic by extracting the 2.5th and 97.5th percentile of the test distribution.
We used the Kolmogorov-Smirnov test statistic for continuous measures, which evaluates the difference between empirical cumulative distribution functions (CDF), using rank-based comparisons.This is appropriate considering variance and synchrony measures demonstrated non-normal distributions through inspection of q-q plots.The Kolmogorov-Smirnov test identifies small differences in large datasets as significant even when differences seem trivial (Riffenburgh, 2012); therefore, permutation samples of size 500 were taken, in accordance with sample size recommendations (Abadie, 2002).Chi-squared tests were used as the test statistic for auto-correlation function width, which is a discrete measure.
We also investigated whether measures of synchrony would change with photic stimulation.The relative phase measure contained significant degrees of freedom (unique channel pairs = 171, frequency bands = 5, resulting in 855 potential comparisons), so we investigated whether factor analysis could reduce the parameter space by identifying channels that demonstrated the greatest variation.Factor analyses of the GGE (PPR) stimulation epochs (using factoran from SML, with 5 component factors to generate the maximum likelihood estimate) were carried out.This analysis was performed for each frequency band separately to investigate whether variation was also frequency-band dependent.
To further investigate whether measures varied by channel, and whether this is apparent in posterior channels reflecting underlying posterior pathophysiology, we compared posterior and global synchrony distributions to the corresponding control distributions.Differences were evaluated using the same permutation based hypothesis tests described above, using the two-sample Kolmogorov-Smirnov test as the test statistic, since synchrony values are continuous.
Bonferroni correction, considered the most conservative multiple-comparison correction, was carried out (the standard pvalue was divided by 56: a total of 12 comparisons for variance and autocorrelation distributions, and 5 frequency bands x 4 distributions x 2 sets of channel groups in the synchrony analyses), requiring a p-value less than 0.00096 for significance.
We also investigated whether ASM, syndrome, and seizure type differed between the two GGE (with and without PPR) groups.We did not reject the null hypothesis of independence between group and ASM [χ 2 (5) = 3.81, p = 0.58], syndrome [χ 2 (4) = 3.53, p = 0.47], or seizure type [χ 2 (5) = 1.87, p = 0.87]; therefore, these groups did not differ in the spread in these variables (see Table 1 for a breakdown of categories).
. The photo-paroxysmal response PPR featured spike-wave morphology in all channels and lasted between 0.6-28 seconds (see Figure 3 for an example).As can be seen in Table 2, GGE patients with PPR exhibited individualized responses to photic stimulation, including differences in their photosensitivity ranges as well as PPR length.
In four participants, spontaneous seizures were observed hours after photic stimulation.We compared the transition to PPR events and spontaneous seizures by their time-varying frequency spectra by visual inspection (Figure 4).Though spectra varied according to the individual, we observed changes to the distributions of power per frequency bin at the transition, including increased power especially in delta and theta frequencies.The spectral dynamics of PPR and seizure transitions appear similar.
We inspected spectrograms over participants' stimulation trains and often observed increased power in the base frequency of stimulation, as well as its harmonics (see Supplementary Figure 2).We observed this in all participant groups during photic stimulation.Comparison is difficult due to varying data quality and different stimulation frequency ordering in the photosensitive group.The visibility of these bursts of increased power varied substantially by participant, though we report that it was most visible in the occipital channels.Interestingly, they are found to exist even during PPR (see Supplementary Figure 2A).

. Measures
We examined whether measures could track individualized responses to perturbations by examining the timeseries of measure values from GGE (PPR) participants as a PPR approaches.

. . Variance and autocorrelation
From inspecting individual stimulation trials, we observed an increase in variance during stimulation epochs close to a PPR (see Figure 5 and Supplementary Figure 3 for further examples).
We analyse whether these measures may track cortical excitability using group-wise comparisons.To do this, we compared the distributions of variance and autocorrelation function widths (ACFW) during photic stimulation epochs between groups.Values were pooled over all 19 electrode channels and the 10 participants in each of the four groups.We compared the cumulative distribution functions of these measures to resting healthy controls, and observed a greater relative proportion of high values in all groups during photic stimulation (see Figure 6).Significant differences in empirical cumulative distribution functions (CDF) was found compared to controls in all groups during stimulation, as well as in the the GGE (PPR) group just prior to their first PPR ("pre-PPR"), where cortical excitability was hypothesized to be greatest (see Table 3 for a breakdown of statistical test outputs).There was no significant differences between the photosensitive GGE (PPR) group and the non-photosenstive epilepsy groups.

. . Synchrony . . . Relative phase
To investigate whether changes in synchrony may indicate increased cortical excitability, we constructed phase diagrams of relative phase values in stimulation epochs leading to a PPR.In some cases, we observed greater clustering in epochs close to a PPR, indicating that a particular channel pair was becoming more synchronous (marked also by longer mean vector lengths, see Figure 7).However, clustering magnitudes varied considerably between individuals and depending on combinations of factors like channel pair, frequency band, and time to PPR, suggesting considerable individual variation.Magnitudes tended to correlate positively between frequency bands, though lower frequency bands (delta, theta, alpha) tended to have ranges extending to higher values, compared to beta and gamma, which tended to have lower values overall, likely related to different degrees of freedom from bandwidth ranges.
To investigate whether certain channels may demonstrate greater change in synchrony, a factor analysis was performed on relative phase values during stimulation from the GGE (PPR) participants.In this analysis, the data from each EEG channel are assumed to depend on a linear combination of latent components, which attempt to capture the variance of the data along their axes, (similar to Principal Component Analysis).Each channel's factor loadings for each component axis indicate the amount of variance in the data that was captured along that component.The three dominant component factors separated a subset of posterior channels, which demonstrated the greatest variance along these factor axes (see Figure 8).This result was relatively consistent over all frequency bands (see Supplementary Figure 4).

. . . Global and posterior synchrony
To investigate changes in synchrony in broader brain regions, we compared global and posterior synchrony in groups during photic stimulation compared to controls.Significant differences in empirical cumulative distribution functions (CDF) was found compared to controls in some frequency bands and groups for both global and posterior synchrony (see Figure 9 for an illustration of the CDFs, and Table 4 for a breakdown of statistical test outputs).For global synchrony, stimulation groups differed significantly to controls in some frequency bands, including gamma in all groups, with an overall reduced probability of high coherence values (see Figure 9).In comparison, posterior distributions showed increased probability of high coherence values in stimulation groups compared to controls, reaching significance for all frequency bands except the GGE (PPR) stimulation gamma distribution.There was no significant differences between the pre-PPR distribution and stimulation distribution for GGE (PPR) group for global or posterior synchrony.
Distinction between each group's stimulation distribution is not seperable from inspection of CDF distributions, though the presence of differences relative to controls across stimulation groups suggests a group-independent effect of the photic stimulation on relative global and posterior synchrony, that is not specific to a particular group.

Summary
The primary finding of this study was that the response to perturbation by photic stimulation was observable via changes in EEG measures.Specifically, distributions of variance and autocorrelation function width (ACFW) shifted to higher values just before photo-paroxysmal response (PPR), indicating greater cortical excitability close to seizure transitions, as hypothesized.In all epilepsy groups, phase coherence distributions differed to controls.As expected, posterior channels demonstrated the greatest variation in synchrony.Greater changes were sometimes observed in epochs close to PPR for some individuals, but there was no significant difference at a group level.This is not surprising considering the patient-specific nature of responses, and we suggest future research apply individualized analyses utilizing repeated-trials.We have demonstrated measurable changes of excitability across different groups of epilepsy patients using photic stimulation, particularly photosensitive patients in the lead up to seizure transitions consistent with previous literature (Scheffer et al., 2009;Meisel et al., 2015;Maturana et al., 2020).Therefore, we argue that PPR is a valid proxy for studying seizure transitions.This human model is a safer alternative to inducing seizures, and can be reliably induced via photic stimulation.We argue this model has great potential for developing new diagnostic methods and evaluating treatments for epilepsy. .

Relevance and significance of findings
At a group level, we observed greater proportions of higher variance and autocorrelation distributions in all group distributions during photic stimulation compared to resting controls (see Figure 6).Individually, patients demonstrated considerable variation in response values, particularly in the photosensitive group, motivating our group-level analysis to understand overall trends.Values tended to be highest in the photosensitive group, compared to the non-photosensitive generalized epilepsy group, and patients with psychogenic non-epileptic seizures displayed the lowest values; however, this difference was not found to be significant in the group-level analysis.As expected, the photosensitive group demonstrated further increases just before photo-paroxysmal response.
We interpret the shift in distributions as a data-driven illustration of underlying dynamical regime changes close to seizure transitions.These observations are consistent with predictions of increased system summary statistics before state transitions in dynamical systems, together with increased responses to perturbation (Scheffer et al., 2009).These early warning signs of transitions have been observed before seizure transitions (McSharry et al., 2003;Maturana et al., 2020).Though, our data differs in that it reflects changes to the actively evoked perturbation response, instead of just passive observations of change.Studies have elicited active perturbation using invasive EEG for focal epilepsy and observed similar warning signs (Meisel et al., 2015(Meisel et al., , 2016)), while we report it here using non-invasive perturbation methods in generalized epilepsy.We argue that the lead-up to PPR demonstrates similar changes in dynamics to seizure transitions, making it a good proxy for studying seizure transitions.Further, PPR demonstrates characteristic spike-wave morphology, and similar spectral content to spontaneous seizures visible in spectrograms.Increased power was evident in a broad range of frequencies including the 2-5Hz activity characteristic of seizure activity (Kasteleijn-Nolst Trenité et al., 2012).
We speculate that the increased variance before PPR, is associated with an increased firing activity of neurons in synchrony, considered to underpin seizure activity.Increased autocorrelation width suggests longer lasting increases to self-similarity (Scheffer et al., 2009), changes to the periodicity in the signal, decreased dimensionality (Elger and Lehnertz, 1998), and possibly increased  variance.We also understand that for cyclo-stationary signals, the autocorrelation function is equal to the absolute value of the square of its Fourier transform, in which case increased width of the major peak may reflect changes to the ratio of power in various frequency ranges.We observed distinct changes in power at the frequency of stimulation and their harmonics, which supports this explanation (see Supplementary Figure 2).It could also reflect a breakdown in the signals' stationarity, which is to be expected during a change in dynamical regime.In sum, it is difficult to disentangle the exact mechanism  Significance is evaluated through permutation (n = 10, 000) hypothesis tests of the test statistic where '*' denotes significance (adjusted for multiple comparisons using bonferroni correction).NS denotes no significant difference.D is the Kolmogorov-Smirnov test statistic.χ 2 is the Chi-squared test statistic.95%CI obtained through the 2.5th and 97.5th percentiles of the test statistic distribution constructed through permutation re-sampling (n = 10, 000).behind these changes, although they do reflect changes in brain dynamics.
Compared to resting controls, all groups demonstrated significantly reduced gamma coherence when averaged over all channels (see Figure 9) but increased coherence when averaged over only the posterior subgroup in nearly all frequency bands (see Figure 8).We also found the majority of change in phase clustering (cosine of the relative phase) values occurred in the  Significance is evaluated through permutation (n = 10, 000) hypothesis tests of the test statistic where "*" denotes significance (adjusted for multiple comparisons using bonferroni correction).NS denotes no significant difference.D is the Kolmogorov-Smirnov test statistic.95%CI obtained through the 2.5th and 97.5th percentiles of the test statistic distribution constructed through permutation re-sampling (n = 10, 000).
posterior subgroup via a factor analysis.This evidence of increased posterior synchrony is well established in the literature (Porciatti et al., 2000;Kalitzin et al., 2002;Parra et al., 2003;Fisher et al., 2022), including in healthy individuals (Salchow et al., 2016), suggesting recruitment of the visual system that is arguably enhanced to flicker in particular.It is likely that subsequent aberrant propagation of posterior synchronization due to abnormal gain control results in PPR (Fisher et al., 2022).However, we were not able to find substantial group differences between photosensitive and nonphotosensitive distributions during stimulation, or between the photosensitive group prior to PPR and earlier epochs, which suggests the measure was not able to indicate increased cortical excitability, as was found in ECoG data (Meisel et al., 2015).We did, however, observe substantial individual variation, and multiple examples of increased coherence close to PPR in some individuals, channels, and frequency bands (see Figure 7).Therefore, we encourage future work to conduct individualized, frequencyspecific evaluations of multiple stimulation trials to elucidate enhanced responses in photosensitive individuals.

. Limitations
In this retrospective study, there were several limitations.This study was deliberately conceived using retrospective data from routine EEG databases to investigate its limits as a non-invasive perturbation dataset.Due to the retrospective study design, there were several limitations.Because of this, we did not have appropriate resting state recordings, instead compared data to resting controls.However, ideally, each individual would be compared to their own resting or baseline distribution to ascertain the true response to perturbation.Another issue was the considerable within-group variability that we observed, suggesting that the population distribution for each epilepsy group was itself heterogeneous and would overlap each other considerably.We observed frequency-dependent differences in the responses perturbation, consistent with previous research that generally found marked responses in the 10-20Hz range (Kasteleijn-Nolst Trenité et al., 2012), and in harmonic frequencies of the PPR-evoking frequency (activation frequency) (Parra et al., 2003).However, it would not have been appropriate to perform an individualized frequency-dependent analysis, since only a single trial was available for each participant, and because of confounders.This included ordering effects due to variation in stimulation frequency order (which could be mitigated through randomization), possible refractory effects from PPR events, and individualized responses to perturbation.Our work highlights the limitations of current data acquisition methods in routine EEG that severely limits individualized analysis.

. Future work
The applications for photic stimulation are three-fold: (1) it can be used to diagnose photosensitive epilepsy; (2) it can be used as a therapeutic tool for epilepsy patients to monitor treatment efficacy; and (3) it can be used as a perturbation paradigm for studying brain dynamics.
Photic stimulation may continue to be used for its current diagnostic purposes for epilepsy.However, we recommend for research purposes, that multiple trials within a single session may be necessary to obtain sensitive results.Multiple trials will also enable the construction of robust, patient-specific profiles of frequencydependent responses to perturbations.These profiles would be clinically significant in that changes to the initial measurements could be used to evaluate the efficacy of medications, which are understood to alter cortical excitability as their mode of action (Premoli et al., 2017).They could also be used to evaluate the safety of medication withdrawal, and create standardized clinical guidelines, by indicating whether drug removal has shifted patients to regimes of higher cortical excitability and, therefore, seizure risk (Meisel et al., 2015(Meisel et al., , 2016)), and monitor disease course in the long-term.Also, it should be noted that observing measurable responses to perturbation from photic stimulation is not limited to photosensitive individuals.We observed responses to perturbation in all groups, regardless of epilepsy diagnosis (GGE or PNES), syndrome, seizure type, age, sex, or ASM.We conclude that the response to perturbation may be observable in the general epilepsy population, suggesting that the photic stimulation paradigm offers a great potential framework for studying seizure dynamics at large.Though, due to documented effects of these clinical factors (Fisher et al., 2022), family history (De Kovel et al., 2010), as well as the effect of ASM action on responses (French et al., 2014), we encourage further investigation of these sub-groups with sufficient power to observe group-level differences.
We wish to emphasize the potential for photic stimulation to be used as a perturbation paradigm for studying general brain dynamics.Previous, studies have already shown that greater responses to perturbation are associated with elevated cortical excitability relating to seizure risk in focal epilepsy (Maturana et al., 2020), with cohorts of patients undergoing pre-surgical evaluations (Freestone et al., 2011;Bergey et al., 2015) and with severe treatment-resistant epilepsy (Cook et al., 2013;Oderiz et al., 2019).Like the implanted electrical stimulation in these studies, photic stimulation enacts a perturbation.Compared to these surgical methods, it can be considered a much safer method for studying the transition to seizure, as well as compared to nonsurgical alternatives like seizure provocation during hospitalization via medication withdrawal and sleep deprivation.In fact, this non-invasive method could potentially be used at home, only requiring an EEG headset, appropriate software, and a monitor for administering the appropriate photic stimulation.Changes to an individual's baseline at any frequency could indicate meaningful differences to cortical excitability and seizure risk, regardless of the type of epilepsy.Also, as the response to photic stimulation differs in other diseases including migraine (De Tommaso et al., 2013), schizophrenia (Portnova and Maslennikova, 2023), alzheimer's disease and mild cognitive impairment (Fisher et al., 2022;Kim et al., 2023), its relevance as a research paradigm to studying wider diseases of brain structure or function is significant, in particular cortico-thalamic circuitry dysregulation in which many brain diseases are implicated (Guillery and Sherman, 2002;Sarnthein et al., 2005).Furthermore, photic stimulation could be used to observe changes in cortical excitability in healthy individuals under initiating factors (Kasteleijn-Nolst Trenité et al., 2012;Fisher et al., 2022) like alcohol withdrawal and sleep deprivation, to provide further insight to acute symptomatic seizures of people without epilepsy and human brain dynamics in general.

. Conclusion
We found that the response to photic stimulation was a valid non-invasive perturbation to brain dynamics.This paper represents photic stimulation data from photosensitive generalized epilepsy patients, patients with GGE who are not photosensitive, and patients with psychogenic non-epileptic seizures, compared to resting controls.At a group level, consistent with the hypothesis that PPR induces a change in brain state, we observed clear changes in group signal statistics including the CDF of the variance and autocorrelation.We further investigated changes in synchrony measures, and could not find consistent group-level patterns of change that correlated with the transition.However, we did find posterior increases in synchrony consistent with past research.In conclusion, EEG time-series analysis was able to capture changes in cortical excitability during active perturbation, via biomarkers that track brain dynamics before state transitions.We find that current data acquisition methods limit individualized analysis, and encourage future to extend current routine EEG protocol to include the randomization of frequency order, acquisition of baseline or resting periods, and importantly, multiple stimulation trials.Photic stimulation demonstrates a non-invasive human method for studying seizure transitions, and potentially monitoring and evaluating treatments in epilepsy and other brain diseases.

FIGURE
FIGUREPhoto-paroxysmal response (PPR).The figure shows a typical PPR lasting approximately seconds, visible as spike-wave morphology.Lines are di erent channels as indicated on the y-axis, and time is indicated in seconds on the x-axis.

FIGURE
FIGURETime-varying frequency spectra for PPR (top row) and spontaneous seizure (bottom row) transitions which occur at the s mark in each subplot.Columns are di erent participants, and data is sourced from the "O " channel (similar spectra were observed in other channels).Time in seconds is shown on the x-axis, frequencies up to Hz are shown on the y-axis, and power per frequency bin for a given window of samples is represented as a color, where values correspond to a colormap (see legend), where higher values are bright red, and lowest values are deep blue.In general, increased power can be observed from the midway ( s) point, especially in lower frequencies.Note that vertical lines at Hz is a power line artifact.

FIGURE
FIGUREThe first subplot shows the EEG trace (blue line) of participant in the leadup to a PPR, visible as high amplitude spiking activity after the marked PPR onset boundary (vertical red line).Consecutive red dots mark seconds of photic stimulation for the frequency (labeled below them).Below, the corresponding variance of the signal is displayed, using the shared time axis (y-axis), where di erent electrode channels are displayed as di erent colored lines.Observe the increasing variance close to the boundary.See Supplementary Figure for further participant examples.

FIGURE
FIGURECumulative distribution functions (CDF) for (A) variance (left plots) and (B) autocorrelation function width (ACFW; right plots), where the x-axis is the values for each measure, and the y-axis is the probability of observing a corresponding x-value and all values less than it.Each colored line represents a di erent group's distribution, pooled over all participants and channels for selected epochs (refer to legend).Right-shifted lines indicate distributions with larger proportions of high values of the variance and ACFW.Observe how all groups show right-shifts compared to the control (gray) line.Also observe how, for the GGE (red) group, their distributions shift further to the right, in the epoch just before PPR (dotted line labeled as "pre-PPR"), compared to earlier stimulation epochs [solid line labeled as "GGE (PPR) stim"].

FIGURE
FIGURE Chronological phase histograms during P 's stimulation train.Each plot represents s of binned relative phase di erence values between O and Fz for epochs leading to the first PPR.Each bin has a direction from to pi around the unit circle, and a length indicating the relative frequency of values, ranging from to .The radius of each circle plot varies to assist visualization, and the the mean relative phase amplitude (time-varying synchrony) is indicated by a black line whose length is labeled above each plot and also ranges from to , where indicates greater clustering of the values in the plot.Each row corresponds to relative phase values within di erent frequency bands (see legend).

FIGURE
FIGURE (A) Biplot of factor loadings (blue vectors) on the three dominant components (axes), generated from a factor analysis of cosine of the relative phase values during photic stimulation, pooled over GGE (PPR) participants.Each channel's factor loadings for each component axis indicate the amount of variance in the data that was captured along that component.Channels with the highest factor loadings, or in other words demonstrated the greatest variance, along the component axes are P , T , O , O , P , T (labeled in red text).This suggests that the greatest change in the data was demonstrated by this subset of channels.Similar results were obtained when repeated in other frequency bands (see Supplementary Figure ).(B) Subgraph of posterior channels.

FIGURE
FIGUREEmpirical cumulative distribution functions (CDF) for coherence in di erent groups (di erent colored lines, see legend), within di erent frequency bands (columns).The x-axis is coherence calculated as the mean phase vector magnitude over a set of channels.The y-axis is the cumulative probability for the corresponding x-values.From visual inspection, there is evidence for decreased global synchrony during stimulation in the first row where left-shifts are observed compared to controls (black line), most evident in the gamma band.Increased posterior synchrony is evident in the second row in all frequency bands, with right shifted CDFs for groups during stimulation (colored lines) compared to controls (black line).See Table for a breakdown of significant di erences.(A) CDFs for Global Coherence (mean magnitude over all channels).(B) CDFs for Posterior Coherence (mean magnitude over posterior channels only).
TABLE Summary of clinical data for participants in all groups.

TABLE Participant -
specific information about PPR events.