Tackling the challenges of group network inference from intracranial EEG data

Introduction: Intracranial EEG (iEEG) data is a powerful way to map brain function, characterized by high temporal and spatial resolution, allowing the study of interactions among neuronal populations that orchestrate cognitive processing. However, the statistical inference and analysis of brain networks using iEEG data faces many challenges related to its sparse brain coverage, and its inhomogeneity across patients.

Methods: We review these challenges and develop a methodological pipeline for estimation of network structure not obtainable from any single patient, illustrated on the inference of the interaction among visual streams using a dataset of 27 human iEEG recordings from a visual experiment employing visual scene stimuli. 100 ms sliding window and multiple band-pass filtered signals are used to provide temporal and spectral resolution. For the connectivity analysis we showcase two connectivity measures reflecting different types of interaction between regions of interest (ROI): Phase Locking Value as a symmetric measure of synchrony, and Directed Transfer Function—asymmetric measure describing causal interaction. For each two channels, initial uncorrected significance testing at p < 0.05 for every time-frequency point is carried out by comparison of the data-derived connectivity to a baseline surrogate-based null distribution, providing a binary time-frequency connectivity map. For each ROI pair, a connectivity density map is obtained by averaging across all pairs of channels spanning them, effectively agglomerating data across relevant channels and subjects. Finally, the difference of the mean map value after and before the stimulation is compared to the same statistic in surrogate data to assess link significance.

Results: The analysis confirmed the function of the parieto-medial temporal pathway, mediating visuospatial information between dorsal and ventral visual streams during visual scene analysis. Moreover, we observed the anterior hippocampal connectivity with more posterior areas in the medial temporal lobe, and found the reciprocal information flow between early processing areas and medial place area.

Discussion: To summarize, we developed an approach for estimating network connectivity, dealing with the challenge of sparse individual coverage of intracranial EEG electrodes. Its application provided new insights into the interaction between the dorsal and ventral visual streams, one of the iconic dualities in human cognition.

1. Introduction

Brain network dynamics give rise to the unique ability of (human) brain to carry out complex processing of external inputs and control complex behavioral outputs. To better understand the structure and dynamics of these networks, the tools of network neuroscience (Bullmore and Sporns, 2009) are being extensively developed and applied, ranging from efficient data acquisition and preprocessing to estimation and further analysis of brain connectivity networks. Multiple methods are available for neuroscientists to gain various types of insight into brain connectivity. Structural imaging, in particular diffusion weighted magnetic resonance imaging allows the estimation of the structural connectivity of the human brain using the methods of tractography, while various methods of functional brain imaging, including in particular the functional magnetic resonance imaging (fMRI) and electroencephalography (EEG), enable the estimation of the patterns of statistical dependencies between remote neurophysiological events, known as the functional connectivity (Friston et al., 1993). While the functional connectivity provides in principle only information that the activity of some brain areas is related, a whole plethora of advanced methods for the estimation of the directed, cause-effect relationships between brain areas has been developed, commonly summarized under the term effective connectivity (Friston, 1994).

Inferring the effective connectivity can be seen as providing the ultimate answer to the question of information flows in the brain, however effective connectivity estimation is hindered by a range of technical challenges. Successfully tackling these challenges would provide an appropriate causal description of the brain network, that would lend itself to further processing using the tools of complex network theory.

In this work, we discuss a range of such challenges and possible solutions, presenting an analysis of a particularly problematic dataset using a pipeline implementing a specific combination of solutions. Note that a combinatorial explosion of the number of possible pipelines makes a comparison of all alternatives not feasible—thus the presentation serves as an example rather than a comparison of all tools available. The challenging task that we select for the demonstration is the estimation of causal network from intracranial EEG (iEEG). Intracranial EEG recording as such is a powerful tool for the discovery of brain function. This type of measurement is characterized by both high temporal and spatial resolution, a combination which is impossible for other, non-invasive, neuroimaging techniques. However, the analysis of the iEEG data faces a set of challenges mainly because of in-homogeneous and spatially sparse data acquisition - the measurement sites, corresponding to the implanted electrodes, are few and far between, thus not covering all the areas of interest in any one particular subject. Moreover, the positioning of the implantation differs widely between subjects, thus complicating a group-level analysis. In our paper, we develop the methodological pipeline for estimation of such group network structure, that can't be estimated from any single patient. As an example, we show the connectivity analysis of the visual processing streams.

1.1. Visual streams conundrum

Large body of evidence suggests that visual information is processed in the brain in two separate streams originating from the visual cortex: the dorsal and ventral streams in parietal and temporal lobes respectively. These streams were originally characterized in terms of the “where” and “what” distinction (Ungerleider and Mishkin, 1982), suggesting the role of the dorsal stream in spatial position and motion coding and of the ventral stream in object identity coding. Examination of patients suffering visual agnosia after bilateral ventral lesions, like D.F., in combination with their unimpaired spatial-motoric abilities lead to an alternative model (Goodale and Milner, 1992). In this model, characterized as “how” and “what” model, the dorsal and ventral streams are distinguished not by the type of information processed, but by the results of their visual processing, which is motoric action or conscious perception. According to this view, the dorsal stream processes the visual information to enable movements in space like grasping an object, while the ventral stream processes the same information for perception of objects, their identity and also their spatial relationship. Several brain areas in the ventral stream are specialized for recognition of various classes of visual percepts. For instance, damage to an area in lingual gyrus leads to landmark agnosia (Aguirre and D'Esposito, 1999) and a more anterior parahippocampal place area (PPA) in the posterior parahippocampal/anterior lingual region is one of the brain regions specifically responding to spatial scenes (Epstein and Kanwisher, 1998).

In contrast to the well established role of ventral stream, there are still controversies regarding the dorsal stream function and its connections with the ventral stream. The dorsal stream actually seems to divide into three parallel pathways serving different functions. The most ventral of them, the parieto-medial temporal pathway, has been associated with spatial information transfer from dorsal to ventral stream. The existence of this pathway was suggested based on monkey and human brain anatomy (Kravitz et al., 2011) and is also supported by human resting-state functional connectivity (Boccia et al., 2017). Its functional properties have been suggested based on the role of the included structures like retrosplenial cortex (Byrne et al., 2007). Besides this pathway, the dorsal and ventral streams seem to be interconnected several times, with the dorsal stream providing attention focus and spatial information (Cloutman, 2013). On the other hand, the occipito-temporal cortex in the ventral stream seems to provide object identity information to the dorsal stream as was only indirectly suggested from a human fMRI experiment (Kristensen et al., 2016). Direct evidence about these connections and their functional role is however still scarce.

2. Causal network inference

Many challenges of causal network inference are generic and relevant irrespective of the method used to obtain the time series. For instance, in principle, any estimation method is negatively affected by the presence of measurement noise and estimation from only short observations, issues present to various extent in all neuroimaging methods. Indeed, standard methods such as the common linear implementation of Granger causal analysis, that rely on assumptions concerning the underlying process (namely, that it is a stationary linear autoregressive process), may provide relative advantage in sensitivity for short time series. This advantage against nonlinear methods with less or no assumptions may be practically relevant even for systems that, in principle, have nonlinear character—such as the brain or the Earth's climate (Hlinka et al., 2013).

However, the character of neuroimaging data poses additional challenges beyond its nonlinearity, ranging from some very general to those quite specific to iEEG data. A first commonly faced challenge is the multivariate, rather than bivariate, character of the data. Indeed, in principle one can estimate the causal interactions separately for each pair of variables. The downside of such straightforward approach is the computational demands, but more importantly, the fact that the estimated interaction between any pair of variables might be confounded by not accounting for the effect of other variables in the role of mediators or common drivers. If one believes that the other variables observed may play such a role, it is advantageous to fit a joint multivariate model that takes (all) these variables into account. Of course, such a model has a large number of parameters, and might thus call for even longer time series samples, and pose additional computational problems, particularly when nonparametric methods, such as conditional mutual information, are in use. Then, an alternative approach lies in the stepwise inference of a smaller set of causally relevant “parent” variables for each target variables—see Kořenek and Hlinka (2020) for a systematic comparison of a set of such state-of-art methods.

A second challenge lies in the nonlinear character of the brain dynamics, which calls for the application of nonlinear, potential nonparametric methods, such as the conditional mutual information. However, a range of prior studies have shown, that the nonlinear/nongaussian character of neuroimaging data may be relatively weak, and thus the use of linear methods might be justified in terms of the trade-off between the lost information and higher speed and general reliability of linear approaches—see a discussion for fMRI (Hartman et al., 2011; Hlinka et al., 2011) and electroencephalography (Blinowska and Malinowski, 1991). Therefore, the choice of linear methods is generally suitable.

A third aspect that must be considered with brain signals, particularly with EEG, is that the dynamics and interactions take place across a range of time scales. In fact, specific EEG activity frequencies have been shown to be related to specific types of cognitive or arousal states. Therefore, the use of frequency-specific estimates of causal interactions is common. The downside of this approach is that the inferred object further increases dimension—causality is resolved not only across all pairs of observed channel signals, but also across frequencies. We shall discuss the consequences of this aspect for overall statistical inference later.

Albeit the multivariate character of the observed data allows (and calls for) taking into account indirect causality mediated/driven by any of the observed variables, the effect of any other, directly unobserved, variables is not accounted for in our methodology, and this omission may lead to inference of spurious causal links. It is important to highlight that this is a generic effect that affects most of the available methods for causal inference. While some methods have been proposed to suppress the effect of such unobserved (latent) variables (such as those based on the utilization of time reversal Haufe et al., 2012; Winkler et al., 2016), their assumptions are not universally fulfilled, and particularly in multivariate setting can lead to wrong conclusions (Korenek and Hlinka, 2021). Although other methods are available that attempt to suppress the effect of latent variables and non-neuronal interactions, such as the volume conduction in tissue, the problem is far from ultimately solved, and one should always interpret the results of causal analysis with caution in terms of the potential role of any latent intervening variables.

Another complication, which is not entirely specific, but quite common in EEG analysis, is that the systems is deemed to be substantially nonstationary, and thus the structure of causal interactions may meaningfully change in time. A typical engineering solution to such situation is to estimate the causal structure using a (potentially overlapping) sliding window for a sequence of time periods. Note that an alternative exists in using variants of Kalman filtering (Kalman, 1960)—this avoids the problem of selecting the window parameters; however still requires selection of other intrinsic parameters.

The requirement of obtaining time resolved causality leads to estimation from short time series; to compensate for the loss of sample size, in the case of task-based iEEG, one may typically utilize data across many trials for the estimation of the parameters of the underlying vector autoregressive process. The consequence of this is that only a single estimated causal structure is estimated across all trials. The trial dimension can not be thus used in straightforward way further statistical assessment of the inferred causal links, i.e., by evaluating the distribution across trials with respect to the posited null hypotheses. However, alternative ways for establishing empirical null distribution are available, as discussed below, that are based on the use of bootstrapping, surrogate data or permutation testing.

Indeed, the evaluation of the statistical significance of causal analysis results obtained by most methods of causal inference relies on some construction of empirical null distribution of selected test statistics. The distribution of the causal indices under the hypothesis of no interaction is, however, in most cases not known analytically, see e.g., the approach of surrogate testing proposed in Kaminski et al. (2001).

Note that with task-based electrophysiology data, we might not necessarily focus on the structure of causal interactions per se, but on the deviation of the structure from the resting baseline. The rationale is two-fold. First, it can be reasonably expected that in the resting state, most parts of brain are functionally connected and the change in the task scenario might contribute only (potentially minor) alteration to this ongoing process of interactions (Biswal et al., 2010). Second, due to the latent variables as well as technical confounds and volume conduction, the causal structure estimated in task is expected to be biased toward detecting a range of spurious interactions. However, the same set of spurious sources are to affect the estimate of interactions outside the task condition of interest. Therefore, some form of correction for these effects by subtraction or comparison with respect to the baseline estimate of causal structure may be warranted to identify the causal interactions specific to the task condition of interest.

A specific problem mentioned earlier is the high dimensionality of the estimated object. In particular, the causal interaction is estimated between any two channels (taking all other channels into account), for each frequency and time window. Even with a valid statistical significance test, such situation invites a multitude of false positive results due to the multiple testing problem involved. Indeed, one may attempt to apply one of the methods for multiple testing correction, that controls the family-wise error rate (FWE), or, as a less strict option, the false discovery rate (FDR) (see Hochberg 1988; Benjamini and Hochberg 1995). However, the application of these methods may decrease the overall statistical power, lead to false negatives and also create extreme computational demand, as p-value estimates with high precision call for the generation of high number of permutations/surrogates to generate the null distributions. A potential solution to this problem lies in some dimension reduction of the test statistics prior to its statistical evaluation.

Finally, there are challenges that are highly specific to the iEEG data. In particular, each subject has their own set of implanted electrode locations at different brain positions, that is moreover extremely sparse and inhomogeneous. How can one in such a situation reliably estimate the structure of brain network interactions across the whole group of subjects? A reasonable solution seems to be to (at least temporarily) sacrifice some of the spatial resolution (that is anyway somehow artificial with only sparse coverage), and aggregate the data within pre-defined generic brain regions. Note that the data can, with advantage, be aggregated (by averaging or with another approach) also later in the processing, such as at the level of connections, rather than already at the level of the measured signals.

Note that for task activation data, such as evoked potentials, an alternative approach would be to merge all the channels into one joint dataset, to obtain a single “super-subject,” while ignoring that the channels actually come from different subjects. However, in case of analysis of connectivity this approach is not directly applicable, because the analysis of connectivity between channels from different subjects might be detrimentally affected by the lack of synchronization of processes across subjects. Also, for some purposes, the interpretation of the brain network calls for some coarse-graining to the level of reproducible and interpretable brain regions, and thus rendering detailed info on individual channels not necessary.

The collating the brain connections into generic region pairs can in principle entail both integrating the data within subjects, as well as across subjects. Given the heterogeneous sampling over space and subjects, different regions pairs will have contributions from different number of channels across different subjects. Together with additional problem of dependence between channels (and channels pairs) within subjects, this poses another problem of statistical inference concerning the connections observed at the inter-regional level. We believe that the solution to this challenge lies in the construction of a surrogate data that respects the within-subject channel-pair dependence as well the contribution of particular subjects to the particular region connections. Note that an alternative solution has also been applied in the literature, that aimed for analysis of only those pairs of regions, that are sampled (at least by a single pair of channels) in some minimal number of subjects, that allows statistical inference across subjects (of some subject-level test statistics, typically obtained by summing/averaging within subject) (see, for example, Burke et al., 2013, 2014; Bastin et al., 2017).

Below, we develop a pipeline for the analysis of the interaction of the two cortical visual processing streams, applied to a recorded iEEG task-based dataset used by Vlcek et al. (2020). In this analysis, we illustrate the considerations outlined above, that need to be taken into account while analyzing the causal structure of such complex neuroimaging datasets.

3. Data

The current study used a dataset including 27 patients, who voluntarily participated in a visual recognition experiment, as described by Vlcek et al. (2020). All patients had normal or corrected to normal vision. This study was approved by the Ethics Committee of Motol University Hospital and all patients gave their informed consent to participate. Intracranial EEG was recorded in the Motol Epilepsy Center in Prague in pharmacoresistant epilepsy patients using depth platinum electrodes (Dixi Medical), designed for local field potential (LFP) recordings. Electrodes were 0.8 mm in diameter, with 5–18 contacts of 2 mm length, separated by insulator segments of 1.5 mm length. Up to 128 contacts were recorded based on the medical requirements. Wideband signals (0.05–512 Hz) were recorded with a sample rate of 512, 2, 048, or 8, 000 Hz, and later subsampled to 512 Hz. All contacts were localized in the standard Montreal Neurological Institute (MNI) space.

Because the implantation is decided by clinical requirements, it naturally includes many areas not involved in visual processing. Thus, to limit computational costs and improve statistical power, only those channels with significant responses to demonstrated pictures were analyzed (see Vlcek et al., 2020 for statistical analysis details). Moreover, for the purposes of connectivity analysis, only those subjects, who had implanted at least two of the initially defined seven regions of interest (ROIs) relevant for the cognitive tasks studied (described in Section 4), were included in the analysis. This provided a subset of 15 subjects, the data of which were included in the final analysis (8 women, median age 30 years, range 17–48 years, education level: 1 primary school, 11 secondary school, and 3 college-level education).

Each time series was normalized by the baseline standard deviation. As described in Vlcek et al. (2020), from the EEG recording of the whole experiment, bipolar derivations were computed between adjacent electrode contacts to suppress contributions from distant neuronal assemblies and further assumed that the bipolar EEG signals can be considered as originating from a cortical volume centered between the two contacts. Thus, any two contacts adjacent spatially on the same implanted electrode, after omission of rejected contacts, were used to define a bipolar channel by subtracting their signals. We refer to the bipolar contact pairs further as “channels.” To eliminate the non-stationarity of the time series in the sense of turning the signal to wide-sense stationary, we computed the differences in time, i.e., from the signal y_t we computed $y_t^{\prime }=y_t-y_{t-1}$ (Barnett and Seth, 2014). The wide-sense stationarity requires the mean and variance of the signal not to vary with the time, which typically does not hold for iEEG signals. The differencing has the effect of de-trending and stabilizing the data, and the first-order differencing is often sufficiently suppressing slow trends, although of course not guaranteeing perfect stationarity. For the details, please see Seth (2010).

3.1. Visual experiment

During the visual recognition experiment, four categories of pictures were demonstrated to the participants on a computer screen: scenes, small objects of daily life, faces, and finally a control category of edibles (fruits or vegetables). Images of the last category were used to control for a potential attention decrease—subjects were instructed to press a button every time they see a picture of this category. In the case of scenes, faces and objects no action was expected from the participants, button presses for these categories were considered an error.

Each of the categories of interest (scenes, objects, faces) included 100 different pictures. Each of these was shown in two different trials, giving the total of 200 trials per category, in pseudo-random order. The number of different pictures of edibles (fruits/vegetables) was 25, each being shown twice giving rise to 50 trials in total. Stimuli were presented for 300 ms, each followed by 800 ms of a black screen with a white cross in the center. For the analysis, the data were cut so that every trial had 200 ms of baseline (the black screen before the stimulus demonstration), 300 ms of picture demonstration, and 500ms of the following black screen. Trials with fruits/vegetables, as well as with button presses by mistake were excluded from the analysis. Also, trials with epileptogenic activity were excluded. A detailed description of stimuli, the task and the procedure used for removing the data with epileptogenic activity can be found in Vlcek et al. (2020).

3.2. ROI definitions

A previous study on the same dataset (Vlcek et al., 2020) has defined 7 clusters of channels that were of particular interest as they have shown a stronger reaction to the visual stimulation by scenes than by objects. These thus define effective regions of interest (ROIs), the connectivity among which we shall attempt to establish.

However, connectivity analysis on that subset of channels might be too restrictive because of the limited coverage of the visual pathways and consequently a relatively small amount of edges that could be tested. Therefore, we extended the dataset with all the channels that reacted to any type of stimuli, and lied close to the previously defined clusters. In particular, we constructed a 3D convex hull of every cluster; and extended every hull to 0.5–4.5 mm such that corresponding regions have more uniform volumes. Those increased convex hulls then defined, which channel belong to which ROI. Number of channels from implanted electrodes per ROI can be found in Table 1; their position within a brain template is shown in Figure 1. Note that we didn't distinguish the hemispheres, i.e., we consider every ROI symmetric, laying in both hemispheres.

TABLE 1

Table 1. Number of available channels per ROI for analyzed subjects.

FIGURE 1

Figure 1. Channels from implanted electrodes (aggregated across all subjects) with ROIs marked by different colors. Connections between ROIs, available for testing, are shown by light blue lines.

The centroids of these clusters were localized in the following structures:

1. OPA - the posterior angular and middle occipital gyrus (MNI [38, –76, 24], Occipital Place Area),

2. pLG - the posterior collateral sulcus at the junction with the lingual sulcus (MNI [25, –72, –8]),

3. PPA - the lingual and fusiform gyrus along the middle collateral sulcus (MNI [30, –45, –7]], Parahippocampal Place Area),

4. aCOS - the parahippocampal and fusiform gyrus along the anterior collateral sulcus (MNI [32, –26, –22]),

5. PCUN - the precuneus (MNI [15, –61, 27],

6. MPA - the superior part of the lingual gyrus and precuneus next to the retrosplenial region (MNI [16, –53, 12], Medial Place Area),

7. HIP - the anterior hippocampus (MNI [25, –12, –20]).

4. Hypothesis

For the purpose of our current analysis, we developed a hypothesis about the connectivity between seven ROIs (Figure 2), based on our previous iEEG study (Vlcek et al., 2020) and other electrophysiological and anatomical published results (Cavanna and Trimble, 2006; Byrne et al., 2007; Kravitz et al., 2011; Bastin et al., 2013; Baldassano et al., 2016; Boccia et al., 2017; Epstein and Baker, 2019; Henriksson et al., 2019). First, the visual cortex areas project to parieto-occipital areas such as the Occipital Place Area (OPA, ROI 1) around the transverse occipital sulcus (Kravitz et al., 2011). The documented responses as early as 60 ms after scene presentation (Henriksson et al., 2019) support its position early in the visual scene analysis network. As our data sample did not contain any channels located in primary visual cortex, we considered the OPA to be the entry point of the studied network of 7 regions.

FIGURE 2

Figure 2. Hypothesis about visual signal spreading via ventral/dorsal visual pathways.

Second, the OPA is functionally connected with the Parahippocampal Place Area (PPA, ROI 3) in the posterior temporal cortex (mainly its posterior part Baldassano et al., 2016). Several studies documented an early scene selectivity in the PPA at around 150–200 ms (Bastin et al., 2013; Vlcek et al., 2020). Based on the results of Vlcek et al. (2020), we distinguished another early scene selective area in the lingual gyrus (pLG, ROI 2). Because of its large selective responses to scenes but a lack of published information about its role in scene analysis, we located it in our hypothetical network in parallel to the PPA area. Therefore, both PPA (ROI 3) and pLG (ROI 2) are hypothesized to receive input from ROI 1, with unclear interaction between ROI 2 and 3.

Next, the resting state fMRI data suggest the PPA area connects with more anterior temporal regions including the anterior hippocampus (ROI 7) (Boccia et al., 2017). In our previous study (Vlcek et al., 2020), we documented late anterior hippocampal scene selectivity starting at around 300–350 ms and even earlier responses at about 200–250 ms around anterior collateral sulcus, which is positioned anatomically between the PPA and hippocampus. We added this region, located between ROI 3 and 7, as ROI 4.

Previously, we documented the scene selectivity also in the dorsal visual stream. In the posterior precuneus, the scene selectivity appeared around 200–250 ms (Vlcek et al., 2020), which is in agreement with its role in visual imagery and spatial attention (Cavanna and Trimble, 2006). Due to its location in the parietal cortex, we added it to our hypothetical network as ROI 5, parallel to ROIs 2 and 3.

Finally, the ROI 6 in our hypothesis corresponds to the Medial Place Area (MPA, also known as Retrosplenial Complex), a region with well documented scene selectivity also from fMRI studies (along with OPA and PPA, see Epstein and Baker (2019) for review). We previously documented scene selectivity in this region starting at around 300–350 ms, later than in the posterior precuneus, but similarly late to the hippocampal responses. As this area was documented to be part of a parieto-medial temporal pathway (Kravitz et al., 2011), translating visuo-spatial information between the dorsal and ventral streams (Byrne et al., 2007), we added it to our hypothetical network, potentially mediating the information flow between the ROI 5 in the dorsal stream and ROI 7 in the ventral stream.

5. Methods

5.1. Connectivity analysis: general approach

For the connectivity analysis purposes, we chose two different connectivity measures. The Phase Locking Value (PLV, Lachaux et al., 1999, see Section 5.2) is an example a undirected measure of synchrony between two time series, while the Directed Transfer Function (DTF, Kaminski and Blinowska, 1991, see Section 5.3) is one of most commonly used directed connectivity measures, i.e., it estimates directed interactions between channels.

As discussed earlier, implantation according to the clinical needs implies limitations on the collected data. For instance, connectivity between two regions can be naturally computed only for those patients, who have electrodes in both regions. Note that in principle one may attempt assessing dependence between signals from different regions in two different subjects and estimate thus some shared, population-level, component of the functional or even effective connectivity, and ultimately of the information flows, it is not entirely straightforward what interpretation would be suitable for such connectivity measure in view of inter-subject variability and noise. While at least at the level of correlations, inter-subject dependencies (even inter-species correlation) has been used e.g., in natural viewing data to identify homological areas across between monkey and human (Mantini et al., 2012), due to the complicated interpretability, we leave this as a matter of future research and restrict ourselves to pairs of channels observed within any single subject.

Figure 3 shows the number of implanted pairs of ROIs and corresponding number of channel pairs. Ideally, we would like to analyze the data with all 7 ROIs sampled in the same patient, but we don't have such subjects—at maximum, 3 ROIs were measured together. Also, as shown in Figure 3, not all links between ROIs were measured. Note that this is a typical situation, rather than an exception, in case of human iEEG studies. Thus, one is limited to analyzing only the observable connections between the regions by calculating connectivity per subject and then collected links per pair of ROIs to form the group-level network connectivity graph.

FIGURE 3

Figure 3. Information about the number of implanted patients (black numbers) and pairs of channels (blue numbers) for each pairs of ROIs. In the (A), ij-th element of the matrix corresponds to the number of patients, that have implanted both ROI i and j (black value) and available number of pairs of channels spanning ROIs i and j (blue value). Color-scale corresponds to the number of patients with implanted pair of ROIs, i.e., it corresponds to the values in black. In the (B), the graph of the hypotheses is combined with the graph of the connections that can be assessed from the data. The hypotheses are shown as blue dashed arrows. The implanted ROI pairs are connected with a black dotted edge, the width corresponds to the number of patients who have that pair of ROIs implanted. For example, for the connection from 3 to 4, the number of implanted subjects and pairs of channels can be found on the (A) (5 subjects, 20 pairs of channels); corresponding hypothesis and pairs of ROIs with implanted electrodes can be found on the (B).

The scheme of forming the group-level network connectivity graph from the subject-specific graphs is described in detail in Section 5.4. For a particular pair of channels (patient-specific), we compute the binary map of significance of the time- and frequency-resolved connectivity (see Figure 4). Significance testing for every time-frequency point is carried out using the multivariate Fourier surrogates (Prichard and Theiler, 1994), estimated from the baseline time series segment; the significance testing thus provides a binary time-frequency map of connectivity between the two channels. Note that the surrogates used are an extension of the phase-randomized Fourier-transform algorithm for generating surrogate data to multivariate time series. These surrogates are produced from the original multivariate data by applying a Fourier Transform, adding a random phase shift to each frequency (the same shift to all variables in a given surrogate realization), and applying inverse Fourier Transform to convert the randomized signals back to the time domain. As shown by its authors, such surrogate data sets mimic not only the auto- correlations of each of the variables in the original data set, but they mimic also the cross correlations (at each frequency, i.e., in principle the full cross-spectrum) between all pairs of the variables. Thus, such surrogates entail a realization of a stationary process with exactly the same (linear) properties as the original data. Surely, for shorter segments of data, the estimates of DTF and PLV naturally differ between the surrogates, providing thus an efficient empirical estimate of the expected spread of these statistics under the hypothesis of no true temporal dynamics (nonstationarity) in the coupling structure.

FIGURE 4

Figure 4. Diagram of the pipeline showing the subsequent steps starting from the time series measured from a single patient, through computation of dynamic connectivity, binary maps, heatmaps and test statistics. DTF is used as an example of connectivity measure.

Then, binary maps for all pairs of channels belonging to a given pair of ROIs are summed up (across all relevant channel pairs within a subject, and across all subjects) to form a ROI-pair specific connectivity density (“heatmap,” see Figure 4). While this heatmap gives a good qualitative overview of which connections are active at which time and frequency, carrying out multiple testing comparison across the many dimensions of the problem (time, frequency, channel pairs, subjects) might lead to severely decreasing the statistical power. To tackle this problem, as the connectivity can be expected to be relatively smooth over time and frequency (which is confirmed upon visual inspection), we chose to further proceed by conservatively evaluating the overall strength of the link between a pair of ROIs by computing an average heatmap value from the period after the stimulation.

The average heatmap thus captures the percentage of channel pairs, that showed significant connectivity at a given time and frequency bin. Notably, even in the case of no true connectivity, and in particular in the case of no difference from the baseline resting activity, one may expect some significant results. In theory, this percentage should correspond to the significance level selected for the channel-pair test. However, even though the standard approach of multivariate Fourier surrogates is used, the baseline data may violate its assumptions—particularly stationarity (and gaussianity), and thus having false positive rate increased above the nominal significance level even in the baseline segment. To further conservatively correct for this effect, we corrected the average heatmap value of the response period by subtracting the average heatmap value during the baseline. To provide a statistical inference regarding this final value of connectivity strength between ROIs, it is finally compared to a histograms of values of the same statistic obtained from multivariate surrogates of baseline data.

5.2. Connectivity measure: PLV

Phase Locking Value (PLV) was defined in Lachaux et al. (1999) and it is a time-dependent connectivity measure, that can be used to investigate task-induced (changes in) synchronization between time series. PLV is defined as

$PL{V}_{ij}(t)=\frac{1}{N}\left|{\displaystyle \sum _{n=1}^N{e}^{ı\theta (t,n)}}\right|\in [0,1],$

where N is a number of trials, θ(t, n) = ϕ_i(t, n)−ϕ_j(t, n), and ϕ_i(t, n) and ϕ_j(t, n) are the instantaneous phases for signals s_i and s_j in trial n at time t. In our study, the instantaneous phase of the signal was computed using a Hilbert transform. Higher values of PLV, with maximum of 1, mean higher synchronization between the time series (in terms of similar phase-lag between two signals across trials), while zero PLV corresponds to independent instantaneous phases. Before computing the PLV, the measured signal was filtered in several frequency bands, using a band-pass filter. Detailed pipeline is described in Section 5.4.

5.3. Connectivity measure: DTF

Directed Transfer Function (DTF) is a connectivity measure based on VAR(p) (Vector Autoregressive Process of the order p, see Lütkepohl (2005) for a general overview).

Let $\left\{X_t\in {ℝ}^n:t=1,\dots ,T\right\}$ denote a multivariate time series. The VAR(p) is defined by

$X_t={\displaystyle \sum _{k=1}^p{A}_k}{X}_{t-k}+{\epsilon }_t,$

where $X_t\in {ℝ}^n$, for t = 1, …, T, each A_k represents a coefficient matrix of dimension n×n, ε_t is a zero mean Gaussian random vector of size n (white noise), and p denotes the maximal lag length. This model can also be rewritten in the form:

${\epsilon }_t={\displaystyle \sum _{k=0}^p{\tilde{A}}_k}{X}_{t-k},$

where Ã₀ is an identity matrix and Ã_k = −A_k for k>0. This can be transformed into the frequency domain:

$X(f)={\tilde{A}}^{-1}(f)E(f)=H(f)E(f),$

where H(f) is a non-symmetric transfer matrix of the system and E(f) is a noise component in a frequency domain. Then, Directed Transfer Function (see Kaminski and Blinowska, 1991; Dauwels et al., 2010) can be defined in terms of H:

${\gamma }_{ij}^2(f)=\frac{|H_{ij}(f){|}^2}{{\displaystyle {\sum }_{j=1}^m|}{H}_{ij}(f){|}^2}\in [0,1],$

where the frequency-dependent normalization is chosen so that ${\gamma }_{ij}^2(f)$ quantifies the fraction of (total) inflow to channel i that arrives from channel j. To obtain one value per link j→i, we averaged the ${\gamma }_{ij}^2(f)$ values through the frequencies f = 1, …, 256. DTF was computed in a sliding window of 100 ms length with 50 ms shift. Similarly to PLV, before computing the DTF connectivity, the signal was band-pass filtered in the same frequency bands. Model order p for VAR(p) was put to 10 for all the patients, in order to decrease heterogeneity in the results.

5.4. Statistical analysis/detailed pipeline

The general diagram of the analysis pipeline is shown in Figure 4, detailed description of the pipeline is below. We denote:

• T_DTF = {[-200:-100]ms, [-150:-50]ms, … , [700:800]ms} is a moving window for DTF computation; baseline and reaction subsets for DTF are T_bsl = {[-200:-100]ms, [-150:-50]ms, … , [-50:50]ms} and T_reac = {[0:100]ms, … , [700:800]ms},

• T_PLV= [1,…,510] is the time index for PLV; baseline and reaction subsets for PLV are T_bsl = [1, …, 102] and T_reac = [103, …, 512],

• $\mathscr{H}=\{$[1:2]Hz, [2:4]Hz, [4:8]Hz, [8:16]Hz, [16:32]Hz, [32:64]Hz, [64:128]Hz} is a set of frequency bands,

• N^pat is the number of patients included in the analysis,

• p = 1, …, N^pat = 15 is a patient number,

• N(p) is a number of analyzed channels for patient p,

• i, j are channel numbers, i, j = 1, …, N(p),

• n(r₁, r₂) is a number of connections between regions r₁ and r₂,

• r(.)∈{1, …, 7} is a function of channel number, providing the number of the ROI to which the channel belongs,

• s = s(t, i)[p] is the measured signal, set of multivariate time series for every patient p; t = 1, …, 512 is time index, and i is the index of the channel, i∈{1, …, N(p)}.

• s_f = s_f(t, i)[p] - measured signal s(t, i)[p], band-pass filtered to frequency band $f\in {H}$. Butterworth filter of order 4 was used for that purpose.

For the measured data, we compute time-frequency resolved connectivity per patient between all channels F(t, f)[i, j, p] on the filtered signal s_f

$F(t,f)[i,j,p]=PL{V}_{ij}(t),t\in T_{PLV},f\in \mathscr{H},$

or

$F(t,f)[i,j,p]={\displaystyle \sum _{\text{f}=1}^{256}{\gamma }_{ij}^2}(\text{f}),t\in T_{DTF},f\in \mathscr{H},$

where p is a index of patient, i, j, i≠j are indices of channels i, j = 1, …, N(p) and f denotes the frequency band. Then, for signal of every subject s[p] we generate 100 multivariate Fourier surrogates from the baseline, i.e., the interval [-200:0] ms. Such surrogate data preserve the signal smoothness and correlation structure present in the (subject-, channel- and trial- specific) baseline. From the described surrogates, we computed the same time-frequency resolved connectivity F_sim(t, f, k)[i, j, p] as we did for the data, k = 1, …, 100 is a surrogate index.

Based on the computed connectivity for data and surrogates, we created the time-frequency binary map of significance per pair of channels for data and corresponding surrogates B(t, f)[i, j, p] and B_sim(t, f)[i, j, p]:

$\begin{array}{l}B(t,f)[i,j,p]=𝕀\{F(t,f)[i,j,p]>q(F_{sim}(t,f)[i,j,p],0.95)\},\\ B_{sim}(t,f,k)[i,j,p]=\\ 𝕀\{F_{sim}(t,f,k)[i,j,p]>q(F_{sim}(t,f)[i,j,p],0.95)\},\end{array}$

where 𝕀 is an indicator function, q(x, p) is a p-level quantile of sample x.

Then, for every pair of ROIs, we created a connectivity density (”heatmap”) for data and surrogates by averaging all the binary maps for pair of ROIs:

$\begin{array}{l}S(t,f)[r_1,r_2]=\frac{1}{n(r_1,r_2)}{\displaystyle \sum _{p=1}^{N^{pat}}{\displaystyle \sum _{i,j=1}^{N(p)}B}}(t,f)[i,j,p]\\ 𝕀\{r(i)=r_1,r(j)=r_2\},\\ S_{sim}(t,f,l)[r_1,r_2]=\frac{1}{n(r_1,r_2)}{\displaystyle \sum _{p=1}^{N^{pat}}{\displaystyle \sum _{i,j=1}^{N(p)}{B}_{sim}}}(t,f,{\xi }_{l,p})\\ [i,j,p]𝕀\{r(i)=r_1,r(j)=r_2\},\end{array}$

where 𝕀 is an indicator function, ξ_{l, p} is a random number indicating which of the 100 surrogate data for patient p should be used for the l-th simulated group heatmap, r(i) is number of ROI including channel i, n(r₁, r₂) is number of heatmaps for ROI pair r₁ and r₂, t∈T_DTF or t∈T_PLV, $f\in {H}$, and l = 1, …, 1000 is number of simulated heatmaps.

Further, we computed mean reaction value l(r₁, r₂) as average heatmap in reaction window minus average heatmap in the baseline:

$\begin{array}{l}(r_1,r_2)=\frac{1}{\left|\mathscr{H}\right|}{\displaystyle \sum _{f\in \mathscr{H}}(\frac{1}{\left|T_{reac}\right|}}{\displaystyle \sum _{t\in T_{reac}}S}(t,f)[r_1,r_2]\\ -\frac{1}{\left|T_{bsl}\right|}{\displaystyle \sum _{t\in T_{bsl}}S}(t,f)[r_1,r_2]),\\ l_{sim}(r_1,r_2,k)=\frac{1}{\left|\mathscr{H}\right|}{\displaystyle \sum _{f\in \mathscr{H}}(\frac{1}{\left|T_{reac}\right|}}{\displaystyle \sum _{t\in T_{reac}}{S}_{sim}}(t,f,k)[r_1,r_2]\\ -\frac{1}{\left|T_{bsl}\right|}{\displaystyle \sum _{t\in T_{bsl}}{S}_{sim}}(t,f,k)[r_1,r_2]),\end{array}$

where T_bsl and T_reac are time points/windows of baseline and reaction, respectively, r₁, r₂ = 1, …, 7 are ROI indices, and k = 1, …, 1000 is number of simulations. Then, p-value for every pair of ROIs was computed by comparing the statistics l with the distribution l_sim, obtained from surrogates as the percentile: l = q(l_sim, P(r₁, r₂)). In the Section 6 we discussed graphs, obtained from the matrix P by applying the Hochberg's procedure (Hochberg, 1988) for controlling the family-wise error rate (FWE) at the 0.05 error level.

5.5. Software

The data preparation and the analysis were performed in MATLAB 2020b (Mathworks, Inc.). DTF was computed using MATLAB implementation (Omidvarnia et al., 2011, Omidvarnia, 2020). Vector autoregression coefficients were computed using the ARFIT package (Neumaier and Schneider, 2001). PLV was computed using MATLAB implementation (Namburi, 2011), based on Lachaux et al. (1999). Pictures of brain were generated via BrainNet (Xia et al., 2013).

6. Results

We computed two general interaction graphs between the considered 7 ROIs: the undirected (based on PLV) and directed (DTF-based) connectivity network. Figure 5 shows the links between 7 ROIs with significant increase after the stimuli. The corresponding p-values can be found in Tables 2, 3.

FIGURE 5

Figure 5. Connectivity graphs between the 7 ROIs, obtained by DTF [(A,C), directed graphs, violet links] and PLV [(B,D), undirected graphs, green links]. Statistical significance was computed for the difference of the average heatmap value after and before the stimulus. Blue dashed lines in the (C,D) show the hypothesis of connections between ROI; black dotted lines show links presented at least in one subject; violet/green thin lines show significant links (uncorrected, p < 0.05); violet/green lines show significant links (FWE-corrected, p < 0.05).

TABLE 2

Table 2. P-values for testing the null hypothesis of no significant increase of the DTF during the reaction interval.

TABLE 3

Table 3. P-values for testing the null hypothesis of no significant increase of the PLV during the reaction interval.

Using PLV, we found a statistically significant increase in synchronization after the stimuli in all the assessed links between ROIs. In contrast, the causal inference created a much sparser network, with only one significant connection, when using the strict multiple comparison correction—The ROI 6 between the dorsal and ventral streams influenced ROI 3, in the ventral stream. More connections were revealed without multiple comparison correction. In particular, in addition to the link from ROI 6 to ROI 3, we found increased connectivity in bidirectional interaction between ROIs 3 and 7, and between ROIs 1 and 6. In addition, the analysis revealed increase in unidirectional connectivity from ROI 6 to 4 and from ROI 4 to 7.

To reiterate, depending on the implantation scheme of each patient, the connectivity of each ROI with other ROIs was based on different channels from different patients. Figure 3 summarizes the number of patients and connections between channels in one or another ROI. To know specifically, from which brain regions the specific ROI connectivity emerged, we investigated in detail the localization of involved channels.

The connectivity from ROI 6 to 3 was based on 9 connections between 10 channels. The six ROI 6 channels were localized mostly on the border between retrosplenial cortex (RSC) and either precuneus or the superior part of the lingual gyrus (sLG), while the four ROI 3 channels were all localized around the border between the anterior and posterior part of the parahippocampal place area (PPA), with the MNI y coordinate between –53 and –31 (–42 being the anterior/posterior boundary according to Baldassano et al., 2016). The interaction between ROIs 3 and 7 is based on data from 17 channels. Of them, the 12 channels in ROI 3 were primarily localized in the fusiform gyrus (FG) or lingual gyrus (LG), with MNI y coordinate mostly (10 of 12) posterior to –42, which again corresponds to posterior PPA. On the other hand, the 5 channels in ROI 7 were localized in the anterior hippocampus (HIP), entorhinal cortex (EC) and temporopolar gyrus, with MNI z coordinate below –19 (–9 being the anterior/posterior hippocampus boundary according to Morgan et al. (2011). The interaction between ROIs 1 and 6 is based on data from one subject, with two channels in ROI 1, both localized between middle occipital gyrus and posterior angular gyrus (in transverse occipital sulcus), and four channels in ROI 6, all localized just between superior lingual gyrus (sLG) and retrosplenial cortex (RSC). The connectivity from ROI 6 to 4 is based on data from two subjects, with four channels in ROI 6, all on the border of RSC and precuneus or LG, and four channels in ROI 4, all in posterior parahippocampal gyrus (PHG), with MNI y coordinate between -27 and –33 (the anterior/posterior border being at –20 according to Pruessner et al. (2002). The connectivity from ROI 4 to 7 is based on data from 27 channels. Of them, 12 channels in ROI 4 were mainly localized in posterior part of PHG and adjacent FG, with MNI y coordinate posterior to –24, and 15 channels in ROI 7 were located primarily in the anterior HIP, EC and amygdala, with MNI z coordinate below –18.

Figures 6, 7 present the heatmaps that were used for computation of graphs in Figure 5. We divided the heatmaps to two classes: connectivity from ROI to itself (Figure 6), and connectivity between the different ROIs (Figure 7). Recall that the heatmap of connectivity from ROI i to itself was estimated from the collection of connections between different channels, both of them laying in ROI i. Heatmap ROI i to ROI j, i≠j was computed from pairs of channels, for which the source channel belonged to ROI i, and the target channel was in ROI j.

FIGURE 6

Figure 6. Time-frequency heatmaps of number of significant links per pair of ROI that are computed for the channels with significant reaction to the stimuli.

FIGURE 7

Figure 7. Time-frequency heatmaps of number of significant links per pair of ROI that are computed for the Channels with significant reaction to the stimuli.

Every heatmap shows the dynamical connectivity for several frequency bands (y-axis), with the time on the x-axis marking the start of sliding window for DTF, and the relevant time instant for PLV. Color of every entry corresponds to the percentage of significant (uncorrected, channel-level) links in particular time-frequency point between investigated pair of ROIs. Higher values of the heatmap (color closer to yellow and white) mark a high proportion of links showing increased connectivity and suggest at which time and frequency an increased communication emerges between the respective regions.

A marked increase in synchronization was often observed at the frequency bands 2 to 4 Hz and 4 to 8 Hz. The particular time interval of this increase is difficult to define and seems to depend on the particular pair of ROIs. DTF increase is not as pronounced as in PLV, but generally, it also showed an increase at 2–4 Hz frequency. DTF to/from ROI 3 and 6 seemed to have an increase at high frequencies, unlike the PLV, which showed a main increase at the low frequencies. Also, both of the connectivity measures showed a stronger increase inside the individual ROIs compared to the connectivity between ROIs. The presented heatmaps, however, are useful only for exploration purposes and should not be over-interpreted. We reiterate that the presented heatmaps element-wise were not associated with the statistical p-values (they correspond to proportion of uncorrected channel-level connection increases), and (corrected) statistical testing was carried out later at a summary level, aggregating results across frequencies and time bins. The heatmaps thus can not be considered as a result of inductive statistics (in the sense of testing a hypothesis), but only rather as descriptive statistics useful for exploration or interpretation of the detected results of higher-level hypothesis testing. Only in this sense we believe that the time-frequency heatmaps can be useful—for neuroscientific interpretation of those connections, that reached global statistical significance, in guiding the interpretation in terms of timing and dominant frequencies behind the effect; or as pure descriptive exploration. For sure, such interpretation is limited—in a similar sense as the interpretation of statistical parametric maps presented in neuroimaging studies. We also analyzed a connectivity decrease after the stimulation but didn't find significant results.

7. Discussion

As stated in the introduction, the inference of brain graph analysis from iEEG data faces an extensive set of challenges. We have outlined these in the introduction, proceeding from those that are shared across a large set of inference situations, to those that are relatively specific to the situation of iEEG—namely the partial, sparse and heterogeneous sampling of the relevant brain network in any single patient. Based on the discussion of the possible ways to tackle these challenges, we have made a range of methodological decisions in order to build up a pipeline allowing robust statistical inference of the alterations of brain connectivity network during processing visual stimuli. The application of the pipeline was demonstrated on an example dataset, giving rise to an group-level network in terms of both functional and effective connectivity.

Note that the functional and effective connectivity correspond to different yet related concepts—the true underlying interactions (in principle captured by the effective connectivity) can be considered as the cause of the synchronization of the brain activity. However, the true interactions that are to be estimated by effective connectivity methods might not necessarily happen exactly between the pair of regions for which synchrony is observed. Indeed, due to common drivers or mediators, two regions can be synchronized in their activity, while not interacting directly. In our particular study, we have observed significant increased of functional connectivity (captured by PLV) in all analyzed connections. However, the application of effective connectivity analysis (DTF) was confirmed the increase a sparser set of directed connections, providing more insight into the stimulus-driven information flows.

Note that for the sake of conserving enough channels for subsequent analysis, we have not excluded the connectivity estimates between pairs of bipolar channels that shared one contact. Indeed, such channels might have biased estimates of PLV or DTF due to sharing one of the two signals in their construction. However, as we are not focusing on estimation of the connectivity strength per se, but only on its change in time, this issue is not so paramount in current analysis. To check the robustness of this analytical choice, we have rerun the analysis while removing the affected channels pairs. Note that this omission affects the channel-pair count quadratically, and leaves less than half of channels pairs in most region pairs for the analysis. Nevertheless, the obtained heatmaps characterizing the connectivity changes over time were qualitatively very similar to the original analysis (controlled vs. original heatmap correlation ranging between 0.6 and 1, on average about 0.8).

Also, while we have made a relatively common choice of connectivity measures from each family, the results could surely be (potentially substantially) different for other choices of connectivity measures. In particular, the functional connectivity measure used (i.e., PLV), is in two important aspects different from the classical counterpart of DTF in the functional connectivity family, that is, from the coherence. In particular, the PLV works with phase only (unlike the DTF and coherence, which depend on both the phase and amplitude of the signal). Secondly, the PLV measures the synchronization across trials, rather than synchronization across time (as would mean phase coherence do). We have had two-fold motivation for this choice of functional connectivity selection: firstly, PLV is an widespread and commonly used method. Secondly, the heterogeneity of connectivity methods used helps to illustrate the wide applicability of the pipeline. Of course, the framework is in principle applicable along with a range of other connectivity methods, including the natural alternative to DTF, i.e., the partial directed coherence (Sameshima and Baccalá, 1999; Baccalá et al., 2001; Baccalá and Sameshima, 2021). Yet, one may surely argue for the application of many alternative methods from the plethora of the existing ones.

Although the specific application was not the main rationale of this study, the utilization of the developed pipeline provided interesting results. In the studied case of network changes in reaction to the visual stimuli, the developed methodology allowed us to robustly detect several interesting phenomena.

First of all, we confirmed the involvement of the parieto-medial temporal pathway in scene perception, translating visuospatial information between dorsal and ventral visual streams during visual scene analysis. This pathway connects the parietal lobe with locations in the medial temporal lobe (ROI 3 and 4), via MPA (ROI 6, also called retrosplenial complex). The existence of this pathway was suggested based on rat studies showing that associative parietal cortex lesions affect the firing of place cells in the medial temporal lobe (Save et al., 2005). The pathway is claimed to include also the retrosplenial region (Byrne et al., 2007). Based on anatomical data, it connects the angular gyrus (homological to caudal inferior parietal lobule in animals) in the parietal lobe to the medial temporal lobe either directly or indirectly, with the indirect connection leading via the retrosplenial cortex and posterior cingulate cortex (Kravitz et al., 2011). Both direct and indirection connections were also demonstrated in humans using resting-state fMRI data (Boccia et al., 2017). Our data show that this indirect pathway (from ROI 6 to ROIs 3 and 4) is active also during static visual scenes processing in humans. In agreement with another connectivity study (Baldassano et al., 2016), we show the MPA connection with the PPA (ROI 3), around the border between its anterior and posterior parts.

Second, we confirmed the anterior hippocampal region (ROI 7) connectivity with more posterior areas in the medial temporal lobe: the parahippocampal place area (PPA, ROI 3) and region around the anterior collateral sulcus (ROI 4). The hippocampus (HIP) is a medial temporal structure associated primarily with declarative and spatial memory, but its function and connectivity vary along its long anterior-posterior axis (Strange et al., 2014). While the function of the posterior portion (pHIP) seems to be related to spatial memory, the anterior part (aHIP) has been associated with emotional processing, novelty detection, and semantic memory. Accordingly, there are connectivity differences, with the aHIP being linked, among others, to the amygdala, ventral tegmental area and hypothalamus, while pHIP is coupled to posterior parietal areas and retrosplenial cortex (Poppenk et al., 2013). However, there seems to be some overlap, as the aHIP shows connection also to more posterior areas, like anterior PPA and MPA (Boccia et al., 2017), in one study even larger than pHIP (Baldassano et al., 2016). Extending these studies, our current results document a connection of aHIP (ROI 7) with even posterior part of PPA (ROI 3). The hippocampal connections to the parahippocampal gyrus (PHG) also differ between its anterior to posterior part. While aHIP shows preferential connectivity with the perirhinal cortex (anterior part of PHG), pHIP is preferentially connected to the parahippocampal cortex (posterior part of PHG) (Libby et al., 2012). However, this preferential connectivity pattern does not exclude the functional coupling of other areas. Our data show information flow from posterior PHG (ROI 4) to aHIP (ROI 7), but our collection of channels did not contain any channel located in pHIP. Therefore, for both posterior PPA and parahippocampal cortex, we could not exclude their indirect association with aHIP via pHIP.

Third, we found the reciprocal information flow between ROI 1 in OPA and ROI 6 in MPA. At the same time, we hypothesized to find connectivity between ROI 1 and PPA (ROI 3), but we found no patients with both ROIs implanted in our dataset. Related to that, a recent resting-state fMRI connectivity study (Baldassano et al., 2016) distinguished two scene processing networks. The anterior one connects angular gyrus with retrosplenial cortex, anterior portion of PPA and anterior HIP. In contrast, the posterior one consists of OPA and posterior portion of PPA. In agreement with this study, the two channels in our ROI 1, showing connectivity with ROI 6 in our data, seem to be a part of the anterior scene processing network, which is related more to memory and navigation. In fact, these two channels in ROI 1 were localized still in medial occipital gyrus, but deep at the anterior end of transverse occipital sulcus, thus on the posterior margin of angular gyrus. Therefore, our data support the separation of anterior and posterior scene processing networks (Baldassano et al., 2016). Finally, we have obtained inconclusive evidence concerning the link between ROI 1 and ROI 2: although descriptively (Figure 7) there seems to be a clear increase in both PLV and DTF following the stimulus at low frequencies (2–4–8 Hz), only the PLV coupling was strong enough to reach significance, see Figure 5. This is likely to the very low number of pairs of channels spanning these two ROIs in the current dataset.

Albeit the methodology was tailored to the nature of the analyzed data, the results achieved are still clearly limited by some aspects of the data. For instance, the low spatial coverage, in particular the lack of implantation of posterior occipital cortex, prevented us from tracking the early development of the visual signal processing. Low number of patients involved not only affected the statistical power of the inference, but indirectly also the amount of pairs of regions that we were able to test for interactions, and together with the relative sparseness of the implantation it generally decreased the spatial resolution of the estimated network. Indeed, the regions were arguably large enough so that they would be functionally heterogeneous. On the other side, should one decide to work with smaller regions, there might not be sufficient number of subjects with implantation spanning any given pair of regions, leaving thus each connection inference specific to a single subject.

Needless to say, the results in principle might be affected by any residual differences in the cognitive processing, and, more importantly, brain dynamics in the epilepsy patients compared to healthy subjects. Although the volunteering patients all had normal vision, and any segments with detected interictal epileptiform activity (or seizures) were removed from the analyzed data, one can never exclude the possibility of some systematic bias in neuronal dynamics and functional neuroanatomy due to the chronic brain disease. Indeed, the selection of the regions here was based on a data driven approach used in a previously published study; an alternative approach used in the field is to apply some standard anatomical atlas.

Data availability statement

Pipeline scripts (in MATLAB) and the estimated connectivity matrices are available upon request from the corresponding author. Requests to access these datasets should be directed to JHl, aGxpbmthJiN4MDAwNDA7Y3MuY2FzLmN6.

Ethics statement

This study was approved by the Ethics Committee of Motol University Hospital and all patients gave their informed consent to participate. Written informed consent to participate in this study was provided by the participant or participants' legal guardian/next of kin.

Author contributions

AP, KV, and JHl conceptualized, validated the study, and wrote the original draft. AP, AK, JHa, KV, and JHl curated the data. AP and JHl performed the formal analysis. PM, KV, and JHl provided the funding acquisition. AP, JHa, PS, KV, and JHl designed the methodology. AP, PM, KV, and JHl did the administration of the project. PM, AK, KV, and JHl provided the resources. AP prepared the software and visualization. JHl supervised the project. All authors contributed to the investigation, writing, editing, and review of the manuscript.

Funding

This study was supported by the Czech Science Foundation project No. 19-11753S.

Acknowledgments

We thank all patients for their participation, and Nad'a Bednárová, Helena Buchtová, Lucie Paterová, Iveta Fajnerová, Tereza Nekovářová, and Jana Kalinová for their help with collecting the iEEG data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Aguirre, G. K., and D'Esposito, M. (1999). Topographical disorientation: a synthesis and taxonomy. Brain 122, 1613–1628. doi: 10.1093/brain/122.9.1613

PubMed Abstract | CrossRef Full Text | Google Scholar

Baccalá, L., Sameshima, K., Baccala, L. A., and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biol. Cybern. 84, 463–474. doi: 10.1007/PL00007990

PubMed Abstract | CrossRef Full Text | Google Scholar

Baccalá, L. A., and Sameshima, K. (2021). Partial directed coherence: twenty years on some history and an appraisal. Biol Cybern. 115, 195–204. doi: 10.1007/s00422-021-00880-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Baldassano, C., Esteva, A., Beck, D., and Fei-Fei, L. (2016). Two distinct scene processing networks connecting vision and memory. J. Vis. 3, 1–16. doi: 10.1523/ENEURO.0178-16.2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnett, L., and Seth, A. K. (2014). The mvgc multivariate Granger causality toolbox: a new approach to Granger-causal inference. J. Neurosci. Methods 223, 50–68. doi: 10.1016/j.jneumeth.2013.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastin, J., Deman, P., David, O., Gueguen, M., Benis, D., Minotti, L., et al. (2017). Direct recordings from human anterior insula reveal its leading role within the error-monitoring network. Cereb. Cortex 27, 1545–1557. doi: 10.1093/cercor/bhv352

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastin, J., Vidal, J. R., Bouvier, S., Perrone-Bertolotti, M., Bénis, D., Kahane, P., et al. (2013). Temporal components in the parahippocampal place area revealed by human intracerebral recordings. J. Neurosci. 33, 10123–10131. doi: 10.1523/JNEUROSCI.4646-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Biswal, B. B., Mennes, M., Zuo, X.-N., Gohel, S., Kelly, C., Smith, S. M., et al. (2010). Toward discovery science of human brain function. Proc. Natl. Acad. Sci. U.S.A. 107, 4734–4739. doi: 10.1073/pnas.0911855107

PubMed Abstract | CrossRef Full Text | Google Scholar

Blinowska, K. J., and Malinowski, M. (1991). Non-linear and linear forecasting of the EEG time series. Biol. Cybern. 66, 159–165. doi: 10.1007/BF00243291

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccia, M., Sulpizio, V., Nemmi, F., Guariglia, C., and Galati, G. (2017). Direct and indirect parieto-medial temporal pathways for spatial navigation in humans: evidence from resting-state functional connectivity. Brain Struct. Funct. 222, 1945–1957. doi: 10.1007/s00429-016-1318-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, J. F., Sharan, A. D., Sperling, M. R., Ramayya, A. G., Evans, J. J., Healey, M. K., et al. (2014). Theta and high-frequency activity mark spontaneous recall of episodic memories. J. Neurosci. 34, 11355–11365. doi: 10.1523/JNEUROSCI.2654-13.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, J. F., Zaghlou, K. A., Jacobs, J., Williams, R. B., Sperling, M. R., Sharan, A. D., et al. (2013). Synchronous and asynchronous theta and gamma activity during episodic memory formation. J. Neurosci. 33, 292–304. doi: 10.1523/JNEUROSCI.2057-12.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Byrne, P., Becker, S., and Burgess, N. (2007). Remembering the past and imagining the future: a neural model of spatial memory and imagery. Psychol. Rev. 114, 340–375. doi: 10.1037/0033-295X.114.2.340

PubMed Abstract | CrossRef Full Text | Google Scholar

Cavanna, A. E., and Trimble, M. R. (2006). The precuneus: a review of its functional anatomy and behavioural correlates. Brain 129, 564–583. doi: 10.1093/brain/awl004

PubMed Abstract | CrossRef Full Text | Google Scholar

Cloutman, L. L. (2013). Interaction between dorsal and ventral processing streams: where, when and how? Brain Lang. 127, 251–263. doi: 10.1016/j.bandl.2012.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Dauwels, J., Vialatte, F., Musha, T., and Cichocki, A. (2010). A comparative study of synchrony measures for the early diagnosis of Alzheimer's disease based on EEG. Neuroimage 49, 668–693. doi: 10.1016/j.neuroimage.2009.06.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, R., and Kanwisher, N. (1998). A cortical representation of the local visual environment. Nature 392, 598–601. doi: 10.1038/33402

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, R. A., and Baker, C. I. (2019). Scene perception in the human brain. Ann. Rev. 5, 373–397. doi: 10.1146/annurev-vision-091718-014809

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J. (1994). Functional and effective connectivity in neuroimaging: a synthesis. Hum. Brain Mapp. 2, 56–78. doi: 10.1002/hbm.460020107

CrossRef Full Text | Google Scholar

Friston, K. J., Frith, C. D., Liddle, P. F., and Frackowiak, R. S. J. (1993). Functional connectivity-the principal-component analysis of large (PET) data sets. J. Cereb. Blood Flow Metab. 13, 5–14. doi: 10.1038/jcbfm.1993.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodale, M. A., and Milner, A. D. (1992). Separate visual pathways for perception and action. Trends Neurosci. 15, 20–25. doi: 10.1016/0166-2236(92)90344-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartman, D., Hlinka, J., Palus, M., Mantini, D., and Corbetta, M. (2011). The role of nonlinearity in computing graph-theoretical properties of resting-state functional magnetic resonance imaging brain networks. Chaos 21, 013119. doi: 10.1063/1.3553181

PubMed Abstract | CrossRef Full Text | Google Scholar

Haufe, S., Nikulin, V. V., and Nolte, G. (2012). Alleviating the influence of weak data asymmetries on Granger-causal analyses. Lecture Notes Compute. Sci. 7191 LNCS, 25–33. doi: 10.1007/978-3-642-28551-6_4

CrossRef Full Text | Google Scholar

Henriksson, L., Mur, M., and Kriegeskorte, N. (2019). Rapid invariant encoding of scene layout in human OPA. Neuron 103, 161.e3–171.e3. doi: 10.1016/j.neuron.2019.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Hlinka, J., Hartman, D., Vejmelka, M., Runge, J., Marwan, N., Kurths, J., et al. (2013). Reliability of inference of directed climate networks using conditional mutual information. Entropy 15, 2023–2045. doi: 10.3390/e15062023

CrossRef Full Text | Google Scholar

Hlinka, J., Paluš, M., Vejmelka, M., Mantini, D., and Corbetta, M. (2011). Functional connectivity in resting-state fMRI: Is linear correlation sufficient? Neuroimage 54, 2218–2225. doi: 10.1016/j.neuroimage.2010.08.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika 75, 800–802. doi: 10.1093/biomet/75.4.800

CrossRef Full Text | Google Scholar

Kalman, R. (1960). A new approach to linear filtering and prediction problems. J. Basic Eng. 82, 35–44. doi: 10.1115/1.3662552

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, M., Ding, M., Truccolo, W. A., and Bressler, S. L. (2001). Evaluating causal relations in neural systems: granger causality, directed transfer function and statistical assessment of significance. Biol. Cybern. 85, 145–157. doi: 10.1007/s004220000235

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, M. J., and Blinowska, K. J. (1991). A new method of the description of the information flow in the brain structures. Biol. Cybern. 65, 203–210. doi: 10.1007/BF00198091

PubMed Abstract | CrossRef Full Text | Google Scholar

Kořenek, J., and Hlinka, J. (2020). Causal network discovery by iterative conditioning: comparison of algorithms. Chaos 30, 013117. doi: 10.1063/1.5115267

PubMed Abstract | CrossRef Full Text | Google Scholar

Korenek, J., and Hlinka, J. (2021). Causality in reversed time series: reversed or conserved? Entropy 23, 1–22. doi: 10.3390/e23081067

PubMed Abstract | CrossRef Full Text | Google Scholar

Kravitz, D. J., Saleem, K. S., Baker, C. I., and Mishkin, M. (2011). A new neural framework for visuospatial processing. Nat. Rev. Neurosci. 12, 217–230. doi: 10.1038/nrn3008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristensen, S., Garcea, F. E., Mahon, B. Z., and Almeida, J. (2016). Temporal frequency tuning reveals interactions between the dorsal and ventral visual streams. J. Cogn. Neurosci. 28, 1295–1302. doi: 10.1162/jocn_a_00969

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachaux, J.-,p., Rodriguez, E., Martinerie, J., and Varela, F. J. (1999). Measuring phase synchrony in brain signals. Hum. Brain Mapp. 208, 194–208. doi: 10.1002/(SICI)1097-0193(1999)8:4<194::AID-HBM4>3.0.CO;2-C

PubMed Abstract | CrossRef Full Text | Google Scholar

Libby, L. A., Ekstrom, A. D., Ragland, J. D., and Ranganath, C. (2012). Differential connectivity of perirhinal and parahippocampal cortices within human hippocampal subregions revealed by high-resolution functional imaging. J. Neurosci. 32, 6550–6560. doi: 10.1523/JNEUROSCI.3711-11.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. Berlin; Heidelberg: Springer-Verlag Berlin Heidelberg.

Google Scholar

Mantini, D., Hasson, U., Betti, V., Perrucci, M. G., Romani, G. L., Corbetta, M., et al. (2012). Interspecies activity correlations reveal functional correspondence between monkey and human brain areas. Nat. Methods 9, 277-U85. doi: 10.1038/nmeth.1868

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgan, L. K., MacEvoy, S. P., Aguirre, G. K., and Epstein, R. A. (2011). Distances between real-world locations are represented in the human hippocampus. J. Neurosci. 31, 1238–1245. doi: 10.1523/JNEUROSCI.4667-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Namburi, P. (2011). Phase Locking Value. Available online at: https://se.mathworks.com/matlabcentral/fileexchange/31600-phase-locking-value.

Google Scholar

Neumaier, A., and Schneider, T. (2001). Estimation of parameters and eigenmodes of multivariate autoregressive models. ACM Trans. Math. Softw. 27, 27–57. doi: 10.1145/382043.382304

CrossRef Full Text | Google Scholar

Omidvarnia, A. (2020). Time-Varying EEG Connectivity: A Time-Frequency Approach. Available online at: https://www.mathworks.com/matlabcentral/fileexchange/33721-time-varying-eeg-connectivity-a-time-frequency-approach.

Google Scholar

Omidvarnia, A., Mesbah, M., O'Toole, J., Colditz, P., and Boashash, B. (2011). “Analysis of the time-varying cortical neural connectivity in the newborn EEG: a time-frequency approach,” in 2011 7th International Workshop on Systems, Signal Processing and their Applications (WOSSPA) (Tipaza), 179–182.

Google Scholar

Poppenk, J., Evensmoen, H. R., Moscovitch, M., and Nadel, L. (2013). Long-axis specialization of the human hippocampus. Trends Cogn. Sci. 17, 230–240. doi: 10.1016/j.tics.2013.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Prichard, D., and Theiler, J. (1994). Generating surrogate data for time series with several simultaneously measured variables. Phys. Rev. Lett. 73, 951–954. doi: 10.1103/PhysRevLett.73.951

PubMed Abstract | CrossRef Full Text | Google Scholar

Pruessner, J. C., Köhler, S., Crane, J., Pruessner, M., Lord, C., Byrne, A., et al. (2002). Volumetry of temporopolar, perirhinal, entorhinal and parahippocampal cortex from high-resolution MR images: considering the variability of the collateral sulcus. Cereb. Cortex 12, 1342–1353. doi: 10.1093/cercor/12.12.1342

PubMed Abstract | CrossRef Full Text | Google Scholar

Sameshima, K., and Baccalá, L. A. (1999). Using partial directed coherence to describe neuronal ensemble interactions. J. Neurosci. Methods 94, 93–103. doi: 10.1016/S0165-0270(99)00128-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Save, E., Paz-Villagran, V., Alexinsky, T., and Poucet, B. (2005). Functional interaction between the associative parietal cortex and hippocampal place cell firing in the rat. Eur. J. Neurosci. 21, 522–530. doi: 10.1111/j.1460-9568.2005.03882.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Seth, A. K. (2010). A matlab toolbox for granger causal connectivity analysis. J. Neurosci. Methods 186, 262–273. doi: 10.1016/j.jneumeth.2009.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Strange, B. A., Witter, M. P., Lein, E. S., and Moser, E. I. (2014). Functional organization of the hippocampal longitudinal axis. Nat. Rev. Neurosci. 15, 655–669. doi: 10.1038/nrn3785

PubMed Abstract | CrossRef Full Text | Google Scholar

Ungerleider, L., and Mishkin, M. (1982). “Two visual streams,” in Analysis of Visual Behavior, eds D. Ingle, M. Goodale, and R. Mansfield (Cambridge, MA: MIT Press), 549–586.

Google Scholar

Vlcek, K., Fajnerova, I., Nekovarova, T., Hejtmanek, L., Janca, R., Jezdik, P., et al. (2020). Mapping the scene and object processing networks by intracranial EEG. Front. Hum. Neurosci. 14, 561399. doi: 10.3389/fnhum.2020.561399

PubMed Abstract | CrossRef Full Text | Google Scholar

Winkler, I., Panknin, D., Bartz, D., Muller, K. R., and Haufe, S. (2016). Validity of time reversal for testing granger causality. IEEE Trans. Signal Process. 64, 2746–2760. doi: 10.1109/TSP.2016.2531628

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, M., Wang, J., and He, Y. (2013). BrainNet viewer: a network visualization tool for human brain connectomics. PLoS ONE 8, e68910. doi: 10.1371/journal.pone.0068910

PubMed Abstract | CrossRef Full Text | Google Scholar