1 Introduction

Communication between cortical neuronal networks is fundamental for cognitive (Campo et al. 2015) and sensorimotor (Brovelli et al. 2004; Seidler et al. 2015) functions. A current challenge in neuroscience is how to define the pathways of neuronal communication based on synaptic connections and information flow between areas. The possibility of recording the local field potentials (LFPs) and the activity of hundreds of individual neurons enables the use of statistical methods to infer these communication patterns. Brain connectivity can be described as: (i) structural, referring to anatomical connections between pairs of neurons or brain regions; (ii) functional, denoting the characterization of temporal correlations between electrical signals from brain regions (Goñi et al. 2014); and (iii) effective, corresponding to the inference of information flow direction or causal relationships between recorded brain signals (Rubinov and Sporns 2010; Bullmore and Sporns 2009; Friston 2011).

Granger causality (GC) (Granger 1969) analysis is a well-established method for effective connectivity inference that has been extensively applied in neuroscience (Bernasconi and König 1999; Zhang et al. 2012). A time series X is considered to Granger-cause a time series Y if past information from X helps to predict the future of Y better than when using past information only from Y (Barnett and Seth 2014). The original GC is a time domain method of effective connectivity and does not provide frequency information about the interactions between brain regions. Some frequency domain counterparts of GC are spectral GC (Geweke 1982), directed transfer function (DTF) (Kaminski and Blinowska 1991), and partial directed coherence (PDC) (Baccalá and Sameshima 2001). DTF and spectral GC detect both direct and indirect causal relationships between time series. PDC has the advantage of detecting only direct connections (Baccalá and Sameshima 2001), but its values depend on the signal amplitude. Generalized partial directed coherence (GPDC) (Baccalá et al. 2007) solves this problem by changing the normalization factor from PDC. GPDC has been applied to evaluate information flow from EEG signals (Omidvarnia et al. 2014; Hu et al. 2017), to identify the directional coupling in LFPs from macaque V4 area (Hoerzer et al. 2010) and to verify the effective connectivity from fMRI signals (Sato et al. 2009; Ning et al. 2018).

Several studies evaluated connectivity methods using electrophysiological data (Omidvarnia et al. 2014; Cadotte et al. 2010; Youssofzadeh et al. 2016; Gao et al. 2015; Sato et al. 2009) and autoregressive models (Baccalá and Sameshima 2001; Faes and Nollo 2010; Seth 2010; Papana et al. 2013; Sommariva et al. 2017). On the one hand, gauging the method using electrophysiological data poses a challenge since the actual structural connectivity between areas is unknown. On the other hand, autoregressive models assume data to be linear and stationary, conditions that may be violated when dealing with brain signals.

In fact, very little is known about the performance of the GPDC when it is applied to infer effective connectivity between brain signals. For example, how does the number of false positives depend on the coupling strengths or the number of physical connections? This question is central for the viability of these methods of real brain data analysis.

A possible solution would be using simulations of spiking neuronal networks based on the cortical architecture. They can provide simulated LFP-like signals, permitting the evaluation of these methods when they are applied to signals with characteristics of electrophysiological data. Moreover, since we can control the network behavior and connectivity pattern, we can relate the causal connections inferred by the methods with connectivity characteristics, such as the presence of physical connections and connection strengths.

In this work, we evaluate the use of GPDC to infer effective connections between networks in models of spiking neuronal networks from synthetic LFP signals, using several connectivity configurations. We also analyze the rate of true and false positives for different coupling strengths and different levels of noise in the LFP signals. Finally, we evaluate the relationship between GPDC values and coupling strength.

2 Methods

2.1 Neuronal model

We created the neuronal network model using Izhikevich neurons (Izhikevich 2003), where a neuron i is represented by equations:

$$\begin{aligned}&\displaystyle \dot{v_i}= 0.04v_i^2+5v_i+140-u_i+\sum _{k~\in ~ \mathrm{{presyn}}}I_{i,k}(t), \end{aligned}$$
(1)
$$\begin{aligned}&\displaystyle \dot{u_i} = a(bv_i-u_i), \end{aligned}$$
(2)
$$\begin{aligned}&\displaystyle I_{i,k}(t) = g_{i,k}(t)(E_{k}-v_{i}), \end{aligned}$$
(3)

with the spiking condition that resets \(v(t)\rightarrow c\) and \(u(t)\rightarrow u(t)+d\) when v(t) reaches 30 mV. We modeled excitatory cells as regular spiking (RS) and inhibitory cells as fast spiking (FS). For each neuron i, we used the parameters \((a,b,c,d)=(0.018,0.18,-\,65.2,7.8)+(0.004,0.4,0.4,0.4)\cdot r_i\) for RS cells and \((a,b,c,d)=(0.08,0.18,-\,65.2,1.8)+(0.04,0.4,0.4,0.4)\cdot r_i\) for FS cells, where \(r_i\) is randomly drawn from a uniform distribution from 0 to 1 (Izhikevich 2003; Tomov et al. 2014). \(I_{i,k}\) is the input current from each presynaptic neuron k, \(E_{k}\) the reversal potential (0 mV for excitatory and \(-80\) mV for inhibitory connections), and \(g_{i,k}\) the synaptic conductance, given by:

$$\begin{aligned} g_{i,k}=\frac{g_\mathrm{{max}}}{G_{\mathrm{{peak}}}}G(t) \end{aligned}$$
(4)

where,

$$\begin{aligned} \ddot{G}_{i}+\alpha \dot{G}_{i}+\beta G_{i}=\overline{g}_{k}\sum _{j}x_{ij}(t), \end{aligned}$$
(5)

and

$$\begin{aligned} \alpha =\frac{\tau _{1}+\tau _{2}}{\tau _{1}\tau _{2}}, \qquad \beta =\frac{1}{\tau _{1}\tau _{2}}. \end{aligned}$$
(6)

G(t) is the synaptic conductance waveform, \(G_\mathrm{{peak}}\) the maximum value of G(t), \(g_\mathrm{{max}} = 1\,\upmu \mathrm{{S}}\) the maximum synaptic conductance value, and \(\tau _{1}\) and \(\tau _{2}\) are the decay and rise time constants, respectively. When \(\tau _{1}=\tau _{2}\), G(t) takes the form of an alpha function, otherwise it has the dual exponential form (Koch and Segev 1988). \(x_{ij}(t)\) represents the arriving spikes at neuron i from each neuron j and behaves like a delta function \(\delta (t-t_{q})\), where \(t_q\) is the arrival time of spike q (Sterratt et al. 2011). The dimensionless conductance increment \(\overline{g}_{k}\) is the synaptic weight and depends on the type of the presynaptic neuron k (excitatory or inhibitory) and is divided by the integration time step for normalization. These equations represent channels as dynamical systems which are perturbed by synaptic inputs, permitting the efficient computation of the synaptic conductances (Koch and Segev 1988). The advantage of this method is that the influence of each spike arriving at the synapse on the postsynaptic neuron needs to be considered only at its arrival time. Table 1 shows the synaptic weights and time constants for each synaptic type. We integrated the equations using Euler method and a time step of 0.05 ms and used Matlab to implement and simulate the model.

Our simulations comprise five networks, each representing local circuitry from different regions of the cortex, with excitatory and inhibitory cells. These local networks are connected among themselves by the so called long-range connections.

In the local network, we simulated connections from excitatory to excitatory neurons, from excitatory to inhibitory neurons, from inhibitory to excitatory neurons and from inhibitory to inhibitory neurons. Each neuron, excitatory or inhibitory, receives an independent random external input, modeled as a 10 Hz independent Poisson spike generator. This input simulates random fluctuations of the membrane potential due to excitatory synapses arriving at each neuron from connections that were not specifically simulated.

Table 1 Synaptic conductance increments and time constants (Eqs. 5 and 6)
Full size table

Synapses among neurons within the same network have random propagation delays extracted from a uniform distribution in the range \({D}_{E\rightarrow E}=[1-10\,{\mathrm{ms}}]\) and \({D}_{E\rightarrow I}=[1--5\,{\mathrm{ms}}] \), and inhibitory connections have fixed delays \({D}_{I\rightarrow E}={D}_{I\rightarrow I}=1\,{\mathrm{ms}}\), where E and I represent excitatory and inhibitory neurons, respectively, and the arrow indicates the direction of the connection. Long-range connections have delays of \(15~{\mathrm{ms}}\) (Izhikevich 2006).

The complete network model consisted of five local networks with 800 excitatory and 200 inhibitory neurons each. Neurons within local networks were randomly connected, with a fixed indegree for each type of the connection (Akam and Kullmann 2010). Table 2 shows the number of connections from presynaptic neurons from the same network for each of the five networks. We choose a different number of connections for each network to create different network behaviors.

Connections between networks were modeled by long-range connections, which were always excitatory (Sancristóbal et al. 2014) and could target both excitatory and inhibitory neurons. Unless otherwise noted, each neuron received 10 connections with synaptic weight 0.015 from random neurons from each source network.

Table 2 Number of presynaptic neurons for each postsynaptic neuron in the networks created
Full size table

2.2 LFP signal

Local field potentials (LFP) are often used to record local cortical activity and are obtained after low-pass filtering extracellular electrical potential signals captured by implanted electrodes (Lindén et al. 2011). Pyramidal neurons contribute predominantly to LFP signals due to the relative open-field geometrical arrangement of their long apical dendrites, with excitatory and inhibitory currents substantially spaced along these dendrites (Mazzoni et al. 2011; Cavallari et al. 2014; Buzsáki et al. 2012). There are several ways to simulate LFP signals in networks of phenomenological neuronal models (Mazzoni et al. 2015). We approximated the LFP signals as the resultant from current dipoles originated from synaptic currents flowing into/from pyramidal neurons. We used the sum of absolute values of AMPA and GABA currents acting on excitatory neurons (Sancristóbal et al. 2014; Mazzoni et al. 2008):

$$\begin{aligned} {\mathrm{LFP}}=\frac{R_\mathrm{{e}}}{N_\mathrm{{E}}}\left( \sum _{i}^{N_\mathrm{{E}}}{\left| I_{i,\mathrm{{AMPA}}}(t)\right| } +\sum _{i}^{N_\mathrm{{E}}}{\left| I_{i,\mathrm{{GABA}}}(t)\right| }\right) , \end{aligned}$$
(7)

where \(N_\mathrm{{E}}\) is the number of excitatory neurons in the network and \(R_\mathrm{{e}}\) denotes the resistance of an electrode used for recording extracellular activity, defined as \(1~M\varOmega \). The terms \(\left| I_{i,\mathrm{{AMPA}}}(t)\right| \) and \(\left| I_{i,\mathrm{{GABA}}}(t)\right| \) are the total AMPA excitatory and GABA inhibitory currents, respectively. The sample average of the signal generated was subtracted from the original signal generating a zero-mean time series. The resultant signal was low-pass filtered at 200 Hz to mitigate aliasing effects caused by the downsampling used to lower the sample rate from 20 kHz to 200 Hz. Finally, we applied a linear detrending to the signal (Matias et al. 2014).

To compute the power spectral density (PSD), we used the multitaper method (Thomson 1982), implemented by Matlab function pmtm. In this method, the PSD is the average of \(K=7\) modified power spectra obtained through orthogonal Slepian taper functions as the window. We considered samples of LFP with 1000 data points, representing 5 seconds of activity in the neuronal network.

2.3 Partial directed coherence

Partial directed coherence (PDC), introduced by Baccalá and Sameshima (Baccalá and Sameshima 2001), is a multivariate frequency domain measure of the directed relationship between pairs of time series. PDC is not scaling invariant and changes in time series amplitude can lead to substantial changes in PDC values (Baccalá and Sameshima 2001). To overcome this deficiency, the generalized PDC (GPDC) was introduced (Baccalá et al. 2007).

Let \({\mathbf {x}}(n)=[x_{1}(n)\cdots x_{N}(n)]^{T}\) be a set of simultaneously acquired time series. The vector autoregressive (VAR) model for \({\mathbf {x}}(n)\) is defined as:

$$\begin{aligned} {\mathbf {x}}(n)=\sum _{k=1}^{p}{\mathbf {A}}_{k}{\mathbf {x}}(n-k)+{\mathbf {w}}(n) \end{aligned}$$
(8)

where p is the VAR model order. \({\mathbf {A}}_k\) are coefficient matrices, where element \(a_{ij}^{(k)}\) describes the influence of \(x_{j}(n-k)\) on \(x_{i}(n)\). \({\mathbf {w}}(n)\) is a zero-mean innovation process composed of white uncorrelated noises with covariance matrix \({\varvec{{\Sigma }}}\). The choice of p is set according to Akaike’s information criterion (AIC) (Akaike 1974), given by

$$\begin{aligned} {\mathrm{AIC}}(p)=n \log (\det (\varSigma ))+2pN^2, \end{aligned}$$
(9)

where n is the number of data points and N is the number of time series (Kamiński and Liang 2005). The GPDC from the time series \(x_j\) to the time series \(x_i\) at frequency \(\lambda \) is defined as,

$$\begin{aligned} \overline{\pi }_{ij}(\lambda )\triangleq \frac{\frac{1}{\sigma _{i}}\overline{A}_{ij}(\lambda )}{\sqrt{\sum _{k=1}^{N}\frac{1}{\sigma _{k}^2} \overline{A}_{kj}(\lambda )\overline{A}_{kj}^{*}(\lambda )}}, \end{aligned}$$
(10)

where

$$\begin{aligned} \overline{A}_{ij}(\lambda ) = {\left\{ \begin{array}{ll} 1- \mathop {\sum }\nolimits _{k=1}^{p}A_{ij,k}e^{-2\pi \lambda k {\mathbf {i}}}, &{} \text{ if } \quad i=j \\ - \mathop {\sum }\nolimits _{k=1}^{p}A_{ij,k}e^{-2\pi \lambda k {\mathbf {i}}}, &{} \text{ if } \quad i\ne j, \end{array}\right. } \end{aligned}$$
(11)

where \({-2\pi \lambda k {\mathbf {i}}}\) is a complex number and \(\sigma _{i}^{2}\) refers to the variance of the innovation process \(w_{i}(n)\) (Baccalá et al. 2007).

The VAR model was estimated by the method of ordinary least squares (OLS) (Hamilton 1994). AIC indicated that the best model order p that was less than or equal to 10.

2.4 Statistical significance

All simulations used five networks, generating five LFP signals. We executed 10 trials for each simulated network. We computed the effective connectivity between neuronal networks using the average GPDC over trials, for all frequencies.

To obtain a threshold for statistical significance of the average GPDC, we used a nonparametric bootstrap algorithm. With the VAR models for the LFP signals for each network, we resampled the residuals. To test the influence from network j to network i, we made the coefficients \(A_{i,j,k}\) from the VAR models equal to zero, for all values of k, while keeping the other coefficients unchanged. With the resampled residuals and the coefficients, we simulated a bootstrapped time series under null hypothesis of no causality from network j to network i. To generate a bootstrap sample, we repeated this process for each trial and computed their average GPDC. We generated 10000 bootstrap samples for each evaluated connection (Sato et al. 2009).

The critical value of the GPDC was defined as the \((1-\alpha )\) quantile of the bootstrap samples, where we used \(\alpha =0.05\) as the significance level adjusted by Bonferroni correction (Maris and Oostenveld 2007). In Figs. 2c, 3c, and 4c, the correction was made for 100 frequencies. In Figs. 5 and 6, we corrected for 25 frequencies. We considered the GPDC values as statistically significant when their average values minus the standard error of the mean was higher than the threshold of statistical significance.

2.5 Signal-to-noise ratio (SNR)

We evaluated the performance of GPDC to detect the connectivity when different levels of noise were applied to the synthetic LFP signals. We used the Matlab function awgn to add white Gaussian noise to the LFP signals according to a desired signal-to-noise ratio (SNR) (Shakil et al. 2016). The SNR in decibels is given by,

$$\begin{aligned} {\mathrm{SNR}}=10\log _{10}\left( \frac{P_\mathrm{{SN}}}{P_\mathrm{{N}}}\right) \end{aligned}$$
(12)

where \(P_\mathrm{{SN}}\) denotes the power of the signal with noise and \(P_\mathrm{{N}}\) denotes the power of the white Gaussian noise (Lowet et al. 2016). The power of the signals is defined as,

$$\begin{aligned} P=\sum _{i}^{n}\frac{|S_{i}|^{2}}{n} \end{aligned}$$
(13)

\(S_{i}\) is each i data-point of the signal and n is the length of the signal.

2.6 ROC curve

We generated a ROC curve (Ito et al. 2011a; Garofalo et al. 2009) by plotting the true positives rate (TPR) versus the false positive rate (FPR) to different values of the significance level\(\alpha \). TPR represents the rate of correct detection, while FPR represents the rate of false alarms. These two rates are given by:

$$\begin{aligned} {\mathrm{TPR}}= & {} \frac{{\mathrm{TP}}}{{\mathrm{TP}}+{\mathrm{FN}}}, \end{aligned}$$
(14)
$$\begin{aligned} {\mathrm{FPR}}= & {} \frac{{\mathrm{FP}}}{{\mathrm{FP}}+{\mathrm{TN}}}. \end{aligned}$$
(15)

A connection identified by GPDC can be a true positive (TP) if it represents a synaptic connection among the neuronal networks, or a false positive (FP) if the connection does not exist. Similarly, connections not detected by the GPDC can be true negative (TN), when the connection in the simulation does not exist, and false negative (FN) if the connection actually exists. We considered that GPDC inferred a connection when the mean of GPDC subtracted by the standard error of the mean is higher of the bootstrap threshold for some frequencies. We performed 100 simulated experiments each one composed of five trials. The significance levels ranged from 0 to 1 in steps of 0.01, adjusted by Bonferroni correction.

3 Results

We considered three models, each composed of five local neuronal networks, with different connectivity patterns between them (Figs 2a, 3a, 4a). Long-range connections between networks were excitatory and could target either inhibitory (open circles) or excitatory neurons (filled circles) and their overall effect on the postsynaptic network was inhibitory and excitatory, respectively. Hence, for simplicity, we just used the short-hand terms “excitatory” and “inhibitory” to characterize the networks and their projections.

The intrinsic dynamics of each local neuronal network without long-range connections exhibit relatively uncorrelated spikes and some transient periods of oscillations in the LFP signals (Fig. 1).

Fig. 1
figure 1

Raster plot and LFP signals for 1 s of simulation. a Raster plot for the five standalone networks used. Neurons 1 to 800 are excitatory and 801–1000 inhibitory. b LFP signals after demeaning, detrending, and low-pass filtering and downsampling to 200 Hz

Full size image

3.1 Inference of effective connectivity

The first model (Fig. 2a) simulates a scenario that includes two simple connections, one excitatory from network 1 to network 2 and another inhibitory, from 1 to 3. It also includes a reciprocal connection pattern, with network 4 inhibiting network 5 and network 5 exciting network 4. Finally, network 4 is also inhibited by network 1, permitting the evaluation of the interactions between the two connections arriving at network 4.

The isolated neuronal networks possess different power spectral densities (PSD), except for networks 3 and 4, which have approximately the same spectrum of frequencies (Fig. 2b). Enabling long-range connections between networks result in small changes in PSD, more notably for excitatory interactions, such as the increase at 20Hz in network 2, which is the peak power from network 1, and the increase at 35Hz of network 4, probably caused by network 5. The effects of inhibitory interactions on firing rates are also present (Table 3), but are less evident than for the excitatory ones. The interactions among networks caused limited changes in their behavior, which is the ideal scenario for testing the effectiveness of causality methods.

The GPDC method identified all interactions between networks. The first column in Fig. 2c (\(j=1\)) shows that the method recognized connections from network 1 to networks 2, 3, 4 (solid lines). Other connections from network 4 to 5 and from network 5 to 4 were also correctly identified. The GPDC value for all other connections was below the significance curve (dashed line), resulting in no false positives. More importantly, it accurately did not detect the indirect connection \(1 \rightarrow 4 \rightarrow 5\), marking the connection from network 1 to 5 as non-existing. The spectral coherence (gray lines) also identified all interactions, but their bidirectional nature provides no information on the direction of the interaction. Also, indirect connections could not be discriminated from direct ones using spectral coherence.

Fig. 2
figure 2

Model 1. a Directed graph representing the connectivity model between neuronal networks. Each node represents a network identified according to Table 2. Edges describe excitatory long-range connections targeting inhibitory (open circle) and excitatory (closed circle) neurons. b Average power spectral density (PSD) between trials with long-range connections disabled (uncoupled) and enabled (coupled). Shaded gray areas represents the standard error of the mean (SEM). The numbers inside the plots indicate the corresponding neuronal networks. c Average GPDC (black solid lines) and average coherence between trials (gray solid lines) for all the connections between networks, where j presents the source of connection and i the target. Shaded areas represent the SEM of the GPDC and dotted lines the threshold of statistical significance

Full size image

The second model (Fig. 3a) contains excitatory connections from networks 2 to 3, 5 to 4, and 5 to 1, and inhibitory connections from networks 1 to 2, 3 to 4, and 4 to 5. This model has two loops: a long one, encompassing all networks, and a smaller one, involving only networks 4 and 5. Similarly to the first model, the PSD of the networks shows small changes when coupled (Fig. 3b) and changes in the average firing rate (Table 3) are in accordance with the connection type. Using GPDC, we could infer all the connections of the model, with no false positives, including the bidirectional connections between 4 and 5 (Fig. 3c). Spectral coherence between interacting networks also detect the connections, showing a peak in the frequency near the maximum GPDC value, but without information on connection direction.

Finally, the third model (Fig. 4a) contains two pathways from network 1 to 3, one direct, and another indirect, passing through network 2. Connections from network 1 to 2, 2 to 3, 3 to 4, and 4 to 5 are inhibitory, and connections from 1 to 3 and 4 to 5 are excitatory. Differently from the previously analyzed models, there are two consecutive inhibitory interactions, from network 1 to 2 and 2 to 3. The PSDs again showed only small variations after enabling the long-range connections (Fig. 4b) and the average firing rate also changed according to connection type. The GPDC method identified correctly all connections established in the model and did not find false positives (Fig. 4c). Although the GPDC values are smaller, they are above the significance threshold for a range of frequency values. As the previous models, the spectral coherence presents a peak in the same frequency of the peak of GPDC.

Table 3 Average firing rate over trials
Full size table

3.2 Effects of coupling strength

We investigated the effects of coupling strength, characterized by synaptic weights and number of long-range connections, on GPDC, using network model 3 (Fig. 4a).

Fig. 3
figure 3

Same as Fig. 2, applied to model 2

Full size image

We compared the true positive rate (TPR) and false positive rate (FPR) for different synaptic strengths, number of connections, and levels of noise applied to the LFP signal. The area under the ROC curve (AUC) is larger for higher synaptic weights (Fig. 5a) and number of long-range connections (Fig. 5b). With 10 long-range synapses per target neuron and a synaptic weight of 0.015, the AUC is almost 1, which means that causal relationships were always detected with no false positive. With smaller weights, the detection is less reliable, with lower TPR and higher FPR values, until reducing toward chance detection (Fig. 5a). Similarly, when using fewer afferent connections per neuron, the AUC is reduced almost to chance level (Fig. 5b). Consequently, GPDC can reliably detect causal relationships with stronger network couplings, but it is less reliable as the coupling strength decreases.

The detection rate of GPDC also changed when different levels of white Gaussian noise are applied to the LFP signals (Fig. 6). Increasing noise level caused a decrease in the connection detection rate. However, even in the case where the SNR is negative, i.e., the power of the noise is higher than the power of the LFP with noise, the detection rate was still above chance level.

Fig. 4
figure 4

Same as Fig. 2, applied to model 3

Full size image
Fig. 5
figure 5

ROC curves from model 3. a Thick curves represent different synaptic weights, considering 10 long-range synapses per target neuron. b Thick curves represent different numbers of long-range synapses per neuron, with synaptic weights of \(1,5\times 10^{-2}\). Each point represents a difference significance threshold. The dotted diagonal thin line in both plots represents the chance level. We considered 100 GPDC values, where each is the average over 5 trials. TPR: true positive rate. FPR: false positive rate

Full size image
Fig. 6
figure 6

ROC curves from model 3. Thick curves represent different SNR values in decibels. We considered 10 long-range synapses per target neuron with synaptic weights of \(1,5\,\times \,10^{-2}\). Each point represents a difference significance threshold. The dotted diagonal thin line represents the chance level. We considered 100 GPDC values, where each is the average over 5 trials of the GPDC values. TPR: true positive rate. FPR: false positive rate

Full size image

We also investigated the relationship between GPDC values and coupling strength. We found a linear relationship between the maximum GPDC values and the number of long-range synapses that arrived to each neuron in the target network (Fig. 7a), with coefficients of determination (\(R^2\)) between 0.64 and 0.77. We found a similar relationship between GPDC values and long-range synaptic weights (Fig. 7b), with \(R^2\) between 0.58 and 0.75. In both cases, coupled networks resulted in p values below \(10^{- 3}\). For uncoupled networks, \(R^2\) was always below 0.02. These results indicate that the GPDC value could be used to estimate network coupling strength.

Fig. 7
figure 7

a GPDC versus the number of long-range synapses per neuron (2,4,6,8, and 10) and b GPDC versus long-range synapse weights (5,8,10,12, and 15 (\(\times \,10^{-3}\))). Each point is the average over 5 trials of the maximum GPDC values. We used 100 GPDC values for each coupling strength and pair of networks from model 3

Full size image

4 Discussion

In this study, we evaluated the GPDC as a technique to infer effective connections between neuronal networks, using simulated LFP signals from three models containing coupled networks of Izhikevich spiking neurons (Izhikevich 2003). The objective was to verify if GPDC could correctly infer connections from signals that are not linear autoregressive time series and have characteristics similar to electrophysiological signals.

We used network models that try to capture general properties of the cortex, such as the distribution of excitatory and inhibitory neurons, but without matching them to any specific cortical area, with the use of random connections. We considered several scenarios with both inhibitory and excitatory interactions between networks, including one source for several targets, multiples sources for a single target, recurrent interactions, indirect interactions, among others. Our results showed that in every model evaluated, the average GPDC was statistically significant only when there was a connection between the networks. Moreover, the PSD of the networks and average firing rates showed small differences when the networks were coupled and uncoupled, suggesting that the coupling among networks caused small changes in network behavior, but that was large enough to be detected by the GPDC.

Different from PDC, GPDC is scale invariant, meaning that it yields values that are independent of the dynamic range of the time series analyzed (Baccalá et al. 2007). But some disadvantages of PDC remain in GPDC. These two methods are normalized with respect to the number of target networks (Baccalá and Sameshima 2001) so adding more target networks influenced by the same source decreases the magnitude of both PDC and GPDC  (Schelter et al. 2009). This seems to limit the use of the GPDC values to infer the coupling strength. Some methods were proposed to overcome these limitations such as a renormalized PDC (Schelter et al. 2009) and isolated effective coherence (iCoh) (Pascual-Marqui et al. 2014). Renormalized PDC introduced a normalization using the variance of the influences between the signals (Pascual-Marqui et al. 2014). The iCoh method estimates the partial coherence under a multivariate autoregressive model, and directional paths of interest are kept while irrelevant connection is set to zero. Some methods proposed before the GPDC, such as directed transfer function (DTF) (Kaminski and Blinowska 1991) and directed coherence (DC) (Saito and Harashima 1981), are normalized with respect to a number of series influencing the target variables.

These methods are all based on the assumption of linearity and stationarity. Although they are shown to work with autoregressive time series (Baccalá et al. 2007; Baccalá and Sameshima 2001; Takahashi et al. 2008) and simple nonlinear time series (Wang et al. 2014; Massaroppe and Baccalá 2015; Papana et al. 2013), it is not clear that when applied to neurophysiological signals (Hoerzer et al. 2010; Sato et al. 2009; Shim et al. 2013; Rodrigues and Baccalá 2016), they provided meaningful results. The evaluation of these models using LFP signals generated from simulated spiking networks fills the gap between the hypothesis of the methods and the nonlinearity and nonstationarity of biological LFP signals. The main advantage of simulations is that we have access to all parameters in the model, which allows us to adequately evaluate the GPDC connectivity predictions against the actual network connections.

We considered that the LFP is the sum of input currents in excitatory neurons (Mazzoni et al. 2008, 2015). This assumption is based on the predominant spatial contribution from apical dendrites of pyramidal neurons to the extracellular open-field approximation (Einevoll et al. 2013). However, in physiological experiments, recorded LFP signals are the result of several biophysical mechanisms, such as ionic processes, calcium spikes and neuron-glia interactions (Buzsáki et al. 2012). Surrounding regions may also contribute to LFP signals and, in some cases, even dominate the LFP generated by a local population, such as when surrounding populations receive correlated synaptic inputs (Lindén et al. 2011).

Connectivity between single neurons was already analyzed using Granger causality in simulated neuronal networks to evaluate structural connections (Wu et al. 2011b), synaptic weights (Shao et al. 2015), synaptic plasticity (Cadotte et al. 2008), and point process models  (Kim et al. 2011). Measures based on information theory were applied to infer connections between pairs of neurons in simulated cortical networks (Garofalo et al. 2009; Ito et al. 2011b). Sameshima and Baccalá (1999) applied PDC on neural spiking data from models and rats, after performing a convolution of the spike impulse trains with a given kernel. Our work complements these simulation studies on single cell connectivity by providing an analysis over population LFP signals.

The use of spiking neural networks permits a direct linkage between neuron states and oscillations in the simulated LFP signal. Oscillations seem to influence the flow of information between structurally connected networks (Akam and Kullmann 2014, 2010; Bastos et al. 2015), and we consider that it is fundamental to consider these dynamics when evaluating connectivity methods. Some simulation studies evaluated methods of connectivity applied to simulated electrophysiological data considering mass models of neurons (Lopes et al. 1974; David and Friston 2003; David et al. 2004; Rodrigues and Baccalá 2016). In neural mass models, the dynamic variables represent the average activity of a population of neurons, but in these models, it is not possible to generate states of the cortical networks like synchronization and irregular or regular spike activity of the neurons.

Performance of GPDC was evaluated using ROC curves. We verified that the effectiveness of the method depends on the number of synapses and synaptic weights. Weakly connected networks are difficult to detect using GPDC. This limitation of the method is relevant, given that in biological scenarios some cognitive functions are supported by weak long-range connections (Santarnecchi et al. 2014), which may not be detected by GPDC. Moreover, in neurophysiological experiments, researchers normally define an arbitrary threshold assuming that lower GPDC values may refer to spurious connections (Rubinov and Sporns 2010). But in some scenarios, such as in Fig. 4, small GPDC values may represent coupled networks which are well above the significance threshold. Our results indicate that GPDC is reliable to detect stronger interactions between coupled networks, but may fail to detect weaker connections.

Neurophysiological data are a mixture of signal of interest with an extrinsic noise unrelated to the investigated process (Lowet et al. 2016; Bastos and Schoffelen 2016). We verified that decreasing the SNR degraded the detection rate of GPDC. This result shows a limitation in the connectivity estimation that is common in several connectivity methods (Wu et al. 2011a; Wang et al. 2014) – the presence of noise generally decreases the accuracy of connection inference.

We also found a linear relationship between the coupling strength and GPDC values, with coefficients of determination between 0.58 and 0.77 for the network interactions from model 3. This shows that most of the variation of GPDC values can be explained from the coupling strength, which shows that the GPDC value can be a rough approximation to coupling strength. But as noted before, the GPDC values are influenced by other factors, such as the number of outbound connections. For instance, GPDC values for low coupling strengths are similar to that of uncoupled networks, which explains why GPDC performs poorly in this case.

We conclude that GPDC seems to be reliable in the presence of stronger connections and its magnitude can be used as a rough estimation of connection strength, at least for our simulated models. But with weaker connections or noisier signals, the ROC curves move toward chance level and GPDC connectivity measures should be interpreted with caution.