Topological Analysis of Seizure-Induced Changes in Brain Hierarchy Through Effective Connectivity

Abstract

Traditional Topological Data Analysis (TDA) methods, such as Persistent Homology (PH), rely on distance measures (e.g., cross-correlation, partial correlation, coherence, and partial coherence) that are symmetric by definition. While useful for studying topological patterns in functional brain connectivity, the main limitation of these methods is their inability to capture the directional dynamics - which is crucial for understanding effective brain connectivity. We propose the Causality-Based Topological Ranking (CBTR) method, which integrates Causal Inference (CI) to assess effective brain connectivity with Hodge Decomposition (HD) to rank brain regions based on their mutual influence. Our simulations confirm that the CBTR method accurately and consistently identifies hierarchical structures in multivariate time series data. Moreover, this method effectively identifies brain regions showing the most significant interaction changes with other regions during seizures using electroencephalogram (EEG) data. These results provide novel insights into the brain’s hierarchical organization and illuminate the impact of seizures on its dynamics.

1 Introduction

Over the past two decades, Topological Data Analysis (TDA) has become a valuable tool in various fields, providing new insights into complex data structures [1, 5, 6, 8, 9]. In biology, TDA has helped clarify genetic and evolutionary processes [19], while in materials science, it has improved predictions of properties in crystalline compounds [16]. In finance, TDA has offered a novel perspective on the dynamics of financial markets during crises [14]. In neuroscience, the use of Persistent Homology (PH), a key technique within TDA, has revealed new aspects of the topological organization of brain dependence networks. By analyzing multivariate brain signals, researchers have uncovered complex patterns and structures that deepen our understanding of brain connectivity and functionality [2, 10, 17].

However, traditional uses of PH often involve Vietoris-Rips filtrations from undirected brain dependence networks. A significant limitation of this approach is its inability to capture the directionality of connectivity between nodes in a brain network [11]. It is well-established that node dependencies within a brain network exhibit asymmetry, where the influence exerted by node $A$ on node $B$ can differ significantly from the influence exerted by node $B$ on node $A$. This asymmetry is crucial for understanding effective brain connectivity, which involves directional information flow, from one region into another, that often varies markedly in its reverse path, thus delineating specific neural pathways [12, 13, 20, 23].

In the study of epilepsy, understanding changes in brain connectivity during seizures is crucial. Research indicates that epileptic seizures is associated with significant alterations in interactions between brain regions in a network. These changes vary depending on whether the seizure is generalized or originates from a specific focal region and spreads across the brain’s hemispheres [4, 7]. Traditional TDA methods, which primarily assess functional connectivity, often fall short in accurately characterizing these dynamic changes in effective brain connectivity.

The concept of causal inference refers to learning a causal model from observational data. In cognitive neuroscience, this process is fundamental for understanding how different brain regions interact, and how these interactions may be altered due to brain disorders  [25]. Experimental interventions (e.g., optogenetics) are considered the gold standard in causal discovery [18, 25], as they provide a way to experimentally manipulate activity of a neuron (or subpopulation of neurons) and investigate how these manipulations change the interaction of a neuron with another neuron (or another subpopulation of neurons). Through these experimental manipulations, one can gain definitive evidence of causal physiological relationships. While such interventions and manipulations are conducted in non-human animal studies, these are ethically and practically challenging when studying brain connectivity in humans. Consequently, we will only focus on identifying causal associations and quantifying their strength from observational data.

To gain better understanding of the dynamic changes in brain connectivity during epileptic seizures, we have developed the Causality-Based Topological Ranking (CBTR) approach. This novel method integrates Causal Inference (CI) and Hodge Decomposition (HD), allowing for the identification of causal pathways and the hierarchical ranking of brain regions based on their mutual influence over time. CBTR assesses dependencies conditioned on parent nodes, thereby avoiding confounding effects among observed variables. Our analysis of EEG data from subjects who experienced seizures helped us identify the brain regions most significantly impacted by seizures for each subject.

2 Methods

The goal is to develop a method for ranking brain regions (derived from the network of electroencehphalogram data) according to effective connectivity using a two-step process. First, causal inference techniques are employed to map the hierarchical structure of brain connectivity, with a focus on the effects of seizures. Next, Hodge Decomposition is applied to net flow of information between brain regions, quantifying their mutual influence.

2.1 Causal Inference of the Brain’s Hierarchical Structure

Understanding the dynamics of the hierarchical structure of the brain’s effective connectivity is crucial for comprehending the impact of seizures. In an EEG network, a channel is considered to be at the top of the hierarchy if it tends to exert greater influence on other nodes, whereas a brain channel is at the bottom if it tends to be influenced more by other nodes. This hierarchical structure is determined by analyzing the causal relationships and time-delayed dependencies among the channels. This is illustrated in Figure 1. Assessing changes in brain hierarchy allows us to quantify the effects of seizures on each specific EEG channel, providing deeper insights into how seizures alter brain connectivity.

Figure 1: Potential hierarchy of the brain in normal (Left) and seizure (Right) conditions. In the normal state, the hierarchy is O1 (left occipital), O2 (right occipital) followed by T3 (left temporal), Fp2 (right frontal pole), and then T4 (right temporal), Fp1 (left frontal pole). During a seizure, the order changes to T3, O1, Fp1 followed by T4, O2, Fp2.

To estimate effective brain connectivity, we employ the PCMCI (Peter and Clark Momentary Conditional Independence) method, a robust causal inference technique well-suited for high-dimensional time series data [21, 22]. Following the PCMCI method, we conducted a two-stage process. First, we used the PC algorithm to identify potential causal parents for each variable, reducing dimensionality and producing a preliminary causal graph. Then, we applied the Momentary Conditional Independence (MCI) test to rigorously assess these dependencies’ significance, considering both contemporaneous and time-lagged relationships while controlling for false positives. This combination of dimension reduction and rigorous testing allowed PCMCI to effectively estimate the connectivity between observed time series components, providing a detailed map of directional dependencies in the brain.

The outputs of the PCMCI method include a causal graph that highlights the most significant causal links between brain regions at various lags, with each link annotated by its corresponding statistical significance and strength. This causal graph forms the foundation for our subsequent analyses. We estimate effective connectivity as a weighted directed network, where the weight of each link ($\mathbf{W}_{p,q}$) between regions $p$ and $q$ is derived from the significance level (p-value) associated with the causal interaction from $q$ to $p$, as defined in Equations 1, 2.

$$\mathbf{W}_{p,q} =1-\overline{p}_{p,q},$$ (1)

$$\overline{p}_{p,q} =\frac{1}{K}\sum_{k=0}^{K}p_{p,q}(k),$$ (2)

where $p_{p,q}(k)$ is the p-value corresponding to the MCI test from variable $q$ to $p$ at lag $k$. The estimated effective brain connectivity network is decomposed into symmetric ($\mathbf{W}_{s}$) and anti-symmetric ($\mathbf{W}_{a}$) components.

$$\mathbf{W} =\mathbf{W}_{s}+\mathbf{W}_{a}$$ (3)

where $\mathbf{W}_{s}=\frac{1}{2}(\mathbf{W}+\mathbf{W}^{\prime})$ represents mutual causal influence and $\mathbf{W}_{a}=\frac{1}{2}(\mathbf{W}-\mathbf{W}^{\prime})$ represents net causal influence. Since brain dependencies can be asymmetric, this decomposition is relevant as it separates the connectivity network into meaningful symmetric and anti-symmetric components. The symmetric component, $\mathbf{W}_{s}$, captures bidirectional interactions, reflecting mutual influences critical for cooperative neural processes. The anti-symmetric component, $\mathbf{W}_{a}$, highlights the predominant direction of causal influence, revealing the hierarchical structure and directional information flow within the brain network. This approach provides a nuanced understanding of neural activities and identifies key regions driving effective connectivity. This decomposition is attractive to neuroscientists - it is an essential tool to deepen our understanding of the hierarchical structure of the brain. The next step involves building the ranking of brain regions using Hodge Decomposition, which will allow us to quantify and rank the hierarchical influence of each region within the network.

2.2 Our Proposed Method for Ranking Brain Regions

Once effective connectivity is estimated and the net influence $\mathbf{W}_{a}$ is derived, our proposed method will rank brain regions by assigning a score to each node or brain channel based on their influence within the network. Our method is based on the HodgeRank approach, as described in [15], which has previously been utilized to rank entities such as movies and football teams. A higher ranking score indicates a region that influences other regions more than it is influenced, while a lower score suggests a region that is more influenced by others. This concept is illustrated in Figure 2.

Figure 2: Node ranking in a graph involves a gradient operator that maps ranking scores on the left to gradient flow on the right (dependence structures). The adjoint, gradient star, reverses this mapping from right to left.

To achieve the most meaningful ranking (i.e., one that minimizes violations of mutual net causal influence), we apply Hodge Decomposition to the anti-symmetric component of the connectivity network. Since $\mathbf{W}_{p,q}=-\mathbf{W}_{q,p}$, $\mathbf{W}_{a}$ is referred to as an alternating function. In the following, we first provide definitions for the gradient (grad) and curl operators, and then introduce the Hodge Decomposition Theorem.

Definition 1 (Gradient Operator)

The ’grad’ operator measures the flow of influence between nodes. For any edge $e=(p,q)\in E$, it is defined as:

$$\text{grad}(s)=s_{q}-s_{p},$$ (4)

where $s_{p}$ and $s_{q}$ represent the scores or potentials of nodes $p$ and $q$. This operator reflects the change in influence from node $p$ to node $q$.

Definition 2 (Curl Operator)

The ’curl’ operator measures the rotation strength of mutual influence among nodes. For any face $t=(p,q,r)\in T$, it is defined as:

$$\text{curl}(\phi)=\phi_{p,q}+\phi_{q,r}+\phi_{r,p},$$ (5)

where $\phi_{p,q}$, $\phi_{q,r}$, and $\phi_{r,p}$ are the edge flows between nodes $p,q$, and $r$. This operator evaluates the total rotation around a triangular face in the graph.

Given the definitions above, we can observe that the image of the gradient operator is included in the kernel of the curl operator. This relationship leads to the Helmholtz-Hodge Decomposition theorem [3].

Theorem 2.1 (Helmholtz-Hodge Decomposition)

The space of anti-symmetric matrices, $\mathbb{R}^{E}$, can be decomposed into three mutually orthogonal subspaces:

$$\mathbb{R}^{E}=\text{Im}(\text{grad})\oplus\text{Im}(\text{curl}^{*})\oplus\text{Ker}(\text{grad}^{*})\cap\text{Ker}(\text{curl}).$$

The anti-symmetric component $\mathbf{W}_{a}$ can exhibit cyclic influence in the network, which indicates that no perfect ranking can be achieved without violating at least one mutual influence relationship within each cycle. According to the Helmholtz-Hodge Decomposition theorem, any edge flow $\mathbf{W}_{a}$ can be uniquely decomposed as:

$$\mathbf{W}_{a}=-\text{grad}(s)+\text{curl}^{*}(\phi)+\mathbf{h},$$ (6)

where $\text{grad}(s)$ is the gradient component capturing hierarchical relationships, $\text{curl}^{*}(\phi)$ is the curl component indicating local cyclic influences, and $\mathbf{h}$ is the harmonic component indicating non-local or global cyclic influences.

To estimate these components, we use the least squares method for the gradient and curl components, and obtain the harmonic component as the residual.

$$\hat{s} =\arg\min_{s}\|\mathbf{W}_{a}+\text{grad}(s)\|_{F}^{2},$$ (7)

$$\hat{\phi} =\arg\min_{\phi}\|\mathbf{W}_{a}-\text{curl}^{*}(\phi)\|_{F}^{2},$$ (8)

$$\hat{h} =\mathbf{W}_{a}-\text{grad}(s)-\text{curl}^{*}(\phi).$$ (9)

The estimated score vector $\widehat{s}=[s_{1},s_{2},\ldots,s_{N}]^{\prime}$ (where $N$ is the number of nodes or time series components) derived from $\mathbf{W}_{a}$ is used to establish the topological ranking of the EEG channels (or brain regions). By minimizing the difference between $-\text{grad}(\widehat{s})$ and $\mathbf{W}_{a}$, this ranking accurately reflects net information flow patterns while avoiding inconsistencies, as $\text{curl}(\text{grad}(\widehat{s}))=0$. If $s_{p}>s_{q}$, the $p$-th time series tends to lead other time series more than the $q$-th time series, indicating that $p$ is higher in the hierarchy. Thus, $\widehat{s}$ captures the brain hierarchy, as illustrated in Figure 2. This is crucial for understanding brain function and the impact of dysfunctions such as seizures.

Our methodology involves a series of steps designed to rank brain regions based on their mutual influence. Initially, we employ the PCMCI method to estimate effective brain connectivity, utilizing causal inference techniques. Following this, the connectivity network is decomposed into symmetric and anti-symmetric components. We then apply Hodge Decomposition to the anti-symmetric component, extracting the gradient, which is essential for establishing the ranking. This allows us to derive the topological ranking of brain regions based on the gradient component. The entire process is visually represented in Figure 3.

Figure 3: Hodge decomposition pipeline for topological ranking of brain regions.

3 Experiments

To evaluate our methodology, we propose six distinct simulation scenarios, each representing different hierarchical dependency structures. These scenarios range from linear hierarchies to complex cyclic dependencies, testing the robustness and versatility of our ranking mechanism across various topologies, as illustrated in Figure 4.

• Scenarios 1 and 4: Exhibit a linear hierarchical dependence structure, allowing for clear node ranking (scenario 4 reverses scenario 1).

• Scenarios 2 and 5: Feature cyclic dependence structures, challenging meaningful ranking due to non-transitivity (scenario 5 reverses scenario 2).

• Scenarios 3 and 6: Present partial rankings with localized cycles, suggesting intermediate dependence structures (scenario 6 reverses scenario 3).

Figure 4: Visual representation of the designed dependency networks for six simulation scenarios. Each arrow indicates a dependency direction at lag one.

We generate multivariate time series data using the mixture model in Equation 10, with $P=5$ dimensions and $T=30,000$ observations, divided into 500-observation epochs to analyze temporal stability and variability in ranking. The first 5,000 observations correspond to scenario one, the next 5,000 to scenario two, and so on.

$$Y(t)=\Phi(t)Z(t)+E(t),\quad\Phi(t)=\begin{bmatrix}0&0.4&0.5&0&0\\ 0&0&0.4&0.5&0\\ 0&0&0&0.4&0.5\\ 0&0&0&0&0.4\\ 0&0&0&0&0\\ \end{bmatrix},$$ (10)

where $Y(t)$ is the observed multivariate time series, $Z(t)$ is a multivariate standard Gaussian process, and $E(t)$ is random noise, also modeled as a standard Gaussian process. The mixing matrix $\Phi(t)$ reflects each scenario’s dependence structure.

Our method follows an epoch-based approach to dynamically estimate causal relationships within the multivariate time series data in order to capture the directional dependencies and hierarchical variations. Figure 5 showcases the causal connections identified for one epoch from scenario one, highlighting both the causal pathways and the strength of these interactions.

Figure 5: Visual representation of the causal interactions within the network for scenario one as identified by PCMCI, showcasing the algorithm’s capacity to accurately capture both the direction and magnitude of causal dependencies across time lags.

By applying our proposed ranking procedure, we are able to identify and highlight nodes within the network that exert significant influence over others, providing a clear understanding of the hierarchical patterns and their dynamics, as can be seen in Figure 6. Our methodology effectively ranks the simulated time series, consistently placing node A at the top and node E at the bottom, as indicated by their respective curves, in scenario one, with the inverse pattern observed in scenario four. A node is ranked above another if its curve is consistently above the other’s. Scenarios two and five exhibit fluctuating rankings due to inherent cyclic dependencies, while scenarios three and six show partial rankings with notable positions for node D, aligning with the designed simulation scenarios.

Figure 6: Aggregate causality-based ranking scores, delineating the hierarchical structure within each scenario. The distinct patterns observed across scenarios validate the robustness of our ranking methodology.

These results validate our approach for ranking multivariate time series data based on net influence, as it robustly identifies the hierarchical structures in scenarios one and three. In contrast, scenarios two and four, characterized by cyclic dependencies, exhibit inconsistent rankings. Through these simulations, we confirm that our CBTR method effectively uncovers the hierarchical structure within dependence networks, distinguishing various dependency patterns.

4 Results

In this section, we focus on applying our methodology to a dataset of neonatal EEG recordings [24]. Neonatal seizures are a critical concern within the neonatal intensive care unit (NICU), raising numerous questions about their temporal and spatial characteristics. The dataset comprises 19-channel EEG recordings from 79 full-term infants at the Helsinki University Hospital NICU, with seizure annotations provided by three expert clinicians for each recording, spanning a median duration of 74 minutes. We concentrate our study on a few subjects, specifically analyzing epochs consistently classified as either seizure or seizure-free by all three experts. Our initial goal is to explore how seizures influence effective brain connectivity in neonates, serving as a foundational step for future detailed research on the nature of seizures and their developmental effects.

In many instances, neonatal EEG recordings reveal distinct seizure patterns, as shown in Figure 7. However, deciphering the interconnections among the 19 channels, particularly the influence of seizures, presents a complex challenge. Our approach quantifies the impact of seizures on the neonatal brain’s effective connectivity at the channel level, aiming to uncover how seizures modify brain connectivity and hierarchical organization in each subject.

Figure 7: Sample of EEG signals from channels F3, F7, Fz, and C4, showcasing epochs (1 second each) with and without seizure. Original EEG recordings are depicted in black, while the incremental changes between successive EEG readings are highlighted in cyan.

The EEG recordings typically exhibit smooth variations at the millisecond scale, leading to high autocorrelation in lagged observations, which can increase spurious causal relationships. To mitigate this, we analyze effective brain connectivity through successive differences or increments in EEG recordings, reducing the influence of lagged dependencies, as illustrated in Figure 7.

Following the approach presented in Section 2.1, we estimated the causal influences among brain channels for each subject across different time lags during seizure and non-seizure periods. Following the procedure in Section 2.2, we evaluate the changes in causal influence during seizure time by estimating ranking scores for all channels across 50 seizure epochs and compare them with scores from 50 seizure-free epochs. We used a t-test to determine the statistical significance of the ranking score differences. This analysis revealed significant changes in causality-based ranking scores for several channels for each examined subject, evidenced by very small p-values. For instance, in subjects 1 and 5, we highlight the EEG channels most significantly impacted by seizures. These findings, depicted in Figures 8 and 9, illuminate how seizures alter the hierarchical architecture of brain connectivity, offering deeper insights into their impacts and pinpointing the brain regions where seizures exert the greatest influence. These results demonstrate the efficacy of our CBTR approach in uncovering the nuanced dynamics of effective connectivity influenced by seizure events.

Figure 8: Channels with most significant shifts in the ranking scores between seizure free epochs and seizure epochs. For subject 1.

Figure 9: Channels with most significant shifts in the ranking scores between seizure free epochs and seizure epochs. For subject 5.

5 Conclusion and Discussion

In this paper, we introduced a novel framework to investigate effective brain connectivity — a domain traditionally explored through the lens of functional connectivity via persistence homology. Our proposed approach is unique in a sense that it integrates causal inference methods with Hodge decomposition, enabling the ranking of brain regions according to their mutual causal influences. This approach offers a novel perspective on the complex dynamics with effective brain connectivity.

Our simulation study demonstrates the robustness and accuracy of our methodology in uncovering and describing various hierarchical structures within multivariate dependence networks. Additionally, the application of our method to a dataset of neonatal EEG recordings highlights its practical value. Our analysis revealed novel insights into the impact of seizures on neonatal brain connectivity, specifically identifying changes in the hierarchical positions of specific brain regions in response to seizure activity.

Our future work entails deeper investigations on epileptic seizure using our proposed method. Given its potential to enhance the localization of seizure origins, this future direction holds promise for making significant contributions to seizure. By refining our understanding of seizure-induced alterations in brain connectivity, we anticipate that our methodology will offer valuable tools for clinical diagnosis and the development of targeted interventions.