A Framework for Evaluating Dynamic Directed Brain Connectivity Estimation Methods Using Synthetic EEG Signal Generation ^†

Abstract

This study presents a method for generating synthetic electroencephalography (EEG) signals to test dynamic directed brain connectivity estimation methods. Current methods for evaluating dynamic brain connectivity estimation techniques face challenges due to the lack of ground truth in real EEG signals. To address this, we propose a framework for generating synthetic EEG signals with predefined dynamic connectivity changes. Our approach allows for evaluating and optimizing dynamic connectivity estimation methods, particularly Granger causality (GC). We demonstrate the framework’s utility by identifying optimal window sizes and regression orders for GC analysis. The findings could guide the development of more accurate dynamic connectivity techniques.

1. Introduction

The human brain consists of neurons connected by synapses, which are organized across various spatial regions and engaged in functional interactions across diverse temporal frames [1]. Electroencephalography (EEG) is a widely used method for measuring brain activity due to its excellent temporal resolution [2]. Brain connectivity analysis can be broadly divided into two main categories: structural and functional. Structural connectivity involves tracking the direction of fibers among brain regions, typically using methods like magnetic resonance imaging (MRI) or diffusion tensor imaging (DTI) [3]. Functional connectivity, on the other hand, examines the information exchange between brain regions or within a single region and can be categorized into undirected and directed measures. Undirected measures gauge the level of connectivity [4,5,6,7,8,9,10], while directed measures assess both the strength and direction of connectivity. This study focuses on directed connectivity measures.

Several methods are used for directed connectivity analysis of EEG data, including the Granger causality (GC) [11,12,13,14,15], Phase Slope Index [16,17], transfer entropy [18], and partial directed coherence [19]. Among these, GC is the most commonly employed method for both static and dynamic directed functional connectivity analysis. For dynamic analysis, GC computes connectivity based on autoregressive models using a temporal window, with the size of this window limiting the ability to describe connectivity dynamics comprehensively. Assessing dynamic capabilities is particularly challenging due to the lack of ground truth for real signals.

In this study, connectivity is defined as the directional relationship between two time-series signals, quantified using GC. GC assesses whether past values of one signal improve the prediction of another signal’s future values, beyond what could be predicted using its own past values. The assumption of GC is that the system is stationary and that a sufficient number of samples can be acquired to ensure reliable estimation.

An autoregressive model was also used in [20], where a novel bivariate irregular autoregressive (BIAR) model was proposed. The BIAR model assumes an autoregressive structure on each time series, it is stationary, and it allows the autocorrelation to be estimated, i.e., the cross-correlation and the contemporary correlation between two unequally spaced time series. It was shown that if the two signals are highly correlated, the model provides more accurate forecasts and predictions using the bivariate model compared to similar methods that use only univariate information. This feature could enable future application of the BIAR model in connectivity analyses of EEG signals, similar to GC. However, the problem of generating synthetic EEG signals with known connectivity remains unresolved in this case as well.

To facilitate the evaluation of connectivity methods for assessing dynamic connectivity [21,22], and specifically the capability of measuring changes in connectivity over time, a method for generating synthetic EEG signals is proposed. Synthetic signals can provide known reference values for directed connectivity, enabling the validation and improvement of dynamic connectivity estimation methods.

Numerous methods exist for generating synthetic EEG signals, each with its own strengths and limitations. Sinusoidal models generate signals using sine and cosine functions to mimic the rhythmic activity of brain waves. This approach is straightforward and computationally efficient, making it a popular choice for basic simulations. However, these models often lack the complexity of real EEG signals, such as non-stationarities and noise characteristics. Consequently, they may not adequately represent the true dynamics of brain connectivity. Studies indicate that while sinusoidal models can capture the basic oscillatory nature of EEG, they fail to reproduce the intricate signal characteristics observed in empirical data [23,24]. Stochastic models use random processes to generate EEG signals, often incorporating elements like autoregressive processes or Gaussian noise. These models can capture more of the variability seen in real EEG data, providing a more realistic approximation of EEG signal dynamics. However, accurately estimating the parameters of these models can be challenging. Additionally, while they can replicate certain aspects of EEG variability, they might still lack some of the intricate patterns present in real signals. Research has shown that stochastic models, such as those based on autoregressive processes, can be quite effective but may still fall short in capturing all the nuances of EEG data [23,24]. Physiologically inspired models, such as neural mass models or neural field models, aim to replicate the physiological processes underlying EEG signals. By incorporating biological constraints and mechanisms, these models offer a more realistic representation of EEG activity. However, their complexity requires extensive computational resources and detailed knowledge of the underlying physiology, which can be a significant drawback. Studies have demonstrated that physiologically inspired models can provide a closer approximation to real EEG signals by simulating neural dynamics, but they are often computationally intensive and complex to implement [23,24]. The mentioned methods do not allow us to generate signals with varying directed connectivity over time, nor do they provide the ability to control the level of such connectivity.

The proposed method aims to generate synthetic signals with realistic properties and predefined connectivity changes over time. This capability allows for the testing and development of improved methods for dynamic analysis of functional connectivity. By addressing the limitations of previous methods, the proposed approach seeks to create more accurate and reliable synthetic EEG signals, facilitating better evaluation and enhancement of connectivity estimation techniques.

The primary contribution of this paper is the development of a method to generate realistic synthetic EEG signals for testing dynamic directed connectivity estimation techniques. Our approach overcomes the limitations of existing signal generation methods by enabling control over time-varying connectivity. The paper is structured as follows: Section 2 describes the methodology for generating synthetic signals and assessing connectivity. Section 3 presents the results of applying the method to a real dataset. Section 4 discusses the findings, limitations, and future directions.

2. Materials and Methods

In this section, the GC method is reviewed, and the proposed procedure for generating realistic synthetic signals suitable for testing the dynamics of connectivity estimation methods is described.

2.1. Granger Causality

Granger causality is a statistical method employed to evaluate the causal link between two time series. It functions on the premise that if a time series, labeled as S (source or Granger-causative), influences another time series T (target), then historical values of S should contain information that aids in forecasting future values of T, beyond what can be predicted solely by past values of T. Essentially, Granger causality prediction determines whether previous values of one time series improve the forecast of another time series. When predicting the next sample in an observed time series by solely using its past values, a univariate autoregressive model is utilized. The univariate autoregressive model is characterized as

$$T\left(n\right)=\sum _{i=1}^M{\omega }_{Ti}T(n-i)+e_T\left(n\right).$$

(1)

The coefficients ${\omega }_{Ti}$ represent the autoregressive model parameters for sequence T(n), where M denotes the regression order, and e_T represents the autoregression error.

Extending this idea, bivariate autoregressive models are formulated to examine how one sequence causally influences another, defined as

$$T\left(n\right)=\sum _{i=1}^M{\omega }_{Ti}T(n-i)+\sum _{i=1}^M{\omega }_{Si}S(n-i)+e_{TS}\left(n\right),$$

(2)

where e_TS is the multivariate regression error. The coefficients $\omega$ of the model are computed by minimizing errors across N provided samples, where N exceeds M. Granger causality (GC) is subsequently defined as

$$GC=log{\displaystyle \frac{Var\left[e_T\left(n\right)\right]}{Var\left[e_{TS}\left(n\right)\right]}}.$$

(3)

In the analysis of directed functional connectivity, predicting Granger causality serves as a robust method to deduce the directional influences among various brain regions using their time-series data. This approach evaluates whether prior activity in one brain area (the source) offers predictive insights into future activity in another (the target), surpassing what can be anticipated from the source’s own past activity alone. Utilizing Granger causality prediction provides valuable insights into both the static and dynamic interactions and causal links within intricate systems, facilitating the comprehension of predictive associations and decision-making mechanisms.

2.2. Generation of Synthetic Signals

Once the regression models are understood, they can be employed to transform real signals into synthetic ones. This allows us to regulate the influence of source signals by introducing an additional mixing parameter, which can vary arbitrarily over time and thus govern the dynamics of connectivity. Before the analysis, the stationarity of the EEG signals was verified using the augmented Dickey–Fuller (ADF) test. The ADF test results indicated that the signals (from electrodes FC1 and F1 after offline preprocessing) were stationary (p-value for both is equal 0.0010), meeting the requirements for GC analysis. However, a limitation lies in the framework’s assumption of stationarity within selected intervals, which may reduce accuracy when simulating highly non-stationary or rapidly fluctuating EEG signals.

Initially, a method for generating signals needs to be established in which the maximum connectivity equals that of the source signals, derived from independent sources. For example, consider two actual source signals, T and S, and construct their regression models to estimate their respective univariate and bivariate regression errors, denoted as $e_T$ and $e_{TS}$, alongside their connectivity measure GC_TS. In this study, two signal segments are extracted from different time periods to ensure minimal correlation, providing a suitable baseline for generating synthetic signals with controlled connectivity. These subsections, named T′ and S′, are expected to exhibit minimal correlation and low Granger causality. Using these sequences, a new sequence R₁ is generated that resembles signal T′ but also incorporates aspects of S′, employing the following regression model:

$$R_1\left(n\right)=\sum _{i=1}^M{\omega }_{Ti}{T}^{\prime }(n-i)+\sum _{i=1}^M{\omega }_{Si}{S}^{\prime }(n-i)+e_{TS}^{\prime }\left(n\right).$$

(4)

Incorporating the regression error ${e^{\prime }}_{TS}$, which equals $e_{TS}$ in the subsection of signal T′, enhances the realism of the generated sequence in terms of predictability. The resultant signal R₁ is anticipated to exhibit Granger causality (GC_R1S) with respect to S, similar to GC_TS.

However, the drawback of this model lies in its inability to dynamically adjust the actual connectivity. To obtain a signal that is unrelated to S, one could employ a univariate regression model. Nevertheless, the resulting outcome would mirror the target signal T′:

$$R_2\left(n\right)=\sum _{i=1}^M{\omega }_i{T}^{\prime }(n-i)+e_T^{\prime }\left(n\right)=T^{\prime }\left(n\right).$$

(5)

Here, ${e^{\prime }}_T$ equals $e_{T^{\prime }}$. By combining realistic signals resembling both high and low connectivity, such as R₁ and R₂, the connectivity can be adjusted dynamically by mixing them together:

$$R^{\prime }\left(n\right)=K\left(n\right)·R_1\left(n\right)+\left(1-K\left(n\right)\right)·R_2\left(n\right)$$

(6)

The parameter K represents the level of connectivity between the synthetic signal and the source. When K = 0, no connectivity is imposed (GC = 0), whereas K = 1 corresponds to the maximum connectivity level derived from the original signals.

Therefore, utilizing the pair {R′, S′} will allow us to explore dynamic connectivity methods in future studies.

3. Results

The proposed method for generating signals is illustrated using the EEG Motor Movement/Imagery Dataset [25]. Specifically, the S001R01 baseline recording with eyes open is utilized. These data were sampled at 160 Hz (f_s) and spanned a duration of 61 s. The selected order for autoregressive models, M, was 19, which is equivalent to 120 ms, a time horizon relevant for EEG connectivity analysis. This selection was based on prior studies indicating its relevance for capturing EEG connectivity [26,27,28]. We also experimented with other regression orders but found that 19 provided a good balance between model accuracy and computational efficiency. The GC connectivity matrix was computed (see Figure 1), revealing a high GC connectivity of 0.92 between electrodes FC1 (T) and F1 (S).

To generate signals, subsections from electrodes FC1 and F1 were selected: from 0 to 14.4 s for T′ (recorded by FC1) and from 28.8 to 43.2 s for S′ (recorded by F1). The GC value between these unrelated signal segments was calculated to be 0.008.

The temporal connectivity parameter K was set to 1 during the first interval of R′ (0 to 4.8 s) and the third interval (9.61 to 14.4 s), while it was set to 0 in the second interval (4.81 to 9.6 s). Using the proposed method (Equation (6)) and the previously obtained autoregressive model parameters, a synthetic signal was generated. This synthetic signal R′ and the source signal S′ were then used for functional connectivity analysis using GC (see Figure 2).

The GC values obtained for the three defined intervals were 0.77, 0.032, and 0.77, respectively (referred to below as reference functional connectivity (RFC) values). As anticipated, the RFC value for the second interval is low, whereas those for the first and third intervals are high. However, these values are still lower than the GC of the original signals (calculated for the entire period of the signal).

To demonstrate the applicability of the signal generation method for analyzing dynamic connectivity estimation methods, dynamic GC [28,30,31] with a sliding window approach was employed. Here, GC was estimated for the generated signal, which mimics the modified FC1 signal relative to the source signal F1. Figure 3a and Figure 3b depict the results obtained using window sizes of 2 s and 400 ms, respectively.

It is evident that the estimated connectivity varies over time, showing expected deviations from the RFC, particularly during transitions. Figure 3, shows that there are different limitations in the estimation of the dynamic GC for both wide and narrow window sizes. For wider windows, there is a good agreement between RFC and the estimated values when there is constant reference connectivity. When an abrupt change occurs in RFC, the estimated dynamic GC slowly transitions to a new value, as expected. The wide window limits the rate of change, as observed in Figure 3a.

For narrow windows, a noisier estimation of dynamic GC is expected, along with a faster response to abrupt changes in RFC. This underscores the limited capability of GC measures in identifying the precise times of connectivity changes and highlights the necessity for further research in this direction.

To further investigate the aforementioned limitations, dynamic Granger causality (GC) was calculated for various window sizes, ranging from a minimal window size to 3000 ms with a step of 25 ms, and for different regression orders (during estimation), M_GC, ranging from 5 to 35 with a step of 1. The minimal window size equals (M_GC + 1)/f_s. During signal generation, a regression order of 19 was used, and the temporal connectivity parameter K was set to 1. The root mean square error (${\epsilon }_{RMS}$) was employed as the measure of goodness of fit between the RFC and estimated dynamic GC:

$${\epsilon }_{RMS}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{(RFC-GC)}^2},$$

(7)

where RFC represents the reference Granger causality value calculated over a predefined time period, and GC represents the estimated dynamic GC value. The root mean square error (${\epsilon }_{RMS}$) was chosen as the primary metric for evaluating dynamic GC estimates because of its sensitivity to deviations of varying magnitudes from reference connectivity values, which offers a more comprehensive measure of estimation accuracy. While we also considered mean absolute error (MAE), ${\epsilon }_{RMS}$ provided better differentiation between the tested methods due to its squared error component that amplifies larger discrepancies. The performance of the model improves with increasing regression order up to a certain point, as higher orders can capture more complex relationships in the data. However, beyond order 9, we observed diminishing returns and an increased risk of overfitting, which negatively affected the generalizability of the model. The selection of order 9 as optimal was based on minimizing the root mean square error across different window sizes, balancing model complexity with estimation accuracy. The change in ${\epsilon }_{RMS}$ error values with respect to different window sizes and regression orders is presented in Figure 4. High ${\epsilon }_{RMS}$ values are observed when using narrower windows and lower regression orders. Figure 5 illustrates the behavior of the estimated dynamic GC value for the minimum ${\epsilon }_{RMS}$ value (window size of 875 ms, M_GC of 9).

Figure 6 shows the optimum window size considering ${\epsilon }_{RMS}$ for different M_GC and respective ${\epsilon }_{RMS}$. It can be seen that the optimum window size increases as the order M_GC is increased. Moreover, the minimum value of ${\epsilon }_{RMS}$ is obtained for order 9, and only slightly increases for further increases in the order.

The transition time (TT) metric is calculated as the duration required for the dynamic GC measure to change from 10% to 90% (rise time) of a step change in the reference connectivity value (and the fall time, the time needed to change from 90% to 10%). This metric provides insight into the responsiveness of the connectivity estimation method, which is crucial for real-time monitoring applications, such as brain–computer interfaces (BCIs) and neurofeedback, where rapid detection of changes in brain connectivity can enhance system performance. The change in TT values with respect to different window sizes and regression orders is presented in Figure 7. As before, dynamic GC was calculated for various window sizes ranging from the minimal window size to 3000 ms with a step of 25 ms, and for different regression orders (during estimation) M_GC ranging from 5 to 35.

Examining Figure 7 and analyzing the TT behavior for all analyzed combinations, it is evident that a smaller TT value generally corresponds to a smaller ${\epsilon }_{RMS}$. An illustration of the estimated dynamic GC value for the minimum TT value is given in Figure 8. Figure 9 shows the optimum window size considering TT for different M_GC and respective TT. In this case, it can also be seen that as the order increases beyond 11, the optimum window size also increases as the order M_GC increases. Contrary to Figure 6, TT does not monotonically increase with increasing order M_GC.

Considering Figure 6 and Figure 9, the window size $M_{GC}$ needs to be optimized according to both the error ${\epsilon }_{RMS}$ and the transition time $TT$. We propose the criterion of the product of both ${\epsilon }_{RMS}$ and $TT$. The product with respect to different window sizes and orders $M_{GC}$ is shown in Figure 10.

An illustration of the estimated dynamic $GC$ value for the minimum product $TT$ and ${\epsilon }_{RMS}$ value is given in Figure 11. It can be observed that $TT$ predominantly influences dynamic GC behavior (due to the similarity between Figure 8 and Figure 11). This can be explained by the small variations in ${\epsilon }_{RMS}$ for the optimum window size after order 9. Finally, the optimum window size considering the $TT$ and ${\epsilon }_{RMS}$ product is presented in Figure 12.

4. Conclusions

This study outlines the process of creating a synthetic signal designed for evaluating the dynamic characteristics of directional connectivity techniques. The signal is generated by reconstructing data from real signals using coefficients from regression models. By employing the approach proposed in this article, the impact of source signals can be controlled through the incorporation of a variable temporal connectivity parameter K, which adjusts dynamically over time, thereby defining the fluctuation in actual connectivity dynamics. This ensures that the generated synthetic signals closely mirror real-world scenarios.

During the evaluation of dynamic connectivity methods, the connectivity between the generated signal R′ and the source signal S′ is assessed. These results are compared against the weighting parameter K. This process was demonstrated through dynamic GC estimation using a sliding window approach. The challenge lies in selecting an appropriate window size: larger windows typically provide lower temporal resolution of dynamic connectivity estimates, while narrower windows can compromise the accuracy of autoregressive models, leading to increased variability in estimated dynamics.

Identifying the optimal window size or comparing different methodologies requires knowledge of the genuine changes in connectivity, which are precisely known only in the case of synthetically generated signals. Hence, it is crucial that the generated signals accurately reflect properties obtained from real signals.

The proposed methodology satisfies both criteria, ensuring that the synthetic test signals are realistic and suitable for validating and advancing dynamic connectivity techniques. Therefore, the proposed method for generating synthetic test signals is well suited for the development and evaluation of dynamic connectivity methods.

Additionally, in this study, a comprehensive analysis was conducted of dynamic GC across varying window sizes and regression orders M_GC. By evaluating the root mean square error between reference and estimated dynamic GC values, the optimum window size was identified. Our findings demonstrate that, in our case, the window size increases with order size; the minimum error is obtained for order 9, and it only slightly increases with further increase in the order.

Furthermore, the transition time of GC change was analyzed by calculating rise and fall times between 10% and 90% of connectivity thresholds. The results show that the smallest TT value generally corresponds to higher ${\epsilon }_{RMS}$. Therefore, an optimal solution that minimizes the TT and ${\epsilon }_{RMS}$ product was further determined.

These findings provide a detailed framework for selecting appropriate window sizes and regression orders in dynamic $GC$ analysis, enhancing the reliability of connectivity studies. In addition, this methodology could be applied to applications in other domains where connectivity analysis of signals is relevant. Future work should focus on further refining these parameters and exploring their applicability in various experimental and real-world scenarios. For real-time applications, where responsiveness is critical, smaller windows can provide faster adaptation to abrupt changes, although with potentially increased noise in the estimates. Increased computational complexity and cost of methods based on varying window sizes should be analyzed in more detail for real-time applications.

In addition, in future work we will research the ways to compensate for the poor temporal resolution with a wide window and poor accuracy with a narrow window by using methods for determining adaptive window width, such as intersection of confidence intervals (ICIs), relative intersection of confidence intervals (RICIs) and single-scale time-dependent (SSTD).