2.1 EEG Data recordings
Participants
The real EEG data presented here were obtained from three healthy subjects, males, right-handed, aged 28 for two of them and 30. They were instructed to perform five minutes of RS with their eyes closed, to relax, and to refrain from thinking (Fig. 1a). All subjects were volunteer to participate in this study and provided written informed consent. This study was conducted with the approval of the local ethics committee of Rennes University Hospital, France (CPP ID-RCB: 2019-A00608-49; approval number: 19.03.08.63626).
Pre-processing of HD-EEG recordings
EEG signals were recorded using a 256-channel HD-EEG system (Electrical Geodesic Inc., EGI) with a sampling frequency of 1000 Hz. The channels were placed according to the 10–20 and 10 − 5 systems, and all channel impedances were checked before the beginning of the task. Recordings were allowed if impedances were below 10 kΩ. A total of 199 electrodes were retained, the majority of the jaw and neck electrodes being removed.
The EEG preprocessing was performed manually using the Brainstorm toolbox (Tadel et al. 2011). The preprocessing steps included the removal of DC offsets, the application of a notch filter at 50 Hz, and the implementation of a band-pass filter from 2.5 to 100 Hz (to remove very slow oscillations that sometimes contaminate the baseline). The signals were visually inspected, and any channels of poor quality were removed before being interpolated using Brainstorm's default parameters.
Finally, the recordings were epoched into ten blocks of 30 seconds for the analysis. The reference epoch was selected based on the presence of minimal artefacts, which were identified as being present in the original data. Subsequently, a visual inspection was performed to identify epochs with excessive residual noise leading to exclusion for further analysis.
Inverse Solution
The EEG inverse problem involves the estimation of the unknown parameters of the dipolar source S(t) at the cortical level, including position, orientation, and magnitude, from the measured EEG signals X(t) at the scalp level. In this study, we used the Desikan-Killiany parcellation with 68 regions of interest (ROIs) (Desikan et al. 2006) to locate cortical sources and constrained their orientation normally to the cortical surface (Dale and Sereno 1993). Therefore, the EEG inverse problem was reduced to estimating the magnitude of sources: (1)
The inverse matrix W was computed with the weighted minimum norm estimate (wMNE) method (Lin et al. 2006; Allouch et al. 2022). This method compensates for the tendency of the classical minimum norm estimate to favor weak and surface sources. The following equation was used: (2)
$$W={BG}^{T} {\left(GB{G}^{T}+C\right)}^{-1}$$
where B is the diagonal weighting matrix (inversely proportional to the norm of lead field vectors), G is the lead field matrix computed with the boundary elements method (BEM) in Brainstrom (Gramfort et al. 2010), α is the regularization parameter (based on signal to noise ratio (SNR): α = 1/SNR, and C is the noise covariance matrix (calculated from our 700 ms pre-stimulus baseline). The Brainstorm toolbox was employed to compute wMNE. The SNR was set to 3, and the depth weighting value to 0.5 (default values).
Regions of interest
The inverse problem was performed to reconstruct RS activity on the 68 scouts. Then, we filtered epochs of reconstructed RS activity at source’s level in the frequency bands of interest (FBOIs) for RS, which are alpha waves between 8-12Hz and beta-gamma waves between 12-45Hz, as previously described (Koenig et al. 2002; Khanna et al. 2015). The filtered RS activity at source level was then used to determine the amplitude of the oscillations. Brainstorm was used to measure the dynamics of oscillations amplitude in FBOIs at the source level from RS activity, by creating a timelapse of these dynamics. Then ImageJ software, version 1.53K (Schneider et al. 2012) was employed to superimpose images of the maximal intensity of FBOIs, from the timelapse created. The ROIs were selected based on the highest incidence of activity on the image from ImageJ, for both FBOIs (see Fig. 1a).
2.2 Simulated EEG with eCOALIA model
Whole brain model based on layered neocortical NMM
The simulated cortical-level activity is generated using the updated version of the computational model implemented in the COALIA model (Bensaid et al. 2019).This model has been shown to produce realistic EEG when compared to real EEG recordings obtained in humans during awake and deep sleep, in terms of the morphology, spectral content, and topographical voltage distribution.
More specifically, at the local node level, each NMM is composed of interconnected neural subpopulations of glutamatergic neurons (pyramidal cells [PYR]) and GABAergic neurons (such as perisomatic targeting PV and basal and apical dendritic targeting SST, VIP, and NGFC) (Köksal-Ersöz et al. 2022; Wendling et al. 2024). At the mesoscale level, the model incorporates physiologically based synaptic kinetics (fast/slow) and circuity between these neuronal subpopulations (Tremblay et al. 2016; Wendling et al. 2024).
At the whole-brain level, the large-scale model is constructed with 82 interconnected NMMs (66 NMMs for the cortical areas and 16 NMMs for the subcortical areas). The structural connectivity matrices was obtained from DTI data collected in healthy subjects and reported in the Human Connectome Project (Van Essen et al. 2013). The template brain morphology (Holmes et al. 1998) is employed to spatially distribute neural masses over the neocortex. Each NMM simulates the local field potential (LFP) of one brain region, wherein the activity is presumed to be homogeneous. The time delay was organized as a matrix, with elements representing the Cartesian distance between cortical NMMs divided by the mean velocity of action potential travel, with a mean velocity of 7.5 m/s. The NMMs are connected via long-distance pyramidal glutamatergic projections (Fig. 1b). We used the Desikan-Killiany atlas to define the cortical ROIs (Desikan et al. 2006).
It is important to note that this study aimed to test methods using a previously validated model, rather than introducing or validating a new model. For a more detailed description of the eCOALIA model, readers are directed to published studies (Bensaid et al. 2019; Köksal-Ersöz et al. 2022; Wendling et al. 2024).
Simulation
The eCOALIA model generates simulated EEG with distinct brain oscillations as a function of the relation between the excitation mediated by the PYR subpopulation and the inhibition mediated by the GABAergic subpopulations, or inter-regional interactions. Alpha-rhythms can be created through the PYR SST loop, beta/gamma-rhythms through PYR-PV loop, and delta-rhythms through increased thalamocortical connectivity and disinhibition through VIP-SST connections.
To simulate a RS EEG, the NMMs were modified in the ROI defined in real EEG following resolution of the inverse problem (Fig. 1c). This allowed the FBOI to be replaced in a subject-specific manner. The real EEG was used to create a FBOI near the same frequency peak in the model and the recorded data. Subsequently, we solved the forward problem using the lead field matrix G computed with BEM (Gramfort et al. 2010), to obtain simulated EEG (Fig. 1c). The model was modified until the results of this first step were deemed satisfactory (through visual inspection).
Measuring the similarity between real and simulated EEG
Comparison in time and frequency domains
One of the aims of our study was to create a realistic EEG. To achieve this, we compared the spatial correspondence of the FBOI using a method that involved a direct comparison of simulated and real EEG. To this end, and to ensure that the chosen metric was efficient in classifying EEGs as more or less different, a series of six different simulations were generated, varying in their degree of similarity. The simulations included RS EEG, RS EEG with modified occipito-temporal rhythm from reference EEG (mRS), RS EEG with suprasylvian interictal spikes (IS), RS EEG with modified parameter spikes (mIS), ictal-like activity called “ictal discharge” (ID), and modified ictal-like activity (mID; Fig. 2a). All epileptic simulations were obtained from previous studies conducted by our team (Wendling et al. 2024). Thus, the simulated EEGs were classified according to their theoretical degree of similarity, by simple visual inspection and based on the difference from the chosen simulation parameters. For example, RS was found to be significantly different from IS and completely different from ID, and all modified versions of these EEGs (i.e., mRS, mIS, mID) were found to be slightly different from their original versions.
The comparison was conducted between a 30-second epoch of "reference" EEG (Fig. 2b) and a "compared" EEG that was cut into 10-second epochs used as a sliding window. When the "reference" EEG was simulated, the epoch could be less than 30 seconds due to an artefact (DC shift) created at the beginning of the simulation, during approximately half a second. Additionally, for ID and mID, epochs were trimmed to encompass only the 10-second discharge pattern, as these EEGs did not reflect a stationary state in the same manner as (m)RS and (m)IS.
Comparisons were conducted within a limited frequency range, spanning from 2.5 Hz (to avoid ocular movement artefacts real data) to 45 Hz (to avoid differences resulting from the notch filter in the real data) as this is the interval most often used in clinical practice.
As no clear metrics for comparison have been described in the literature, several metrics were tested. In fact, when comparing simulated and real EEGs, no consistent and repeatable metric was found. We chose metrics already described in connectivity studies. We applied these different metrics in both frequency and time domains. We tested phase locking value (PLV) (Lachaux et al. 1999), phase lag index (PLI) (Stam et al. 2007), cross correlation (CC), Jensen-Shannon distance (JS). In addition, JS and CC were applied in the frequency domain after using Fast Fourier Transformation (referred to as FFT JS and FFT CC), as the time lag between EEG signals is known to degrade the results (Zunino et al. 2022; Eilts and Putze 2022). Comparisons between metrics were made by comparing simulated EEGs with decreased theoretical similarity as follows: ID vs. mID, RS vs. mRS, RS vs. mIS, IS vs. mID, RS vs. mID (box plots of the different metrics are shown in Supplementary Fig. S1). We applied metrics between the "reference" and "compared" EEGs as follows: the "compared" EEG was shifted through the "reference" EEG where each step is corresponding to one sample. This process yielded a specific number of values for each bipolar derivation. For instance, when using a 30-second "reference" EEG and a 10-second "compared" EEG, with a sample frequency of 1000 Hz, we computed 20 000 values. This was done without padding EEG for comparison. We then averaged the resulting values for each line to obtain one value for each bipolar derivation. Figure 2B illustrates the method employed for comparisons conducted in the frequency domain (i.e., FFT JS and FFT CC). Following the extraction of the FFT for each line, a similarity metric was employed. Subsequently, the obtained values were averaged for each bipolar derivation.
The similarity metric that appeared to be the most efficient in terms of theoretical similarity of simulated EEGs was FFT CC. The CC is a linear method to compare different signals, defined by (3):
$${r}^{2}=\frac{{cov}^{2}\left(x,y\right)}{var\left(x\right).var\left(y\right)}$$
where r is the correlation coefficient between x and y. The CC obtained ranges between − 1 (signals that are totally opposite) and 1 (signals that are identical; 0 means that the signals are totally different). CC is known to be sensitive to the lag between phases. To avoid this, we chose to apply CC in the frequency domain, and so the latter metric was computed by applying the FFT method for each bipolar derivation.
We represented box plots of FFT CC obtained from analysis in Fig. 2c. Box plots were realized by using R software (R Core Team 2020).