*3.1. Neural Mass Models*

Neural mass models (NMM) are biologically plausible mathematical descriptions of neural mechanisms [39]. They represent the electrical activity of neural populations at a macroscopic level through a set of stochastic differential equations [40]. NMMs allow generating mildly nonlinear time series with properties that resemble the oscillatory dynamics of electrophysiological signals, such as EEG, and how they change as a result of coupling between different cortical areas. Therefore, NMMs are useful to study the behavior of brain connectivity measures that aim to capture such interactions [8,24,40,41]. Figure 1A shows a schematic representation of an NMM with two interacting cortical areas from which two signals, **x** and **y** (see Figure 1B), can be obtained. The parameters *C*<sup>12</sup> and *C*<sup>21</sup> are known as coupling coefficients, and they determine the strength of the coupling from Area 1 to Area 2, and from Area 2 to Area 1, respectively. The parameter *ν* represents the interaction lag between the two areas, while **p**<sup>1</sup> and **p**<sup>2</sup> are external inputs coming from other cortical regions.

**Figure 1.** (**A**) Schematic representation of a neural mass model. (**B**) 1 s long unidirectionally coupled time series generated by the model. (**C**) Average power spectrums peaking in the *α* and lower *β* frequency bands.

In this work, we use NMMs to generate interacting time series with known oscillatory properties in order to test the performance of the proposed phase TE estimator. In particular, we test our proposal in terms of its ability to detect directed interactions for different levels of coupling strength, under the presence of noise and signal mixing, and for bidirectional narrowband couplings. We proceed as follows: first, we set the model parameters describing Areas 1 and 2 as in [40], so as to generate signals with power spectrums peaking in the *α* (8 Hz–12 Hz) and lower *β* bands (12–20 Hz), as depicted in Figure 1C. Then, in order to generate unidirectionally coupled signals, with interactions from **x** to **y**, we set the parameter *C*<sup>21</sup> to 0 for all simulations. Also, the parameter *ν* is set to 20 ms, and the extrinsic inputs **p**<sup>1</sup> and **p**<sup>2</sup> are modeled as Gaussian noise [40]. Afterward, we generate 50 pairs (trials) of 3 s long signals, using a simulation time step of 1 ms, equivalent to a sampling frequency of 1000 Hz, for each condition in the three scenarios detailed in Sections 3.1.1–3.1.3. Next, we select a 2 s long segment from the signals, from 0.5 s to 2.5 s, and downsample them to 250 Hz. Then, we compute connectivity estimates for the simulated data in the frequency range between 2 Hz and 60 Hz, in steps of 2 Hz. After that, we obtain net connectivity values, defined as

$$
\Delta\lambda(\mathbf{x}, \mathbf{y}, f) = \lambda(\mathbf{x} \to \mathbf{y}, f) - \lambda(\mathbf{y} \to \mathbf{x}, f),
\tag{20}
$$

where *λ* stands for any of the phase-based effective connectivity measures studied, except for the PSI, in which case all subsequent analyses are performed directly on the PSI values. Lastly, for each condition in the three scenarios and at each frequency evaluated, we perform permutation tests based on randomized surrogate trials [34,42] to determine which net couplings or directed connections are statistically significant. The permutation test employed uses the trial structure of the data to generate surrogate datasets for the null hypothesis (absence of directed interactions). It does so by shuffling the data from different trials. The significance level for the tests was set to 3.3 × <sup>10</sup>−<sup>4</sup> after applying the Bonferroni correction to an initial alpha level of 0.01 in order to account for 30 independent tests, one for each evaluated frequency per condition.
