Synthetic neuronal datasets for benchmarking directed connectivity metrics.

Background: Datasets consisting of synthetic neural data generated with quantifiable and controlled parameters are a valuable asset in the process of testing and validating directed connectivity metrics. Considering the recent debate in the neuroimaging community concerning the use of directed functional connectivity metrics for fMRI data, synthetic datasets that emulate BOLD dynamics have played a central role by supporting claims that argue in favor, or against, certain metrics. Generative models often used in studies that simulate neuronal activity, with the aim of gaining insight into specific brain regions and functions, have different requirements from the generative models for benchmarking datasets. Even though the latter must be realistic, there is a tradeoff between realism and computational demand that needs to be contemplated and simulations that efficiently mimic the real behavior of single neurons or neuronal populations are preferred, instead of more cumbersome and marginally precise ones. Methods: this work explores how simple generative models are able to produce neuronal datasets, for benchmarking purposes, that reflect the simulated effective connectivity and, how these can be used to obtain synthetic recordings of EEG and fMRI BOLD. The generative models covered here are AR processes, neural mass models consisting of linear and non-linear stochastic differential equations and populations with thousands of spiking units. Forward models for EEG consist in the simple three-shell head model while fMRI BOLD is modeled with the Balloon-Windkessel model or by convolution with a hemodynamic response function. Results: the simulated datasets are tested for causality with the original spectral formulation for Granger causality. Modeled effective connectivity can be detected in the generated data for varying connection strengths and interaction delays. Discussion: all generative models produce synthetic neuronal data with detectable causal effects although the relation between modeled and detected causality varies and less biophysically realistic models offer more control in causal relations such as modeled strength and frequency location.


Introduction
The field of brain research that studies connectivity relies in an ever evolving array of methods that aim at finding directed or undirected connectivity links from recorded neuronal datasets which pose different challenges such as having non-linear relations, nonstationary dynamics, low signal-to-noise ratio (SNR), linear mixes or slow dynamics.
Therefore, most studies that introduce or discuss connectivity metrics use not only real but also simulated datasets for demonstration and validity purposes on a controlled environment. These simulated datasets are obtained through generative models that simulate the conception of neuronal activity, like local field potentials (LFPs), with a certain degree of realism and can be followed by forward biophysical modeling where LFPs are transformed into other neuronal recording datasets like electroencephalography (EEG), magnetoencephalography (MEG) or functional magnetic resonance imaging (fMRI) blood activation level dependent (BOLD), for example. With these controlled simulations it is possible to isolate the effect of certain parameters and understand how directed connectivity metrics perform in a realistic range of values.
As far as directed connectivity is concerned, Granger causality (GC) based metrics have been in the center of several studies that discussed their plausibility in the analysis of fMRI BOLD datasets due to the inherent low SNR, slow dynamics compared to neuronal activity and confounding effects due to variable vascular latencies across brain regions (Deshpande, Sathian & Hu, 2009;Smith et al., 2010;Valdes-Sosa et al., 2011;Seth, Chorley & Barnett, 2013;Rodrigues & Andrade, 2014). Most of these studies used simulated datasets to support their claims simulated with generative models such as networks of 'on-off' neurons (Smith et al., 2010), AR processes with real LFPs (Deshpande, Sathian & Hu, 2009;Rodrigues & Andrade, 2014) and columns of thousands of spiking units (Seth, Chorley & Barnett, 2013) and BOLD forward models consisting of linear operations or more biophysically realistic models. Other problems such as the effect of volume conduction in directed connectivity has also been addressed with LFPs being modeled as sinusoids and the forward EEG model as a linear mix of different LFPs (Kaminski & Blinowska, 2014 The objective of this work is to present and test the ability of known generative and forward models to produce synthetic neuronal datasets with modeled causal effects, to be used as benchmarks.

Methods
In this section we succinctly present the methodological framework used in this study from the modeling techniques to the connectivity metrics employed. Figure 1 depicts the generative models to simulate LFPs and the forward models for fMRI BOLD signals and EEG signals. The ensuing subsections expand on each technique separately by laying theoretical foundations, mentioning available software and explaining why these simulations have been important in the development of brain connectivity tools. 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 AR models 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 (dDTF) specific sensitivity to directed causal effects was also demonstrated with data from a multivariate AR 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).
AR 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 in which causality was assessed with GC for different experimental parameters (Roebroeck, Formisano & Goebel, 2005). Later, Schippers used the same AR 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 multivariate AR models to study the effectiveness of GC applied to cluster sets obtained with independent component analysis (ICA) (Londei et al., 2006) and principal component analysis (PCA) with partial canonical correlation analysis (pCCA) (Sato et al., 2010) dependence of GC metrics to experimental parameters (Deshpande, Sathian & Hu, 2009;Rodrigues & Andrade, 2014).
(t X is the multivariate time-series, j Φ is the AR coefficient matrix for time-lag j, p is the model order, D is the interaction delay matrix and ) (t w 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. Considering the AR model for variable 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, 1  p  , the spectral peak frequencies have approximately the same values as the arguments of the roots of the AR polynomial. Other properties, specially 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 multivariate AR 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 2 w  ) 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 multivariate AR matrix, transfer function and noise covariance matrix, and solving it for the inter-variable coefficient.

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 directed transfer function (DTF) can be interpreted within the GC framework (Kamiński et al., 2001) by modeling interacting cortical columns with excitatory and inhibitory populations with delayed-coupled nonlinear stochastic differential equations (SDEs). 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: (1) and inhibitory population defined by: Variables n x and n y represent the LFPs of excitatory and inhibitory populations respectively of node n, ei k and ie k their respective coupling coefficients, nm k is the coupling coefficient from node n to node m, S is a sigmoid function, nm d 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,   (Freeman, 1987), (1) and (2) is a coupling between KOi and KOe subsets hence, each node can be seen a KIe,i set. In this work these parameters were used with the same values as in (Freeman, 1987) and (Kamiński et al., 2001): 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 milliseconds, with the noise term being integrated with the Euler method using the same time-steps.
As (Friston et al., 2014) concludes that GC is not appropriate for data generated by delayed-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 modelled 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. (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: Here, n  is the angular frequency of each oscillator (rad/s), nm C 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 )) ( sin( t r n n   of phase variable. The system's numerical solution was approximated with the Euler-Maruyama (Kloeden & Platen, 1992)

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 time-series sampled with decreasing TRs. Besides spiking units, their synapses were also modeled with explicit NMDA, AMPA, GABAa and GABAb conductances and short-term plasticity. 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 modelled with two ordinary differential equations and four parameters: And the additional after-spike resetting when mV 30 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: excitatory pyramidal . 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 ( 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 Short term plasticity (STP) is also modeled in the synaptic weights influenced by presynaptic activity: Here  is the Dirac function, ) (t s is the scale factor that changes the synaptic weight due to STP, determined by state variables ) (t x and ) (t u with baseline levels of 1 and U respectively. Parameters 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 used modified functions from CARLsim 2.0 (Richert et al., 2011) to simulate 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. Instead of a single parameter to model the synaptic strength between columns, like the coupling coefficients or the AR coefficient seen in the previous models, here it is possible to model three parameters related to intercolumn connection strength. These are the percentage of neurons in one column that project  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 1 depicts a typical three layer model with an electrode placed in the scalp with radius r , an intracranial current dipole with radius q r and moment q and respective angles. This experiment uses a setup with two dipoles placed beneath the three layers with 8 cm outer radius ( Here, α is the angle the dipole expresses in relation to its location vector q r , as can be seen in Figure 1, σ is the conductivity of the shell, d is the value of the direct distance between q r and r , and  is the angle between the plane formed by q r and q , and the plane formed by q r and r . Following (Berg & Scherg, 1994;Zhang, 1995) . The μ and  are the "Berg parameters" and are used to create three new dipoles with locations consisting in scaling q r by μ in its radial direction and scaling the moment q by  . The approximation 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;Darvas et al., 2004). The Brainstorm application offers several methods for the EEG forward model among others (Tadel et al., 2011).

BOLD forward modeling
The BOLD 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 (HbO) to deoxygenated hemoglobin (HbR) and BOLD fMRI uses the latter as the contrast agent (Ogawa & Lee, 1990). This activity can peak four seconds after the neuronal event onset though this value varies within and between subjects (Handwerker, Ollinger & D'Esposito, 2004). Hence, the location, dynamics and magnitude of the BOLD activity are vastly influenced by the local vascular bed. This, combined with the fact that fMRI scanners sample an entire volume with time repeat (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 important to answer 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)

Spectral Granger causality
We use 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 an AR model parameters by: Here  is the covariance matrix of the model's errors, H is the transfer matrix, and S is the spectral matrix. The AR 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.

Results and discussion
This section presents the results of applying (8)  Simulations ran on a PC equipped with an Intel® Core™ i7-2600K CPU @3.4GHz, 8 GB of RAM and an NVIDIA® GeForce® GTX 580 GPU with 3 GB graphics memory.

LFP
LFPs represent the local activity of neuronal populations without forward modeling hence, without the confounding effects expected from fMRI BOLD 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
A bivariate AR(2) model was created with variable 1 exerting causal effect in variable 2 at a specific frequency  , with variable intensity 2 1 F X X  and delay 21 d following the expression: Here 1 w and 2 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 in the 33 Hz so, the AR parameters for variable 1 AR process were chosen so the spectral peak occurs at this frequency. Following the relationship between spectral peak frequency and AR parameters found in (Jiru, 2008) for AR (2) Figure 5A. The simulations produced 60 s of data and the first 20 s were discarded to avoid initial conditions, data was generated at 1k Hz and subsampled to 250 Hz. Each 60 s of data required ~1 s of simulation time.
The estimated GGC DOI present in Figure 3A shows, as expected, causal influence on the 33 Hz with the same absolute value from what was modeled. Changing the interaction delay from 4 to 100 milliseconds doesn't affect the estimated causality, as can be seen in Figure 4A, and both BIC and AIC increase linearly with the increase in interaction delay ( Figure 5A). Also, these model orders correspond exactly to the neuronal delay: with a  (1) and (2) Figure 3B shows that there are two frequency bands where connectivity is detected, around 20 Hz and 60 Hz. These causal influences are more expressive for strong couplings (>1) and are inexistent at weak couplings below 0.5. When the interaction delay changes the connectivity also changes frequency ( Figure 4B). This occurs as different delays change the oscillatory behavior of these neuronal populations and also because the AR model order also changes which influences the AR spectrum. In Figure 5B it is possible to see that both BIC and AIC are sensitive to the increase in the interaction delay and overestimate the true model order by approximately 2 lagged observations.

Phase-delayed coupled Kuramoto oscillators
A network of two nodes was simulated 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 (3). 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 to 2. GGC only detects causal relations with positive DOI in the 33 Hz as can be seen in Figure   3C. However, these only present expressive values on strong couplings (>3). Different interaction delays do not affect the GGC intensity or frequency distribution ( Figure 4C) although the model order suggested by BIC and AIC is not linearly related to the interaction delay ( Figure 5C) and is overestimated. can be seen in Figure 3D, results for  Figure 4D and the respective model orders estimated with BIC and AIC is shown in Figure 5D.

Izhikevich columns
From Figures 3D, 3E, 3F it is possible to see that GGC detects connectivity below 20 Hz and that a linear increase only occurs with unit k . Interaction delays change the connectivity frequency distribution ( Figure 4D) as different delays change the oscillatory behavior of the neuronal populations and different AR model order produce different AR spectrums.
Increasing the interaction delay produces a linear increase in model orders estimated with BIC although these are slightly overestimated ( Figure 5D).    Kuramoto oscillators, D) Izhikevich columns. The lagged observations for these interaction delays with 250 HZ sampling rate are [1,5,10,15,20,25].

EEG forward modeling
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 brain r = 7.04 cm, skull  = 0.0042 S/m, skull r = 7.44 cm, scalp  = 0.33 S/m, scalp r = 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. Figure 6 shows the same experiments as Figure 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.    ( Figure 8D) the opposite seems to occur. The overestimation problems in the LFPs from the Kuramoto oscillators are more pronounced after EEG forward modeling ( Figure 8C). observations for these interaction delays with 250 HZ sampling rate are [1,5,10,15,20,25].

BOLD forward modeling
For the BOLD simulations new LFPs were created with all generative models differing from the previous in the interaction delays, length and connectivity frequency. Interaction delays are increased to 200 milliseconds 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 The results in Figure 9 show that GGC is greatly reduced after BOLD forward modeling regardless of the generative or forward models. In BOLD generated with AR modeling causal effects are more visible after a modeled causality of 2.5 ( Figure 9A)   By keeping the coupling strengths at the maximum values in Figure 9 and varying the interaction delay from 100 milliseconds to 300 milliseconds the results in Figure 10 suggest that all generative and forward models benefit from higher delays. With AR modeling, 150 milliseconds is the value when causality is detectable ( Figure 10A) whereas for the remaining generative models this value is at 200 milliseconds. Again, there are no relevant differences between canonical HRF and extended BW models.

Discussion
This study explores four generative models that represent distinct methodologies: multivariate AR modeling, neural mass models and spiking neuron populations. Although most of these models have been used previously in simulation studies that aim at benchmarking connectivity metrics, their capability to reflect causal interactions in their generated neuronal time-series had not been compared yet. Also, these neuronal time-series are not limited to the LFPs, BOLD and EEG covered by this work but these are the most widely used in connectivity studies. these subgroups of generative models. AR modeling and Kuramoto oscillators, as these allow to manipulate the strength and frequency of causal relationships (or even the theoretical GGC in the AR modeling), seem more adequate for initial testing of directed functional connectivity metrics (Baccalá & Sameshima, 2001), studying the effect of post processing or other transformations in data prior to causality inference (Barnett & Seth, 2011) or to compare different metrics performance in an extended benchmark (Rodrigues & Andrade, 2014). Due to their superior neurophysiological realism KIe,i sets and Izhikevich columns are more useful to inquire about the effectiveness of certain connectivity metrics for data recorded in specific brain locations with known dynamics (Friston et al., 2014). For example, in (Freeman, 1987) coupled KI sets are used to simulate chaotic EEG emanating from the olfactory system and in (Richert et al., 2011) Izhikevich columns are used to simulate a large-scale model of cortical areas V1, V4, and middle temporal (MT) with color, orientation and motion selectivity.
All generative models are able to produce LFPs with detectable causal relations however, and these results show the relation between experimental parameters and simulated causality. Only AR modeling allows for a direct specification of causality by solving the analytical equations for GGC applied to the AR coefficients. Concerning the neural mass models, the KIe,i sets show causal effects with lower coupling strength than the Kuramoto oscillators as the first start having causal relations with 21 k = 0.7 while the second required values of k > 3 for identifiable causality. In the Izhikevich columns the percentage of projecting neurons from the source column projection p is the variable with least influence in the observed GGC, except for when it is zero. Increases in both the percentage of target connections per projection branch p and the synaptic strength unit k lead to increases in the observed GGC. Although it wasn't tested in this work it is possible that these values change for different number of neurons per column. In all generative models the modulation of the interaction delay is possible without loss of causal relations although the neurophysiologically realistic models show different frequency spectrums for different interaction delays. This neuronal delay is also detected by the BIC and AIC which suggest higher model orders for higher interaction delays. Forward models lead to a decrease in the absolute value of GGC, specially the forward BOLD model where negative DOIs could be found in worst scenarios. EEG forward modeling reduces the estimated GGC due to the added noise and to the volume conduction effect that produces a "cross-talk" between the two neural populations and electrodes and these effects would be more adverse if more neuronal populations were added inside the three-shell sphere. Nevertheless, these results show that all generative models produce LFPs that, when treated as activity from radial dipoles, preserve the causal relations in the resulting EEG. BOLD forward modeling is more detrimental to the preservation of causal effects and this work did not model hemodynamic variability between brain regions (Handwerker, Ollinger & D'Esposito, 2004) which would further affect causality preservation in BOLD time-series (Deshpande, Sathian & Hu, 2009). Causality was reduced approximately by a factor of 10 in all generative models and there were no concise differences between HRF convolution and BW modeling except for the Kuramoto models where the last leads to smaller values of causality.

Conclusion
This work presented and analyzed different modeling strategies to generate artificial neuronal datasets for benchmarking purposes. LFPs are obtained by generative models and can be used by forward models to produce other recordings of neuronal activity such as the BOLD signal or EEG. All the analyzed models were able to transmit their causal structure ( ) 1 ( 21  , 21 k , k , projection p , branch p , unit k ) into their generated data however with different relations between these and the identified GGC. This study only covered bivariate models but the same analysis could be performed with larger networks with large-scale fluctuations (Cabral et al., 2011). This would be useful to identify the directed functional connectivity metrics most appropriate to analyze large scale data such as fMRI BOLD from resting state networks.