Nine patients with drug-resistant epilepsy (4 women; mean (± SEM) age = 37 ± 3 years, see Table S1 for details) were enrolled at the Motol Epilepsy Center in Prague. They underwent iEEG monitoring for precise localization of the epileptic seizure onset zone before surgery. All the patients signed informed consent to participate, and the study was approved by the Ethics Committee of Motol University Hospital. All the patients had normal or corrected-to-normal vision.

The experimental task aimed to investigate the neural processes involved in memory-guided actions based on different types of visual information. Each trial consisted of fixation, encoding, delay, and action/recall phases. The trial began with a jittered fixation period of 1.9–2.1 s, displaying a white cross on a dark gray screen (Fig. 1A), followed by a 2 s encoding period with a central red cross and two types of objects: either two identical circles ("same" condition) or a square and a triangle ("different" condition). The positions of the objects were randomized for each trial, with one object always positioned closer to the cross than the other, maintaining a consistent distance ratio of 1.5. In the 'different' condition, the triangle was closer in 50% of the trials, and the square was closer in the other 50%. The participant's task was to remember which object was closer to the cross: its position in the 'same' condition and both its position and identity in the 'different' condition. Since the cross was always centered on the screen, corresponding to the participant's midsagittal plane, we assume that distance estimation was performed using egocentric coordinates.

Task design, iEEG channel locations, and exemplary responses. A Experimental design of the delayed action task. Each trial includes a jittered fixation period (1.9 s–2.1 s), encoding (2 s), a jittered delay (3.9 s–4.1 s), and an action/recall phase (2 s). Participants have to remember the object's position (same condition) or both its position and identity (different condition) and reach for it with a joystick after the delay. In the different condition, a two alternative question about the object's identity is followed with a 2 s time-window to answer. Participants respond using the gamepad: green 'A' for a triangle or red 'B' for a square. B The projection of all implanted channels in 3 brain regions in 9 patients (369 channels) from both hemispheres on the standard MNI brain template (adult MNI-ICBM152 head model, [26]; http://www.ucl.ac.uk/dot-hub). Left, sagittal view; right, coronal view. IPL, inferior parietal lobule; VTC, ventral temporal cortex; HIP, hippocampus. C Trial-averaged spectrograms of exemplary individual channels in IPL, VTC, and HIP for the same and different conditions, 0—the start of the encoding phase; vertical lines mark the boundaries of the task periods, and only the last 3.9 s of delay are plotted. Prominent changes during the delay occur at low-frequency power.

After the encoding phase, a white fixation cross reappeared on a dark gray screen during a delay phase, jittered between 3.9 s and 4.1 s to prevent anticipatory responses. Similar jitter ranges have been effectively used in previous studies [9, 10]. Although this range may not entirely eliminate anticipation, it was selected for task feasibility and participants’ comfort, given the clinical settings of our experiment. Following this, a green cross appeared for 2 s during the action/recall phase, prompting participants to respond. Using a gamepad joystick (Xbox wireless controller), they moved the green cross to the memorized position of the nearest object. The subjects made a straight reach from the center of the screen to the memorized object and returned to the center of the screen. A response was considered correct if the minimum distance of their trajectory was within a square area defined as the X,Y coordinates of the correct object plus or minus half the minimum distance between objects over all trials (10.5% of the screen width). Reaction time was measured from the start of the action phase to the time of reaching the correct area. In the ‘different' condition, an additional two-alternative question about the object's identity appeared for 2 s at the end of the trial. Participants responded by pressing a button on the gamepad: the green 'A' button for a triangle or the red 'B' button for a square. The total trial duration, including the pre-trial period, ranged from 9.8 s to 10.2 s for the 'same' condition and from 11.8 s to 12.2 s for the 'different' condition.

The task was counterbalanced with 160 trials equally divided between the 'same' and 'different' conditions. Trials were grouped into blocks of 10, each assigned to one condition, with a subject-controlled break between blocks. The order of blocks was counterbalanced so that the average order of 'same' and 'different' blocks was the same. At the beginning of each block, participants received simple on-screen instructions about the upcoming condition. The experiment began with a brief presentation explaining the task, followed by a training session with shortened blocks of five trials per condition, with feedback provided after each trial. The experiment also included 160 immediate trials, in which participants reached for the object immediately without any delay, again equally divided between the 'same' and 'different' conditions. As the main focus of this study was the analysis of the delayed trials, in particular the functional connectivity between the two streams during the maintenance phase, we report the results of the analysis of the immediate trials in the Supplementary Materials.

The experimental task was implemented in PsychoPy3 v2020.1.3 [27] and run on a 15.6-inch TFT notebook monitor with a refresh rate of 60 Hz. The task and iEEG recording were synchronized by TTL pulses sent to a trigger port of the iEEG recording system at the start of each trial.

The iEEG was recorded with stereotactically implanted multi-contact electrodes, often also referred to as stereo-EEG. Recording sites were selected on an individual basis, strictly according to the medical requirements of the pre-surgical evaluation of epileptic zones, without reference to the present study. Eleven to fifteen semi-rigid electrodes were implanted intracerebrally per patient, depending on the suspected origin of their seizures. Each electrode had a diameter of 0.8 mm and consisted of 8–18 contacts, 2 mm long, and 1.5 mm apart (DIXI Medical Instruments). The iEEG signal was recorded with medical amplifiers (Quantum, NeuroWorks), sampled at 2048 Hz, and later downsampled to 512 Hz to reduce the computational load. The reference contact for each patient was located in the white matter. Post-implantation CT, co-registered with pre-implantation MRI, was used to identify the positions of electrode contacts in each patient [28]. The anatomical locations of the electrode contacts were labeled by an experienced neurologist. The contact positions were then normalized to the Montreal Neurological Institute (MNI) space using standard Statistical Parametric Mapping algorithms (SPM 12) for group-level visualizations.

The downsampled iEEG recordings were first visually inspected, and any bad electrode contacts with obvious artifacts were discarded. Contacts identified as being in the seizure onset zone or heterotopic cortex were also excluded. Bipolar derivations were calculated between adjacent contacts to suppress contributions from distant neuronal assemblies and enhance spatial specificity. Importantly, bipolar re-referencing reduces the effects of volume conduction and common reference artifacts, which can lead to spurious correlations in connectivity analyses due to shared signals [23, 29]. Bipolar iEEG signals were visualized at the center between two contacts. Henceforth, we refer to the bipolar contact pair simply as a 'channel'. When the channel was derived from two contacts in a different brain structure, we labeled it with the structure with a larger unilateral response. Line noise was removed from the iEEG signal with a notch filter (4th order Butterworth stop-band filter of 1 Hz width centered at 50 Hz and harmonics, zero phase shift). Data processing and analysis were performed in MatLab R2018a.

For further analysis, we selected channels located in three brain regions of interest (containing a total of 369 channels in 9 patients, see Table 1 and Fig. 1B): the IPL, including supramarginal and angular gyri and the adjacent lateral wall of the intraparietal sulcus; the VTC, including the inferior temporal, lingual, fusiform and parahippocampal gyri; and the hippocampus (HIP).

Time-frequency analysis was applied for linearly increasing frequencies between 2 and 120 Hz, with a resolution of 1 Hz bins, using a technique called the filter-Hilbert method. Similar to our previous studies [30,31,32], we estimated the power change using the following procedure. First, we band-pass filtered the entire recording dataset (third-order Butterworth filter, zero phase shift) in consecutive non-overlapping 1 Hz frequency bands. For each band, we extracted the amplitude envelope using the Hilbert transform; the obtained envelope was downsampled to 64 Hz, resulting in a time resolution of 15.625 ms. The envelope of each band was normalized by dividing it by its mean value over the entire recording session, channel-wise for each frequency band, effectively whitening the broad frequency band and compensating for the 1/f frequency decay of EEG signals [33]. Then, similar to [10] and based on visual inspection of time-frequency spectrograms, we extracted power for two non-overlapping frequency bands: theta (2–7 Hz) and alpha (8–13 Hz). To do this, the original 1 Hz bands in these theta and alpha ranges were averaged and multiplied by 100 to obtain a single time series of theta and alpha power for each channel expressed as percentages of the mean value. These two-time series signals were then divided into epochs. To ensure that the epochs had the same duration, we left out from the epochs the jittered period from the beginning of the delay phase (0.0s–0.2s). Each epoch was then divided into the following five task periods: (1) baseline (0.5 s, end of the fixation period), (2) encoding (2.0 s), (3) first half of delay (1.9 s, referred to as Delay 1 in the text), (4) second half of delay (1.9 s, referred to as Delay 2) and (5) recall (2.0 s). In the result, the epochs were from −500 ms to 7800 ms relative to the stimulus onset - the beginning of the encoding phase.

The mean of the −500 ms to 0 ms of the fixation period (i.e. the baseline) was subtracted from each epoch to remove signal changes independent of the respective stimulus. We excluded from further analysis epochs in each channel containing interictal epileptiform discharges, which were identified by a spike detector implemented in MatLab [34]. In addition, trials with incorrect behavioral responses—incorrect joystick responses for the object's position and incorrect button responses for the object's identity—were excluded from the iEEG analysis.

To test for the significance of the mean power in a given frequency band relative to the baseline, we used a Wilcoxon signed-rank test, with false discovery rate (FDR) correction across the time samples and channels, similar to our previous studies [30,31,32]. For each channel, we compared the average power for all trials of the respective condition during the pre-stimulus interval (–500 ms–0 ms before encoding) with all the time points during the post-stimulus period (0 ms–7800 ms). As a conservative estimate, we used a sliding window of six samples (93.75 ms) with the highest p-value. If there was a significant difference at any time point relative to the baseline for a selected condition, the channel was considered 'active' in the given frequency band.

Then, for all active channels in a given frequency band (theta or alpha), we averaged the power in each task period: encoding, first and second halves of delay, and recall. These averaged values were submitted to a linear mixed effects model (LMEM) with task period and condition and their interaction as fixed effects and channel and patient as random effects (P <0.008, Bonferroni correction for six models - three brain regions and two frequency bands):

$$ }\sim + }*} + \, \left( |}} \right) \, + \, \left( |}:}} \right) $$

To evaluate the functional connectivity between channel pairs, we calculated the PLV [35]. This was done using a multitaper frequency transformation with two tapers based on the Fourier transform, covering a frequency range of 2–20 Hz with a resolution of 1 Hz, as implemented in the FieldTrip toolbox [36], and similar to other studies [17, 24].

The PLV between channels i and j is defined as:

$$ }_ \left( f \right) = \frac\left| ^ \frac \left( f \right) \cdot \left( \left( f \right)} \right)^ }} \left( f \right)} \right| \cdot \left| \left( f \right)} \right|}} } \right| $$

where N is the number of trials, X(f) is the Fourier transform of signal x(t) (1.9 s task period), and (∙)* represents the complex conjugate. Note that this formula is equivalent to the PLV calculation presented by Lachaux et al. [35] and should not be confused with the coherence equation.

Phase differences were calculated for each channel pair (i,j) between IPL and VTC, and between IPL and HIP, using the spectra of the encoding, the two halves of the delay, and the recall periods to quantify inter-electrode phase coupling. To determine statistically significant differences in the PLV during the encoding, both parts of the delay, and the recall, we compared them to the baseline (1.9 s fixation interval) using permutation statistics with cluster correction. The null distribution was created by calculating the PLV differences between the baseline and each task period (the encoding, the first half of the delay (1.9 s, Delay 1), the second half of the delay (1.9 s, Delay 2), and the recall (0.5 s, the part before the response was made)) for randomized data. Trials were randomly assigned to these two periods, and the PLV differences between them were calculated and repeated 200 times. This randomization preserves any frequency-dependent biases due to temporal smoothness present in both task and baseline periods. Only those frequency bins with PLV above the 95th percentile threshold of the null distribution were considered statistically significant [17, 24]. Cluster correction was applied to account for multiple comparisons across the 2–20 Hz frequency range (P <0.05 was used to obtain null-hypothesis clusters). A channel pair was considered significant if its PLV was significant at any frequency point in the 2–20 Hz range after cluster correction.

We found that different channel pairs exhibited PLV significance at various frequency points within the 2–20 Hz range. We aggregated all channel pairs significant in any of the four task periods (the encoding, two parts of the delay, and recall) from all patients and calculated the ratio of significant pairs (Sig.P.Ratio) in each 1-Hz frequency bin. This was done by dividing the number of significant pairs between two brain regions by the total number of channel pairs between these regions, for all subjects and in each 1-Hz bin independently [37]. This Sig.P.Ratio was calculated separately for IPL-VTC and IPL-HIP pairs. To identify frequency ranges at which the proportion of channel pairs was significantly above the chance level, we applied the binomial test to these Sig.P.Ratio values (one-sided test, P <0.05, FDR-corrected across frequency bins). The median Sig.P.Ratio of each region pair across all frequency bins was used as the chance level for the binomial test [37].

We then averaged the PLV for each task period (baseline, encoding, two parts of the delay, and recall) for each frequency range with a Sig.P.Ratio significantly higher than the median chance level. These averaged values were then analyzed using the LMEM with task period, condition, and their interaction as fixed effects, and channel pair and patient as random effects (P <0.0167, Bonferroni correction for three models, each for one identified frequency range; see Results):

$$ }\sim + }*} + \, \left( |}} \right) \, + \, \left( |}:}} \right) $$

To evaluate the direction of information flow between IPL and VTC, and between IPL and HIP, we applied spectral non-parametric Granger causality (GC) analysis as a measure of directed functional connectivity [23]. Similar to another study [24], we downsampled the iEEG signals to a sampling rate of 40 Hz and calculated the spectral non-parametric GC in the 1–20 Hz frequency range. We transformed the signals into the frequency domain using the multitaper method with two Hann tapers and zero-padding to 20 s to reduce spectral leakage and control frequency smoothing.

In the frequency domain, GC quantifies the extent to which the spectral content of the source signal contributes to the spectral content of the target signal at each frequency. We used a non-parametric spectral approach to compute GC at a given interval time [23]. In this method, the spectral transfer function H(f) and the noise covariance matrix Σ are estimated from the cross-spectral density matrix obtained from the Fourier transforms of the data. The total power spectral density matrix S(f) is computed as S(f) = H(f)ΣH*(f), where H*(f) denotes the complex conjugate transpose of H(f). The GC from signal Y to signal X at frequency f is then calculated as:

$$ }_ \left( f \right) = \ln \left. \left( f \right)}}_ \left( f \right)}}} \right.} \right) $$

Where Sxx(f) is the total power of signal X, and S̃xx(f) is the intrinsic power of X at frequency f when the influence of Y is removed.

We applied the GC analysis using the implementation in the FieldTrip toolbox [36]. GC was calculated for four task periods—baseline (fixation), encoding, delay 1, and delay 2 (each 1.9 s in duration)—for each condition (same and different objects). We did not analyze the recall phase because it was too short for reliable spectral GC analysis; behavioral responses started ~500 ms–600 ms after the beginning of the recall phase, making the effective time window insufficient for accurate spectral estimation. We focused our analysis on the channel pairs between IPL and VTC, and between IPL and HIP, which were found to be significant in the PLV analysis, as these pairs already suggested meaningful interactions. By running the GC analysis only on these pairs, we concentrated on the most relevant connections.

To identify the frequency ranges where the difference in directional influence between the two directions in VTC-IPL connections was the largest, we aggregated the GC data across all subjects, channel pairs, task periods, and conditions by averaging, and computing the difference in GC values (net Granger, similar to [24, 38]) between the two directions for each frequency bin:

$$\Delta }_} \to }}} = }_} \to }}} - }_} \to }}} $$

and similarly for HIP and IPL:

$$\Delta }_} \to }}} = }_} \to }}} - }_} \to }}} $$

We identified the frequency ranges where the net GC showed maximal differences. We then averaged the net GC over these frequency ranges for each task period and condition. These averaged values were analyzed using an LMEM with task period, condition, and their interaction as fixed effects, and channel pair nested within participants as random effects:

$$ }\sim + }*} + \, \left( |}} \right) \, + \, \left( |}:}} \right) $$