*2.1. Simulation Procedure*

The simulation framework developed by [1] was used to generate simulated EEG data originating from three dipoles. This data generating process, as well as the forward and inverse problems, were implemented in MATLAB (2020). Figure 1 shows the data generation procedure. The standard length of each generated series was 1500 time steps for Ground Truth 1 and 2.

**Figure 1.** Simulation Procedure followed by Connectivity Estimation. TCDF = Depthwise Separable 1D Temporal Causal Discovery Framework, LSTM-NUE = Long Short-Term Memory with Non-Uniform Embedding, Conv2D = 2D Convolutional Network, TRGC = Time-Reversed Granger Causality. Coordinates in the Ground Truths denote MNI-coordinates.

In Ground Truth 1, three fixed dipoles were used with the directionality of the connections as well as their strength being imposed (Figure 2), a strategy used before [20–22]:

$$\begin{array}{c} \mathcal{X}\_1(t) = 0.5\mathcal{X}\_1(t-1) - 0.7\mathcal{X}\_1(t-2) + c12(t)\mathcal{X}\_2 + \in\_1(t) \\ \mathcal{X}\_2(t) = 0.7\mathcal{X}\_2(t-1) - 0.5\mathcal{X}\_2(t-2) + 0.2\mathcal{X}\_1 + c23(t)\mathcal{X}\_3(t-1) + \in\_2(t) \\ \mathcal{X}\_3(t) = 0.8\mathcal{X}\_3(t-1) + \in\_3(t) \end{array} \tag{1}$$

with *X*1, *X*<sup>2</sup> and *X*3*,* three electrical sources contributing to the simulated scalp-EEG signals and with:

$$\begin{aligned} c12(t) &= 0.5 \frac{t}{L} \text{ if } t \le \frac{L}{2}, \text{ c12(t)} = 0.5 \frac{L-t}{\frac{L}{2}} \text{ if } t > \frac{L}{2} \\ c23(t) &= 0.4 \text{ if } t < 0.7L, \text{ c23(t)} = 0 \text{ if } t \ge 0.7L \end{aligned} \tag{2}$$

*L* = length of the generated time series (*L* = 1500), *t* = the current time step and E = uncorrelated white noise, varying with time. We further assume an EEG cap with 108 electrodes.

**Figure 2.** Ground Truth 1 with three fixed dipoles.

For Ground Truth 2, we considered two fixed, one moving dipole and only one true connection (Figure 3) and focused on the presence or absence of this connectivity as well as its directionality:

$$
\begin{bmatrix} X\_{\mathfrak{s}}(t) \\ X\_{\mathfrak{r}}(t) \\ X\_{\mathfrak{n}}(t) \end{bmatrix} = \sum\_{p=1}^{p} \begin{bmatrix} a\_{11}(p) & 0 & 0 \\ a\_{21}(p) & a\_{22}(p) & 0 \\ 0 & 0 & a\_{33}(p) \end{bmatrix} \begin{bmatrix} X\_{\mathfrak{s}}(t-p) \\ X\_{\mathfrak{r}}(t-p) \\ X\_{\mathfrak{n}}(t-p) \end{bmatrix} + \begin{bmatrix} \varepsilon\_{1}(t) \\ \varepsilon\_{2}(t) \\ \varepsilon\_{3}(t) \end{bmatrix} \tag{3}
$$

with *Xs*, the moving dipole, as a sender, and two fixed dipoles, with *Xr* the receiver and *Xn* the fixed non-interactive dipole, and *a*ij (*p*), i, j - {1, 2, 3} and *p* - {1, ... , P} the coefficients with *a*<sup>21</sup> the coupling strength between sender and receiver. All *a*ij are randomly picked from the interval [0.3, 1]. Finally, is uncorrelated, biological, white noise.

**Figure 3.** Ground Truth 2 with two fixed, one moving dipole.

The moving dipole (the sender) changes location (far, deep, close, superficial) at every iteration, with a total of 1004 iterations. The maximum time lag t is two. The reason for this ground truth is that the sender can be located at really challenging locations (too close to one of the other dipoles or very deep in the brain).

Two conditions were created for both ground truths: one condition consisted of three superficial dipoles far away from each other, while the other consisted of three dipoles located "deep" in the brain, but each dipole was still positioned far away from the other dipoles. The corresponding MNI-coordinates of the two fixed dipoles that Ground Truth 1 and 2 have in common are depicted in Figure 4. The full set of coordinates of Ground Truth 1 (including the coordinates of the third fixed dipole) is denoted in Figure 1.

In Ground Truth 2, the first two coordinates are the same as in Ground Truth 1, for each dipole condition, while the third dipole moves throughout the brain as described above. The Far–Superficial versus Far–Deep configurations indicate (relative) distances: "deep" denotes a distance from the origin (located at the anterior commissure) <6 cm and "superficial" >6.5 cm. The distance between dipoles is evaluated as "far" if the relative distance to the other dipoles exceeds 8 cm.

**Figure 4.** Locations (MNI-coordinates) of the dipoles the Ground Truths have in common. Left: Far–Superficial fixed dipoles. Right: Far–Deep fixed dipoles.

As an additional check for robustness of source localization, noise sources were added as a background activity. These were modeled using pink noise, also called 1/f noise, and created by scaling the amplitude spectrum of random white Gaussian noise with the factor1/f using the Fourier transform and its inverse.

After generating these noise sources, the forward problem is construed:

$$Y = LX + \mathfrak{e},\tag{4}$$

where *Y* denotes the scalp-recorded potentials, *X* represents the electrical sources in the brain (the dipoles), "*e*" is measurement noise (electrode noise) and *L* is the head volume conductor model (also called the leadfield matrix). The leadfield matrix determines how the activity flows from dipoles to electrodes. In this work, the New York Head model [23] was used.

The pink noise and the source activity are then projected onto the scalp, after which they are summed:

$$Y^{\text{train}}(t) = \gamma \times \frac{Y^{\text{active}}(t)}{||Y^{\text{active}}(t)||\_{FRO}} + (1 - \gamma) \times \frac{Y^{\text{noise}}}{||Y^{\text{noise}}(t)||\_{FRO}} \tag{5}$$

*Yactive* and *Ynoise* refer to the scalp-projected source signals and pink noise activity, respectively; both are scaled by dividing them by their Frobenius norm (||*Yactive* (*t*)||*FRO*, ||*Yactive* (*t*)||*FRO*). The Signal-to-Noise Ratio (SNR) is computed for all dipoles simultaneously and set to 0.9 (γ = 0.9).

Next, white noise (spatially and temporally uncorrelated activity) is added to *Ybrain* to simulate electrode noise, resulting in Equation (6) where *Ymeasurement* represents the simulated EEG signal. Again γ = 0.9 is imposed as Signal-to-Noise Ratio:

$$Y^{\text{measure}, \text{current}}(t) = 0.9 \times \frac{Y^{\text{train}}(t)}{||Y^{\text{train}}(t)||\_{FRO}} + 0.1 \times \frac{Y^{\text{max\\_noise}}}{||Y^{\text{meas\\_noise}}(t)||\_{FRO}} \tag{6}$$

Afterwards, the simulated scalp-EEG data are source-reconstructed using exact lowresolution brain electromagnetic tomography (eLORETA) [24]. There have also been improvements to eLORETA, such as Sparse eLORETA, which uses a masking approach to improve the source localization density [25]. The eLORETA method is a discrete, threedimensional (3D), linear, weighted minimum norm inverse solution [24]. In the absence of noise, an exact zero-error localization accuracy can be obtained with eLORETA, but this does not hold for noisy data, as was shown in a study comparing both scenarios [26]. The MATLAB implementation of the eLORETA algorithm (mkfilt\_eloreta2.m) from which spatial filters are obtained was developed by G. Nolte and is available in the MEG/EEG Toolbox of Hamburg (METH; https://www.uke.de/english/departments-institutes/institutes/ neurophysiology-and-pathophysiology/research/research-groups/index.html, accessed

on 22 December 2021). As the input, it takes the leadfield tensor (i.e., the head model file N\*M\*P containing N channels, M voxels, and P dipole directions) as well as a regularization parameter gamma (set to 0.01); as the output, an N\*M\*P tensor A of spatial filters is returned.
