Synthetic neuronal datasets for benchmarking directed functional connectivity metrics

Autoregressive modeling

Autoregressive (AR) modeling is the simplest and most straightforward method for simulating neuronal datasets. By specifying the MVAR equations and parameters it is possible to simulate datasets from networks with different number of nodes, topology, noise, interaction delays and duration. Although it has been shown that the AR model’s impulse response function has a transfer function that resembles the transfer function of a simple physiological model of EEG generation (Blinowska & Franaszczuk, 1989), as most connectivity metrics are parametric, estimating an MVAR model for datasets generated by linear AR processes might not pose the required challenge expected when presenting a novel connectivity metric. Nevertheless, this method has been used in many well-known works either standalone or as initial validation, followed by real neurobiological or neuroimaging datasets. It has been used to demonstrate that partial directed coherence (PDC) can correctly identify the directed connectivity in multivariate datasets from elaborate networks with reciprocal connections and interaction delays with different magnitudes and durations (Baccalá & Sameshima, 2001). Similarly, the direct directed transfer function specific sensitivity to directed causal effects was also demonstrated with data from a MVAR model with unitary delays between variables (Korzeniewska et al., 2003). Equivalent models have been used to further study metrics like directed transfer function (DTF) (Eichler, 2006), improvements on spectral GC (Chen, Bressler & Ding, 2006) or demonstration of connectivity metrics in a vastly used toolbox (Seth, 2010).

MVAR models have also been key in data simulation for the still open debate (Smith et al., 2010; Valdes-Sosa et al., 2011; Seth, Chorley & Barnett, 2013) concerning the use of directed connectivity metrics to assess causal influences from fMRI BOLD signals. Roebroeck used synthetic data generated from a bi-dimensional first-order AR process as LFPs which were later transformed into BOLD signals in which causality was assessed with GC for different experimental parameters (Roebroeck, Formisano & Goebel, 2005). Later, Schippers used the same MVAR model to understand the effect of hemodynamic lag opposing the interaction neuronal delay in group analysis (Schippers, Renken & Keysers, 2011) and Barnett used a similar first order model to solve GC analytically after digital filtering (Barnett & Seth, 2011). Other studies used MVAR models to study the effectiveness of GC applied to cluster sets obtained with independent component analysis (Londei et al., 2006) and principal component analysis with partial canonical correlation analysis (Sato et al., 2010). Others have used AR networks where real LFP propagate according to predefined weights and delays to study the dependence of GC metrics to experimental parameters (Deshpande, Sathian & Hu, 2009; Rodrigues & Andrade, 2014).

An MVAR model of order p can be defined by $\mathit{X}\left(t\right)=\sum _{j=1}^p{\mathit{\Phi }}_j\mathit{X}\left(t-j-\mathit{D}\right)+\mathit{w}\left(t\right)$ where X(t) is the multivariate time-series, Φ_j is the MVAR coefficient matrix for time-lag j, p is the model order, D is the interaction delay matrix and w(t) is the innovation matrix consisting in independent white noise. With these parameters, it is possible to define all the aforementioned network properties.

Here we also consider the possibility of defining the spectral peak frequency of the AR process that generates the activity of each variable. Considering the AR model for variable n, without accounting for inter-variable interactions, $X_n\left(t\right)=\sum _{j=1}^p{\phi }_{nj}{X}_n\left(t-j\right)+w_n\left(t\right)$ a spectral peak frequency can be identified if the model is causal which requires all the roots of the AR polynomial ${\phi }_{nj}\left(z\right)=1-\sum _{j=1}^p{\phi }_{nj}{z}^j$ to be outside the unit circle. As proved by Jiru (2008), when the absolute value of the AR coefficient for the time-lag p is close to unity, |φ_p| ≈ 1, the spectral peak frequencies have approximately the same values as the arguments of the roots of the AR polynomial. Other properties, especially for models with p = 2, 3 can be found in Jiru (2008). The possibility to manipulate the spectral peak frequencies of each variable in the MVAR model adds up to the previous parameters so now it is possible to generate synthetic data for a given number of variables with a known connectivity pattern (reflected in the inter-variable AR parameters distribution), interaction delays (reflected in D), interaction strength (reflected in the inter-variable AR parameters value), SNR (reflected in ${\sigma }_w^2$) and peak frequency (reflected in the intra-variable parameters).

Although the inter-variable AR coefficients define the interaction strength between variable, it is also possible to compute the required values for these parameters in order to obtain a desired theoretical GC as done in Barnett & Seth (2011). This is done for each variable pair by analytically deriving the expression for their theoretical GC value, based in the MVAR coefficient matrix, transfer function and noise covariance matrix, and solving it for the inter-variable coefficients.

For our simulations, a bivariate AR(2) model was created with variable 1 exerting causal effect in variable 2 at a specific frequency ω, with variable intensity F_X1→X2 and delay d₂₁ following the expression: (1)${\displaystyle \begin{array}{l}{X}_1\left(t\right)={\phi }_{11}\left(2\right)X_1\left(t-2\right)+{\phi }_{11}\left(1\right)X_1\left(t-1\right)+w_1\left(t\right)\hfill \\ X_2\left(t\right)={\phi }_{22}\left(2\right)X_2\left(t-2\right)+{\phi }_{22}\left(1\right)X_2\left(t-1\right)+{\phi }_{21}\left(1\right)X_1\left(t-d_{21}\right)+w_2\left(t\right).\hfill \end{array}}$ Here w₁ and w₂ are the model’s uncorrelated, zero mean, unit variance, white Gaussian innovation processes. Causality was chosen to occur always in the gamma band, more specifically at 33 Hz so, the AR parameters for variable 1 AR process were chosen so the spectral peak occurs at this frequency. Following the relationship ω = ± arccos(φ₁₁(1)(φ₁₁(2) − 1)/(4φ₁₁(2))) between spectral peak frequency and AR parameters found in Jiru (2008) for AR(2) models, variable 1 parameters must be φ₁₁(1) = 1.337 and φ₁₁(2) = − 0.98 for a spectral peak to exist at 33 Hz.

The causal effect occurs due to φ₂₁(1). To establish a relationship between this parameter and the consequent GC value, spectral GC can be solved analytically for the model (1) similarly to what occurs in Barnett & Seth (2011). Skipping intermediate steps this results in: (2)${\displaystyle {\mathrm{F}}_{X_1\to X_2}\left(\omega \right)=ln\left(1+\frac{{\phi }_{21}{\left(2\right)}^2}{{\phi }_{11}{\left(1\right)}^2+{\phi }_{11}{\left(2\right)}^2+2{\phi }_{11}\left(1\right)cos\left(\omega \right)\left({\phi }_{11}\left(2\right)-1\right)-2{\phi }_{11}\left(2\right)cos\left(2\omega \right)+1}\right).}$ Solving (2) for φ₂₁(1) allows an AR(2) model to be built with the desired causality. Although LFPs simulation might seem redundant (GC estimation from data modeled with parameters obtained by the analytical solution of the spectral GC formulation) this data is helpful to analyze the effects of the forward EEG and BOLD models and to benchmark or demonstrate other causality metrics. Varying theoretical causality is modeled (F_X1→X2(33 Hz) = [from 0 to 5 in steps of 0.5], d₂₁ = 20 ms) at 33 Hz for different interaction delays (F_X1→X2(33 Hz) = 5, d₂₁ = [4 ms, and 20 to 100 ms in steps of 20 ms]). Initial values were randomly assigned, the simulations produced 60 s of data and the first 20 s were discarded to remove transient effects, data was generated at 1k Hz and subsampled to 250 Hz. Each 60 s of data required ∼1 s of simulation time.

Coupled oscillators

Coupled oscillator networks are a more realistic way to represent the dynamics between neuronal populations. Each node usually represents a population of excitatory and inhibitory spiking neurons that exhibit oscillations with varying levels of synchrony in specific frequency ranges. The network can comply with one of two following dynamical regimes: a synchronous state with self-sustained oscillations (Börgers & Kopell, 2003) or in an asynchronous state with transient oscillations (Mattia & Del Giudice, 2002). These networks have been used to show that DTF can be interpreted within the GC framework (Kamiński et al., 2001) by modeling interacting cortical columns with excitatory and inhibitory populations with delay-coupled nonlinear stochastic differential equations (SDE). Similar SDE delayed networks, called Wilson–Cohen models, are also used in Deco, Jirsa & McIntosh (2011) to study the dynamics of simulated resting state networks (RSNs). A similar but more detailed dynamic causal model is used to study GC (Friston et al., 2014) where each cortical column is comprised of pyramidal and inhibitory cells from supra-granular layers, excitatory spiny cells in granular layers and deep pyramidal cells in infra-granular layers.

A general network of coupled second order differential equations similar to Kamiński et al. (2001), represents a node as a system with delay-coupled nonlinear SDE for the excitatory population defined by: (3)${\displaystyle {\ddot{x}}_n+\left(a+b\right){\dot{x}}_n+ab{x}_n\left(t\right)=-k_{ei}S\left(y_n\left(t\right)\right)+\sum _{m=1}^N{k}_{nm}S\left(x_m\left(t-d_{nm}\right)\right)+I_n\left(t\right)+w_{x_n}\left(t\right)}$ and inhibitory population defined by: (4)${\displaystyle {\ddot{y}}_n+\left(a+b\right){\dot{y}}_n+ab{y}_n\left(t\right)=k_{ie}S\left(x_n\left(t\right)\right)+w_{y_n}\left(t\right).}$ Variables x_n and y_n represent the LFPs of excitatory and inhibitory populations respectively of node n, k_ei and k_ie their respective coupling coefficients, k_nm is the coupling coefficient from node n to node m, d_nm is the delay from node m to node n, a and b are time constants that define the rate at which activity decays without input and, w_xn and w_yn are the independent white noise processes for nodes n and m respectively. I_n(t) is the external input to the excitatory population and S is the following sigmoidal function for a modulatory parameter S_m: (5)${\displaystyle S\left(x,S_m\right)=\left\{\begin{array}{ll}{S}_m\left(1-{\text{e}}^{-\left({\text{e}}^x-1\right)/S_m}\right)\hfill & \text{if }x>-u_0\hfill \\ -1\hfill & \text{if }x\le -u_0\hfill \end{array}\right.}$ ${\displaystyle u_0=-ln\left(1+ln\left(1+1/S_m\right)\right).}$ Following the notations in Freeman (1987), (3) and (4) is a coupling between KO_i and KO_e subsets; hence, each node can be seen as a KI_e,i set. In this work, these parameters were used with the same values as in Freeman (1987) and Kamiński et al. (2001): a = 0.22/ms, b = 0.72/ms, k_ei = 0.4, k_ie = 0.1, S_m = 5 and the independent Gaussian white noise processes had zero mean and 0.04 variance. The system’s numerical solution was approximated with a fourth order Runge–Kutta method (delays where linearly interpolated for the intermediate increments) using a time-step of 0.1 ms, with the noise term being integrated with the Euler method using the same time-steps. The LFPs were initialized as zero for as long as the longest delay present in the simulation required.

A network of two KI_e,i sets, defined by the interactions in (3) and (4), was simulated for 60 s with varying values of weak and strong coupling (k₂₁ = [0, 0.1, 0.2, 0.5, 0.7, 1, 3, 5, 10, 15, 22, 30], d₂₁ = 20 ms) and interaction delay (k₂₁ = 30, d₂₁ = [4 ms and 20–100 ms in steps of 20 ms]) between sets 1 and 2. The first 20 s were discarded to remove transient effects, and data was generated at 10k Hz and subsampled to 250 Hz. External input I_n(t) lasted 1 ms and had 1% probability of occurrence for each set. Each 60 s of data required ∼181 s of simulation time.

As Friston et al. (2014) concludes that GC is not appropriate for data generated by delay-coupled oscillators with unstable modes, and because self-sustained oscillations occur in large scale simulations (Deco, Jirsa & McIntosh, 2011), we also focus our simulations on networks functioning in the synchronous mode. In this regime, since oscillators show a limit cycle phase space trajectory, this phase can be modeled by a single dynamical variable reducing the former models to a simpler phase oscillators where it is possible to define the oscillating frequency. Cabral et al. (2011) modeled RSNs with several delayed-phase oscillator networks using the Kuramoto model (Kuramoto, 1984), layered for each frequency of interest. Therefore, in this work, we use the same network building scheme consisting of several stacked layers of two coupled oscillators (one layer for each natural frequency) where each phase variable θ_n is governed by the following dynamical equation: (6)${\displaystyle {\dot{\theta }}_n={\omega }_n+k\sum _{m=1}^N{C}_{nm}sin\left({\theta }_m\left(t-d_{nm}\right)-{\theta }_n\left(t\right)\right)+w_{{\theta }_n}\left(t\right).}$ Here ω_n is the angular frequency of each oscillator (rad/s), C_nm is the relative coupling coefficient from node m to node n, k is a global coupling coefficient and the remaining variables represent the same parameters as in the previous dynamic equations. The neuronal activity can be obtained as the firing rate of the population of neurons represented by each oscillator. Following the procedure in Cabral et al. (2011), the firing rate is the sine function r_n = sin(θ_n(t)) of phase variable.

For our simulations, the system’s numerical solution was approximated with the Euler–Maruyama (Kloeden & Platen, 1992) method using a time-step of 0.1 ms. Phases were initialized randomly in the initial instants corresponding to the longest delay. A network of two nodes was simulated for 60 s (first 20 s were discarded to remove transient effects) by having three independent layers of pairs of phase-delayed coupled Kuramoto oscillators where at each layer share the same natural frequency and follow the dynamics in (6). Layers 1, 2 and 3 oscillate at 5, 33 and 60 Hz, respectively, and only layer 2 has a coupling coefficient between oscillators different than 0. Therefore, phase-delayed coupling can only occur at the 33 Hz and from node 1–2. Coupling varied from weak to strong (k = [0, 0.1, 0.2, 0.5, 0.7, 1, 3, 5, 10, 15, 22, 30], d₂₁ = 20 ms) and interaction delays followed the same values as the previous simulations (k = 30, d₂₁ = [4 ms and 20–100 ms in steps of 20 ms]). Each 60 s of data required ∼55 s of simulation time.

Izhikevich columns

The former oscillator models can be simulated in greater detail by modeling their constituting spiking units individually. In the context of testing directed connectivity, this modeling was used in Seth, Chorley & Barnett (2013) to understand how GC’s sensitivity is affected by vascular latencies opposing the neuronal lag in BOLD signal time-series sampled with decreasing time repeat (TR). Besides spiking units, their synapses were also modeled with explicit NMDA, AMPA, GABAa and GABAb conductances and short-term plasticity (STP). The Izhikevich spiking model is able to produce several firing patterns observed in real neurons without the computational demand of more biophysically realistic models like the Hodgkin–Huxley, which makes it appropriate for larger scale simulations (Izhikevich, 2003). This model replaces the bio-physiologic meaning of the Hodgkin–Huxley model’s equations and parameters by a topologically equivalent phase dynamic modeled with two ordinary differential equations and four parameters: (7)${\displaystyle \begin{array}{l}\dot{v}=0.04v^2+5v+140-u+I\hfill \\ \dot{u}=a\left(bv-u\right).\hfill \end{array}}$ Also, the additional after-spike resetting when v ≥ + 30 mV, v = c and u = u + d.

The variable v represents the neuron’s membrane potential and u its recovery variable. Parameters a, b, c and d are defined in order to implement one of the twenty spiking neurons shown in Izhikevich (2004). Similarly to Seth, Chorley & Barnett (2013), in this work we simulate cortical columns with two different neurons: regular spiking excitatory pyramidal neurons (a = 0.02, b = 0.2, c = − 65, d = 8) and fast spiking inhibitory interneurons (a = 0.1, b = 0.2, c = − 65, d = 2). Due to the model’s simplicity, it is also computationally amenable to compute the synaptic input to each neuron with the neurotransmitter conductances for each receptor type (g_AMPA, g_NMDA, g_GABAa and g_GABAb): (8)${\displaystyle I=g_{\mathrm{AMPA}}\left(v-0\right)+g_{\mathrm{NMDA}}\frac{{\left[\left(v+80\right)/60\right]}^2}{1+{\left[\left(v+80\right)/60\right]}^2}\left(v-0\right)+\hspace{0.17em}{g}_{\mathrm{GABAa}}\left(v-70\right)+g_{\mathrm{GABAb}}\left(v-90\right).}$ These conductances are modeled with spike timing dependent plasticity (STDP) by being affected by incoming spike’s origin and timing. Therefore, spikes incoming from excitatory neurons can change g_AMPA and g_NMDA, while spikes incoming from inhibitory neurons can change g_GABAa and g_GABAb. This update depends upon the difference between the timing of the post and pre-synaptic spikes, Δt, and two time constants for slow (τ = 0.1 ms) and fast (τ = 0.01 ms) synapses and is modeled exponentially: increment ∝e^−Δt/τ if Δt > 0 and decrement ∝ − e^Δt/τ if Δt ≤ 0.

STP is also modeled in the synaptic weights, influenced by pre-synaptic activity, with the scale factor s(t): (9)${\displaystyle \begin{array}{l}s\left(t\right)=x\left(t\right)u\left(t\right)\hfill \\ \dot{x}=\frac{1-x\left(t\right)}{t_D}-u\left(t\right)x\left(t\right)\delta \left(t-t_{\mathrm{spike}}\right)\hfill \\ \dot{u}=\frac{U-u\left(t\right)}{t_F}-U\left[1-u\left(t\right)\right]\delta \left(t-t_{\mathrm{spike}}\right).\hfill \end{array}}$ Here δ is the Dirac function, and the state variables x(t) and u(t) have baseline levels of 1 and U, respectively. Parameters t_D and t_F are the depression and facilitating times which govern how fast x(t) and u(t) return to baseline. Excitatory synapses’ STP is modeled with U = 0.5, t_D = 800 and t_F = 1,000 while inhibitory synapses have U = 0.2, t_D = 700 and t_F = 20. A small subgroup of excitatory neurons is stimulated by pre-synaptic spikes distributed randomly in time without STP or STDP.

At each column, LPF is obtained assuming dendritic AMPA currents as a good indicator of this activity by the average AMPA conductance in all afferent excitatory synapses (Seth, Chorley & Barnett, 2013).

In this work, we simulated two columns of 5k randomly distributed Izhikevich neurons (80% excitatory, 20% inhibitory) with STDP and STP between neurons of the same column and without both learning strategies in inter-column connections. Five percent of these excitatory neurons are stimulated by pre-synaptic spikes distributed randomly in time without STP or STDP. Instead of a single parameter to model the synaptic strength between columns, like the coupling coefficients or the inter-variable AR coefficient seen in the previous models, here it is possible to model three parameters related to inter-column connection strength. These are the percentage of neurons in one column that project to another column p_projection, the percentage of connections each branch does in the target column p_branch and the synaptic strength of the individual connections k_unit. Simulations produced 60 s of data (the first 20 s were discarded to remove transient effects) with causal influence from column 1 to column 2. In this case, as three parameters define the coupling strength, the simulations were done for each one individually: for varying p_projection (p_projection = [0 to 1 in steps of 0.2], p_branch = 0.05, k_unit = 1, d₂₁ = 20 ms), for varying p_branch (p_projection = 0.2, p_branch = [0.01, 0.05, 0.1, 0.3, 0.6, 1], k_unit = 1, d₂₁ = 20 ms) and for varying k_unit (p_projection = 0.2, p_branch = 0.05, k_unit = [0.01, 0.05, 0.1, 0.3, 0.6, 1], d₂₁ = 20 ms). In conformity with the previous generative models, simulations for different interaction delays were also performed d₂₁ (p_projection = 0.2, p_branch = 0.05, k_unit = 1, d₂₁ = [4 ms and 20–100 ms in steps of 20 ms] ). The spiking units are ate rest in their initial state. Each 60 s of data required ∼52 s of simulation time.

EEG forward modeling

EEG is simulated with a simple forward model that relates the LFP as the activity of a current dipole to the surface potential measured at the scalp. This relation is obtained by representing the head as a multilayer surface, with each layer having constant isotropic conductivity. Although there are realistically shaped multilayer head models like the finite-element method or the boundary element method (Mosher, Leahy & Lewis, 1999; Fuchs et al., 2002; Darvas et al., 2004), these are more important in the EEG inverse problem solution as they mitigate the distortion produced by simpler models (Ermer & Mosher, 2001). As the scope of this work is to simulate a generic EEG signal, we can adopt a simpler model of the skull like the three layer sphere with isotropic conductivities (Berg & Scherg, 1994) where the problem of volume conduction is observed. Figure 2 depicts a typical three layer model with an electrode placed in the scalp with radius r, an intracranial current dipole with radius r_q and moment q and respective angles. This experiment uses a setup with two dipoles placed beneath the three layers with 8 cm outer radius (σ_brain = 0.33 S/m r_brain = 7.04 cm, σ_skull = 0.0042 S/m, r_skull = 7.44 cm, σ_scalp = 0.33 S/m, r_scalp = 8 cm) spaced by 2 cm and two electrodes are placed in the scalp also spaced by 2 cm.

The solution v¹(r; r_q; q) at radius r for the simplest case of a single spherical layer head model for a current dipole with moment q at radius location r_q can be obtained by the sum of the radial $v_r^1\left(\mathbf{r};{\mathbf{r}}_{\mathbf{q}};\mathbf{q}\right)$ and tangential $v_t^1\left(\mathbf{r};{\mathbf{r}}_{\mathbf{q}};\mathbf{q}\right)$ potentials: (10)${\displaystyle \begin{array}{l}{v}_r^1\left(\mathbf{r};{\mathbf{r}}_{\mathbf{q}};\mathbf{q}\right)=\left(\frac{qcos\alpha }{4\pi \sigma }\right)\left(\frac{2\left(rcos\gamma -r_q\right)}{d^3}+\frac{1}{r_{q}d}-\frac{1}{r{r}_q}\right)\hfill \\ v_t^1\left(\mathbf{r};{\mathbf{r}}_{\mathbf{q}};\mathbf{q}\right)=\left(\frac{qsin\alpha }{4\pi \sigma }\right)cos\beta sin\gamma \left(\frac{2r}{d^3}+\frac{d+r}{rd\left(r-r_{q}cos\gamma +d\right)}\right).\hfill \end{array}}$ Here α is the angle the dipole expresses in relation to its location vector r_q, as can be seen in Fig. 2, σ is the conductivity of the shell, d is the value of the direct distance between r_q and r, and β is the angle between the plane formed by r_q and q, and the plane formed by r_q and r.

Following (Berg & Scherg, 1994; Zhang, 1995) a three layer model can be approximated with good accuracy by single layer spheres with the approximation v³(r; r_q; q) ≅ v¹(r; μ₁r_q; λ₁q) + v¹(r; μ₂r_q; λ₂q) + v¹(r; μ₃r_q; λ₃q). The μ and λ are the “Berg parameters” and are used to create three new dipoles with locations consisting in scaling r_q by μ in its radial direction and scaling the moment q by λ. The approximation v³(r; r_q; q) is computed for each electrode and dipole present in the simulation. Noise is added to the resulting time-series.

Details about the computation of the “Berg parameters” can be found in Zhang (1995) and other methods for EEG forward models in Mosher, Leahy & Lewis (1999), Ermer & Mosher (2001) and Darvas et al. (2004). The Brainstorm application offers several methods for the EEG forward model among others (Tadel et al., 2011).

For our simulations, each generated LFP was fed to an EEG forward model with two radial dipoles spaced by 2 cm placed beneath the three layers with 8 cm outer radius (σ_brain = 0.33 S/m r_brain = 7.04 cm, σ_skull = 0.0042 S/m, r_skull = 7.44 cm, σ_scalp = 0.33 S/m, r_scalp = 8 cm) and two recording electrodes placed in the scalp also spaced by 2 cm. White Gaussian noise was added in order for a linear SNR of 10.

BOLD forward modeling

The BOLD signal time-series is the result of a series of neuronal and vascular events that produce a measurable change in the blood hemoglobin concentration. It is therefore an indirect and noisy observation of the neuronal activity as during neuronal activation local vessels are dilated to increase the blood flow and with it, oxygen and glucose delivery. The increased metabolism results in a localized increase in the conversion of oxygenated hemoglobin to deoxygenated hemoglobin and BOLD fMRI uses the latter as the contrast agent (Ogawa & Lee, 1990). This activity can peak four seconds after the neuronal event onset, although this value varies within and between subjects (Handwerker, Ollinger & D’Esposito, 2004). Hence, the location, dynamics and magnitude of the BOLD signal’s activity are vastly influenced by the local vascular bed. This, combined with the fact that fMRI scanners sample an entire volume with TR in the time scale of a second, raised the question if directed functional connectivity, which aims at detecting temporal precedence between neuronal events in the order of tens to hundreds of milliseconds (Ringo et al., 1994; Formisano et al., 2002; De Pasquale et al., 2010), can offer accurate measures from BOLD signals. The simulation of BOLD signals was an important factor to answering this question by allowing experimental control over neuronal and hemodynamic parameters, and this has been achieved mainly by convolution with a canonical hemodynamic response function (HRF) or by dynamic modeling of the vascular activity with the extended Balloon-Windkessel (BW) model (Friston et al., 2000).

The first approach started being used with one gamma function for the HRF convolution kernel (Goebel, 2003; Roebroeck, Formisano & Goebel, 2005) with the purpose of investigating the effect of filtering, down-sampling and noise in GC estimation. Following simulation studies (Deshpande, Sathian & Hu, 2009; Schippers, Renken & Keysers, 2011; Seth, Chorley & Barnett, 2013; Rodrigues & Andrade, 2014) started using a dual-gamma function as used in SPM software (Friston, Holmes & Ashburner, 1999) with parameters as time to peak, time to undershoot, onset time, dispersion of response, dispersion of undershoot and their ratio following the distributions found in Handwerker, Ollinger & D’Esposito (2004) to study the effects of down-sampling, noise and HRF variability. Compared to the convolution approach, the BW model is more biophysically interpretable and can also present nonlinear neuro-vascular couplings although when used in simulation studies there were no considerable differences in GC estimates between both (Smith et al., 2010; Seth, Chorley & Barnett, 2013).

In this work, both approaches were used with the default parameters offered by SPM12 in the functions spm_hdm_priors for the BW model and spm_hrf for the canonical HRF, with a 0.5 TR and a linear SNR of 10. However, in these simulations, new LFPs had to be created for all generative models with different interaction delays, duration and connectivity frequency. Interaction delays are increased to 200 ms, in order to counteract the reduction in sensitivity due to the low sampling period of typical fMRI TRs and low SNR (Deshpande, Sathian & Hu, 2009), and the LFP length had to be increased due to the low sampling rate; hence, 300 s were simulated. For the generative models where it is possible to control the frequency where the causal influence is exerted these were set, from 33 Hz in the previous analysis, to 0.1 Hz in the Kuramoto oscillators and <1 Hz for the MVAR modeling.

Spectral Granger causality

We used the standard Geweke’s spectral decomposition for GC (GGC) (Geweke, 1982) to infer causality between synthetic time-series. Unlike other spectral directed functional connectivity metrics such as PDC or DTF, GGC is not bounded between 0 and 1, which is useful in this study to see how the increase in the connection strength between variables effects the absolute value of causality across the different simulation methods. GGC decomposes the GC index (GCI) (Granger, 1969) into frequency components additively, meaning that the sum of all the frequency components from zero to the Nyquist frequency result in the GCI. For bivariate time-series, which is the case in this study, GGC can be computed from a MVAR model parameters by: (11)${\displaystyle {\mathrm{F}}_{j\to i}\left(f\right)=ln\frac{S_{ii}\left(f\right)}{S_{ii}\left(f\right)-\left({\Sigma }_{jj}-\frac{{\Sigma }_{ij}^2}{{\Sigma }_{ii}}\right)|H_{ij}\left(f\right){|}^2}j,i\in \left\{1,2\right\}.}$ Here Σ is the covariance matrix of the model’s errors, H is the transfer matrix, and S is the spectral matrix. The MVAR model order can be estimated with the Bayesian information criterion (BIC) (Schwarz, 1978) or with the Akaike information criterion (AIC) (Akaike, 1974).

For multivariate time-series, refer to Chen, Bressler & Ding (2006) for more details.

LFP

LFPs represent the local activity of neuronal populations without forward modeling hence, without the confounding effects expected from fMRI BOLD signals or EEG data and, the only randomness is due to the stochastic parameters from each model. This way, this analysis tests each generative model’s capability to produce synthetic datasets with detectable causal effects.

AR modeling

Results for GGC with varying modeled causality and interaction delays can be seen in Figs. 3A and 4A, respectively, and the model order estimated with BIC and AIC for the different interaction delays are shown in Fig. 5A.

The estimated GGC DOI present in Fig. 3A shows, as expected, causal influence at 33 Hz with the same absolute value from what was modeled. Changing the interaction delay from 4 to 100 ms doesn’t affect the estimated causality, as can be seen in Fig. 4A, and both BIC and AIC model orders increase linearly with the increase in interaction delay (Fig. 5A). Also, these model orders correspond exactly to the neuronal delay; with a sampling rate of 250 Hz, the neuronal delays in d₂₁ correspond to [1, 5–25 in steps of 5] lagged observations which equal the model orders suggested by BIC and AIC (except for the neuronal delay of 4 ms).

EEG forward modeling

Figure 6 shows the same experiments as Fig. 3 only this time, the time-series is the EEG recorded at the scalp electrodes. It is possible to see that, except for the Izhikevich columns, GGC amplitude is reduced.

EEG forward modeling does not affect the detection of interaction delays as can be seen in Fig. 7. However, for the KI_e,i sets, causality around the 20 Hz dissipates and causality at 60 Hz is maintained. At some delays (0.4 and 0.8 ms) the EEG from the Kuramoto oscillators (Fig. 7C) shows GGC values similar to the LFPs.

Overall, model orders are overestimated for EEGs from every generative models. In the EEG from MVAR modeling and KI_e,i sets (Figs. 8A and 8B, respectively) this is more pronounced in the lower interaction delays, and in the EEG from the Izhikevich columns (Fig. 8D) the opposite seems to occur. The overestimation problems in the LFPs from the Kuramoto oscillators are more pronounced after EEG forward modeling (Fig. 8C).

BOLD forward modeling

The results for both the convolution with canonical HRF and the extended BW model are shown in Figs. 9 and 10.

The results in Fig. 9 show that GGC is greatly reduced after BOLD forward modeling regardless of the generative or forward models. In BOLD signals generated with MVAR modeling, causal effects are more visible after a modeled causality of 2.5 (Fig. 9A) although these remain around 0.5. The same values are achieved by KI_e,i sets and Kuramoto oscillators after a coupling of 1 and 10 respectively. The Izhikevich columns only achieved these values for p_projection > 0.4, p_branch > 0.6 or k_unit > 0.3. Both the canonical HRF and the extended BW show similar results.

By keeping the coupling strengths at the maximum values in Fig. 9 and varying the interaction delay from 100 ms to 300 ms, the results in Fig. 10 suggest that all generative and forward models benefit from higher delays. With MVAR modeling, 150 ms is the value when causality is detectable (Fig. 10A) whereas for the remaining generative models this value is at 200 ms. Again, there are no relevant differences between canonical HRF and extended BW models.