Synergy and redundancy in the Granger causal analysis of dynamical networks

Living organisms can be modeled as an ensemble of complex physiological systems, each with its own regulatory mechanism and all continuously interacting between them [1]. Therefore, inferring the interactions within and between these modules is a crucial issue. Over the past few years, the interaction structure of many complex systems has been mapped in terms of networks, which have been successfully studied using tools from statistical physics [2]. Dynamical networks have modeled physiological behavior in many applications; examples range from networks of neurons [3], genetic networks [4], protein interaction nets [5], and metabolic networks [6–8].

The inference of dynamical networks from time series data is related to the estimation of the information flow between variables [9]; see also [10, 11]. Granger causality (GC) [12, 13] has emerged as a major tool to address this issue. This approach is based on prediction: if the prediction error of the first time series is reduced by including measurements from the second one in the linear regression model, then the second time series is said to have a Granger causal influence on the first one. It has been shown that GC is equivalent to transfer entropy [14] in the Gaussian approximation [15] and for other distributions [16]. See [17] for a discussion about the applicability of this notion in neuroscience and [18] for a discussion on the reliability of GC for continuous dynamical processes. It is worth stressing that several forms of coupling may mediate information flow in the brain, see [19, 20]. The combination of GC and complex networks theory is also a promising line of research [21].

The pairwise Granger analysis (PWGC) consists of assessing GC between each pair of variables independently from the rest of the system. It is well known that PWGC cannot disambiguate direct and indirect interactions among variables. The most straightforward extension, the conditioning approach [22], removes indirect influences by evaluating to what extent the predictive power of the driver on the target decreases when the conditioning variable is removed. It must be noted however that, even though its limitations are well known, the PWGC approach is still used in situations where the number of samples is limited and a fully conditioned (CGC) approach is unfeasible. As a convenient alternative to this suboptimal solution, a partially conditioned approach, consisting of conditioning on a small number of variables, chosen as the most informative one for the driver node, has been proposed [23]; this approach leads to results very close to those obtained with a fully conditioned analysis that are even more accurate in the presence of a small number of samples [24]. We remark that the use of partially conditioned GC (PCGC) may also be useful in non-stationary conditions whereas the GC pattern must be estimated on short time windows.

However, sometimes a CGC approach can encounter conceptual limitations in addition to the practical and computational ones: in the presence of redundant variables, the application of the standard analysis leads to an underestimation of influences [25]. Redundancy and synergy are intuitive, yet elusive, concepts that have been investigated in different fields, from pure information theory [26–28], to machine learning [29] and neural systems [30, 31], with definitions that range from purely operative to conceptual. When analyzing interactions in a multivariate time series, redundancy may arise if some channels are all influenced by another signal that is not included in the regression; another source of redundancy may be the emergence of synchronization in a subgroup of variables, without the need for an external influence. Redundancy manifests itself through a high degree of correlation in a group of variables, both for instantaneous and lagged influences. Several approaches have been proposed to reduce dimensionality in multivariate sets that eliminate redundancy, which rely on generalized variance [32], principal components analysis [33], or GC itself [34].

A complementary concept to redundancy is synergy. The synergetic effects that we address here, related to the analysis of dynamical influences in a multivariate time series, are similar to those encountered in sociological and psychological modeling, where suppressors are variables that increase the predictive validity of another variable after its inclusion into a linear regression equation [35]; see [36] for examples of easily explainable suppressor variables in multiple regression research. Redundancy and synergy have been further connected to information transfer in [37], where an expansion of the information flow has been proposed to prove redundant and synergetic multiplets of variables. Other information-based approaches have also addressed the issue of collective influence [28, 38]. The purpose of this paper is to provide evidence that, in addition to the problem related to indirect influence, PWGC shows another relevant pitfall: it fails to detect synergetic effects in the information flow. In other words, it does not account for the presence of subsets of variables that provide some information about the future of a given target only when all the variables are used in the regression model. We remark that, because it processes the system as a whole, CGC evidences synergetic effects; when the number of samples is low, PCGC can also detect synergetic effects after an adequate selection of the conditioning variables.

The paper is organized as follows. In the next section, we briefly recall the concepts of GC and PCGC. In section 3, we describe some toy systems that illustrate how redundancy can affect the results of CGC, whereas indirect interactions and synergy are the main problems inherent in PWGC. In section 4, we provide evidence of synergetic effects in epilepsy: we analyze EEG recordings of an epileptic patient corresponding to ten seconds before the seizure onset; we show that the two contacts that constitute the putative seizure onset act as synergetic variables driving the rest of the system. The pattern of information transfer proves the actual seizure onset only when synergy is correctly considered.

In section 5, we propose an approach that combines PWGC and GCC to prove the pairwise influences are only because of redundancy and not recognized by CGC. The conditioned GC pattern is used to partition the pairwise links in two sets: those that are indirect influences between the measured variables, according to CGC, and those that are not explained as indirect relationships. The unexplained pairwise links, presumably due to redundancy, are thus retained to complement the information transfer pattern discovered by CGC. In cases where the number of samples is so low that a fully multivariate approach is unfeasible, PCGC may be applied instead of CGC. We also address here the issue of variable selection for PCGC, and consider a novel strategy for the selection of variables: for each target variable, one selects the variables sending the highest amount of information to that node, as revealed by a pairwise analysis. By construction, this new selection strategy works more efficiently when the interaction graph has a tree structure: indeed, in this case, conditioning on the parents of the target node ensures that indirect influences will be removed. In the epilepsy example, the selection based on the mutual information with the candidate driver provides results closer to those obtained by CGC. Finally we apply the proposed approach to time series of gene expressions, extracted from a data set from the HeLa culture. Section 6 summarizes our conclusions.

GC is a powerful and widespread data-driven approach used to determine whether and how two time series exert direct dynamical influences on each other [39]. A convenient nonlinear generalization of GC has been implemented in [40], exploiting the kernel trick, which makes computation of dot products in high-dimensional feature spaces possible using simple functions (kernels) defined on pairs of input patterns. This trick allows the formulation of nonlinear variants of any algorithm that can be cast in terms of dot products, for example support vector machines [41]. Hence, in [40], the idea is still to perform linear GC but in a space defined by the nonlinear features of the data. This projection is conveniently and implicitly performed through kernel functions [42] and a statistical procedure is used to avoid over-fitting.

Quantitatively, let us consider n time series ${{\{{{x}_{\alpha }}(t)\}}_{\alpha =1,\ldots ,n}}$; the lagged state vectors are denoted

$$\begin{eqnarray*}&&{{X}_{\alpha }}(t)=\left( {{x}_{\alpha }}(t-m),\ldots ,{{x}_{\alpha }}(t-1) \right),\end{eqnarray*}$$

where m is the order of the model (window length). Let $\epsilon \left( {{x}_{\alpha }}|{\bf X} \right)$ denote the mean squared error prediction of ${{x}_{\alpha }}$ on the basis of all the vectors ${\bf X}=\{{{X}_{\beta }}\}_{\beta =1}^{n}$ (corresponding to the kernel approach described in [43]). The multivariate GC index ${{\delta }_{mv}}(\beta \to \alpha )$ is defined as follows: consider the prediction of ${{x}_{\alpha }}$ on the basis of all variables except ${{X}_{\beta }}$ and the prediction of ${{x}_{\alpha }}$ using all variables, then the GC is the (normalized) variation of the error in the two conditions, i.e.,

$$\begin{eqnarray}&&{{\delta }_{mv}}(\beta \to \alpha )={\rm log} \frac{\epsilon \left( {{x}_{\alpha }}|{\bf X}\backslash {{X}_{\beta }} \right)}{\epsilon \left( {{x}_{\alpha }}|{\bf X} \right)};\end{eqnarray} \tag{ 1 }$$

In [44], it has been shown that not all the kernels are suitable to estimate GC. Two important classes of kernels that can be used to construct nonlinear GC measures are the inhomogeneous polynomial kernel (whose features are all the monomials in the input variables up to the pth degree; p = 1 corresponds to linear GC) and the Gaussian kernel.

The PCGC is given by:

$$\begin{eqnarray}&&{{\delta }_{bv}}(\beta \to \alpha )={\rm log} \frac{\epsilon \left( {{x}_{\alpha }}|{{X}_{\alpha }} \right)}{\epsilon \left( {{x}_{\alpha }}|{{X}_{\alpha }},{{X}_{\beta }} \right)}.\end{eqnarray} \tag{ 2 }$$

The PWGC is defined as follows. Let ${\bf Y}$ denote the variables in ${\bf X}$, excluding ${{X}_{\alpha }}$ and ${{X}_{\beta }}$. Then equation (1) can be written as:

$$\begin{eqnarray}&&\delta _{c}^{{\bf Y}}(\beta \to \alpha )={\rm log} \frac{\epsilon \left( {{x}_{\alpha }}|{{X}_{\alpha }},{\bf Y} \right)}{\epsilon \left( {{x}_{\alpha }}|{{X}_{\alpha }},{{X}_{\beta }},{\bf Y} \right)}.\end{eqnarray} \tag{ 3 }$$

When ${\bf Y}$ is only a subset of the total number of variables in ${\bf X}$ not containing ${{X}_{\alpha }}$ and ${{X}_{\beta }}$, then $\delta _{c}^{{\bf Y}}$ is called the PCGC. In [23], the set ${\bf Y}$ is chosen as the most informative for ${{X}_{\beta }}$. Here, we will also consider an alternative strategy: fixed a small number k, we select ${\bf Y}=\{{{{\bf X}}_{\gamma }}\}_{\gamma =1}^{k}$ as the k variables with the maximal pairwise GC ${{\delta }_{bv}}(\gamma \to \alpha )$ w.r.t. that target node, excluding ${{X}_{\beta }}$.

In this section, we provide some typical examples to note possible problems that pairwise and fully conditioned analysis may encounter.

3.1. Indirect GC among measured variables

We consider the following lattice of ten unidirectionally coupled noisy logistic maps, with

$$\begin{eqnarray}&&{{x}_{1}}(t)=f\left( {{x}_{1}}(t-1) \right)+0.01{{\eta }_{1}}(t),\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{{x}_{i}}(t)=(1-\rho )f\left( {{x}_{i}}(t-1) \right)+\rho f\left( {{x}_{i-1}}(t-1) \right)+0.01{{\eta }_{i}}(t),\end{eqnarray} \tag{ 5 }$$

with $i=2,\ldots ,10$. Variables η are unit variance Gaussian noise terms. The transfer function is given by $f(x)=1-1.8{{x}^{2}}$. In this system, the first map is evolving autonomously, whereas the other maps are influenced by the previous maps with coupling ρ, thus forming a cascade of interactions. In figure 1(a), we plot, as a function of ρ, the number of GC interactions found by PWGC and CGC using the method described in [25] with the inhomogeneous polynomial kernel of degree two. The CGC output is the correct one (nine links), whereas the PWGC output also accounts for indirect influences and therefore fails to provide the underlying network of interactions. In this example we also tested PCGC (see figure 1(b)). We considered only one conditioning variable, chosen according to the two strategies described previously. First, we consider the most informative w.r.t. the candidate driver, as described in [23]; we call this information-based (IB) strategy. Second, we choose the variable characterized by the maximal pairwise influence to the target node, a pairwise–based (PB) rule. The PB strategy yields the correct result in this example, whereas the IB strategy fails when only one conditioning variable is used and requires more than one conditioning variable to provide the correct output. This occurrence is due to the tree topology of the interactions in this example, which favors PB selecting by construction the parents of each node.

Zoom In Zoom Out Reset image size

3.2. Redundancy due to a hidden source

3.3. Synergetic contributions

Let us consider three unit variance i.i.d. Gaussian noise terms x₁, x₂, and x₃. Let

$$\begin{eqnarray*}&&{{x}_{4}}(t)=0.1({{x}_{1}}(t-1)+{{x}_{2}}(t-1))+\rho {{x}_{2}}(t-1){{x}_{3}}(t-1)+0.1\eta (t).\end{eqnarray*}$$

Considering the influence $3\to 4$, the CGC reveals that 3 is influences 4, whereas PWGC fails to detect this causal relationship—see figure 1(d), in which we use the method described in [25] with the inhomogeneous polynomial kernel of degree two, where x₂ is a suppressor variable for x₃ w.r.t. the influence on x₄. This example shows that PWGC fails to detect synergetic contributions. We note that use of nonlinear GC is mandatory in this case to verify the synergy between x₂ and x₃.

3.4. Redundancy due to synchronization

As another example, we consider a toy system made of five variables $\{{{x}_{i}}\}$. The first four variables constitute a multiplet made of a fully coupled lattice of noisy logistic maps with coupling ρ, evolving independently of the fifth variable. The fifth variable is influenced by the mean field from the coupled map lattice. The equations are, for $i=1,\ldots ,4$:

$$\begin{eqnarray}&&{{x}_{i}}(t)=(1-\rho )f({{x}_{i}}(t-1))+\rho \sum _{j=1,j\ne i}^{4}f({{x}_{j}}(t-1))+0.01{{\eta }_{i}}(t),\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{{x}_{5}}(t)=\sum _{i=1}^{4}\frac{{{x}_{i}}(t-1)}{8}+{{\eta }_{5}}(t),\end{eqnarray} \tag{ 7 }$$

where η denotes unit variance Gaussian noise terms. Increasing the coupling ρ among the variables in the multiplet $\{{{x}_{1}},{{x}_{2}},{{x}_{3}},{{x}_{4}}\}$, the synchronization of these variables (measured by Pearson correlations) increases and they become nearly synchronized for ρ greater than 0.1 (perfect synchronization cannot be achieved due to the noise terms); redundancy, in this example, arises due to the complex inherent dynamics of the units. In figure 2, we depict both the causality from one variable in the multiplet (x₁; the same results hold for x₂, x₃, and x₄ ) to x₅, and the causality between pairs of variables in the multiplet: both linear and nonlinear PWGC and CGC are shown for the two quantities.

Zoom In Zoom Out Reset image size

Concerning the causality ${{x}_{1}}\to {{x}_{5}}$, we note that, for low coupling, both PWGC and CGC, with a linear or nonlinear kernel, correctly detect the causal interaction. Around the transition to synchronization, in a window centered at $\rho =0.05$, all algorithms fail to detect the causality ${{x}_{1}}\to {{x}_{5}}$. In the almost synchronized regime, $\rho \gt 0.1$, CGC continues to fail due to redundancy, whereas the PWGC correctly provides the causal influence, using both the linear and the nonlinear algorithm.

As far as the causal interactions within the multiplet are concerned, we note that using the linear approach, we obtain small values of causality just at the transition, whereas we obtain zero values far from the transition. Using the nonlinear algorithm, which is the correct one in this example because the system is nonlinear, we obtain a nonzero causality among the variables in the multiplet, using both PWGC or CGC: the resulting curves are non-monotonous, as one may expect, due to the inherent nonlinear dynamics. For $\rho \gt 1$, nonzero GC is observed because of the noise that prevents the system to go in the complete synchronized state.

This example again shows that, in the presence of redundancy, one should consider both CGC and PWGC results. Moreover, it also shows how nonlinearity may render extremely difficult the inference of interactions: in this system, there is a range of values, corresponding to the onset of synchronization, in which all methods fail to provide the correct causal interaction.

As a real example, we consider intracranial EEG recordings from a patient with drug-resistant epilepsy with an implanted array of 8 × 8 cortical electrodes (CE) and two depth electrodes (DE) with six contacts each. The data are available in [45] and further details on the dataset are given in [46]. Data were sampled at 400 Hz. Here we consider a portion of data recorded in the preictal period, 10 seconds preceding the seizure onset. To handle this data, we use linear GC with m equal to five. In figure 3, we depict the PWGC between DEs (panel a), from DEs to CEs (panel b), between CEs (panel c), and from CEs to DEs (panel d). We note a complex pattern of bivariate interactions among CEs, whereas the first DE seems to be the subcortical drive to the cortex. We note that there is no PWGC from the last two contacts of the second DE (channels 11 and 12) to CEs and neither to the contacts of the first DE. In figure 4, we depict the CGC among DEs (panel a), from DEs to CEs (panel b), among CEs (panel c), and from CEs to DEs (panel d). The scenario in the conditioned case is clear: contacts 11 and 12, from the second DE, are the drivers both for the cortex and for the subcortical region associated to the first DE. These two contacts can be then associated to the seizure onset zone (SOZ). The PWGC strength among CEs is due to redundancy, as these are all driven by the same subcortical source. Because contact 12 is also driving contact 11 (see figure 4(a)), we conclude that contact 12 is the closest to the SOZ and that contact 11 is a suppressor variable for it, because it is necessary to include it in the regression model to verify the influence of contact 12 on the rest of the system. Conversely, contact 12 acts as a suppressor for contact 11. We stress that the influence from contacts 11 and 12 to the rest of the system emerges only in the CGC and is neglected by PWGC: these variables are synergetically influencing the dynamics of the system. To our knowledge, this is the first time that synergetic effects are found in relation with epilepsy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also apply PCGC on this data using one conditioning variable. The results are depicted in figure 5: using the IB strategy, we obtain a pattern very close to the one from CGC, whereas this is not the case of PB. These results seem to suggest that IB works better in the presence of redundancy; however, we do not offer arguments to claim that this a general rule. It is worth mentioning that in the presence of synergy, the selection of variables for partial conditioning is equivalent to the search of suppressor variables.

Zoom In Zoom Out Reset image size

In the preceding sections, we have shown that CGC encounters issues resulting in poor performance in presence of redundancy, and that information about redundancy may be obtained from the PWGC pattern. Here we develop a strategy to combine the two approaches: some links inferred from PWGC are retained and added to those obtained from CGC. The PWGC links that are discarded are those that can be derived as indirect links from the CGC pattern. In the following, we describe the proposed approach in detail.

Let Δ denote the matrix of influences from CGC (or PCGC). Let ${{\Delta }^{*}}$ denote the matrix from PWGC. Nonzero elements of Δ and ${{\Delta }^{*}}$ correspond to the estimated influences. Let these matrices be evaluated using a model of order m.

The matrix

$$\begin{eqnarray*}&&{{M}_{\alpha \beta }}={{\Delta }^{\alpha }}{{\left( {{\Delta }^{\top }} \right)}^{\beta }}\end{eqnarray*}$$

contains paths of length $\alpha +\beta$ with delays in the range $[-\beta m+\alpha ,\ldots ,-\beta +\alpha m]$. Indeed:

$$\begin{eqnarray}&&{{M}_{\alpha \beta }}(i,j)={{\sum }_{{{i}_{1}}}}{{\sum }_{{{i}_{2}}}}\cdots {{\sum }_{{{i}_{\alpha +\beta -1}}}}\Delta (i,{{i}_{1}})\Delta ({{i}_{1}},{{i}_{2}})\cdots \Delta (j,{{i}_{\alpha +\beta -1}});\end{eqnarray} \tag{ 8 }$$

because all elements of Δ are non-negative, it follows that ${{M}_{\alpha \beta }}(i,j)$ is non-vanishing if, and only if it is possible, in matrix Δ, to go from node i to node j moving $\beta$ steps backward and α steps forward, where a step is allowed if the corresponding element of Δ is not zero. Therefore, the nonzero elements of the matrix ${{M}_{\alpha \beta }}$ describe a situation where two nodes receive a common input from a third node that is α steps backward in time from one node and $\beta$ steps backward in time from the other node. In other words, if the element ${{M}_{\alpha \beta }}(i,j)$ does not vanish, then there exists an indirect interaction between nodes i and j due to a common input. The circuit corresponding to ${{M}_{21}}$ is represented in figure 6(a): if the element ${{M}_{2,1}}(i,j)$ is non-vanishing, then i and j are connected as shown in figure 6(a).

Zoom In Zoom Out Reset image size

Because the order of the model is m, a simple comparison between the delays from the common source to i and j demonstrates that the indirect influence corresponding to the nonzero element ${{M}_{\alpha ,\beta }}(i,j)$ might be detected by PWGC only if

$$\begin{eqnarray*}&&[-\beta m+\alpha ,\ldots ,-\beta +\alpha m]\cap [1,\ldots ,m]\ne \varnothing ;\end{eqnarray*}$$

this is equivalent to

$$\begin{eqnarray*}&&\frac{\beta +1}{m}\leqslant \alpha \leqslant (\beta +1)m.\end{eqnarray*}$$

Now, the matrix ${{F}_{\alpha }}={{\Delta }^{\alpha }},$ with $\alpha \geqslant 1$, contains paths of length α with delays in the range $[\alpha ,\ldots ,\alpha m]$, indeed:

$$\begin{eqnarray}&&{{F}_{\alpha }}(i,j)={{\sum }_{{{i}_{1}}}}{{\sum }_{{{i}_{2}}}}\cdots {{\sum }_{{{i}_{\alpha -1}}}}\Delta (i,{{i}_{1}})\Delta ({{i}_{1}},{{i}_{2}})\cdots \Delta ({{i}_{\alpha -1}},j);\end{eqnarray} \tag{ 9 }$$

Any nonzero element of the matrix ${{F}_{\alpha }}(i,j)$ describes an indirect causal interaction between nodes i and j, where i sends information to j through a cascade of α links: i→i₁, ${{i}_{1}}\to {{i}_{2}}$, ..., ${{\imath }_{\alpha -1}}\to j$. The circuit corresponding to F₂ is depicted in figure 6(b). The indirect causal interaction $i\to j$, corresponding to the nonzero element ${{F}_{\alpha }}(i,j)$, might be detected by PWGC if $\alpha \leqslant m$.

Let us now consider the matrix

$$\begin{eqnarray*}&&B={{\sum }_{\alpha ,\beta }}{{M}_{\alpha \beta }}+\sum _{{{\alpha }^{\prime }}=1}^{m}{{F}_{{{\alpha }^{\prime }}}},\end{eqnarray*}$$

where the first sum is over pairs $\{\alpha ,\beta \}$ satisfying $\frac{\beta +1}{m}\leqslant \alpha \leqslant (\beta +1)m.$ If B_ij is non-vanishing, then according to CGC, there is an indirect causal interaction between i and j: therefore, PWGC might misleadingly reveal such interact, considering it a direct one. In the approach previously described, we discard (as indirect) the links found by PWGC for which $B(i,j)\ne 0$. Therefore, in the pairwise matrix ${{\Delta }^{*}}$, we set to zero all the elements such that ${{B}_{ij}}\gt 0$ (pruning). The resulting matrix ${{\Delta }^{*}}$ contains links which cannot be interpreted as indirect links of the multivariate pattern, and will be retained and ascribed to redundancy effects,

For m = 1, we have that the only terms in the first sum are those with $\alpha =\beta +1$, so the first non-trivial terms are

$$\begin{eqnarray*}&&{{B}_{1}}=\Delta +{{\Delta }^{2}}{{\Delta }^{\top }}.\end{eqnarray*}$$

For m = 2, the simplest terms are:

$$\begin{eqnarray*}&&{{B}_{2}}=\Delta +{{\Delta }^{2}}+\Delta {{\Delta }^{\top }}.\end{eqnarray*}$$

Because, due to the finite number of samples, a mediated interaction is more unlikely to be detected (by the pairwise analysis) if it corresponds to a long path, we limit the sum in the matrix B to the simplest terms.

As a toy example to illustrate an application of the proposed approach, we consider a system made of five variables $\{{{x}_{i}}\}$. The first four variables constitute a multiplet made of unidirectionally coupled logistic maps, equations (4)–(5) with i ranging in $\{1,2,3,4\}$, coupling ρ, and interactions $1\to 2$, $2\to 3$, and $3\to 4$. The fifth variable is influenced by the mean field from the coupled map lattice, see equation (7). The four variables in the multiplet become almost synchronized for $\rho \gt 0.4$. In figure 7, we depict both the average influence from the variables in the multiplet to x₅ and the average influence between pairs of variables in the multiplet: both linear and nonlinear PWGC and CGC are shown for the two quantities. Note that only the nonlinear algorithm correctly verifies the causal interactions within the multiplet of four variables, whereas the linear algorithm detects a very low causal interdependency among them. The driving influence from the multiplet to x₅ detected by CGC vanishes at high coupling redundancy, both in the linear and nonlinear approach, due to the redundancy induced by synchronization.

Zoom In Zoom Out Reset image size

To explain how the proposed approach works, we describe two situations corresponding to low and high coupling, respecively. At low coupling, the CGC approach estimates the correct causal pattern in the system, and the nonzero elements of Δ are $1\to 2$, $2\to 3$, $3\to 4$, $1\to 5$, $2\to 5$, $3\to 5$, and $4\to 5$. The nonzero elements of the matrix ${{\Delta }^{*}}$, from PWGC analysis, are the same as Δ plus $1\to 3$, $1\to 4$, $2\to 4$, corresponding to indirect causalities; however, these three interactions lead to nonzero elements of Δ² (and, therefore, of B), hence they must be discarded. It follows that ${{\Delta }^{*}}$ does not provide further information than Δ at low coupling.

On the contrary, at high coupling, due to synchronization, the CGC approach does not reveal the causal interactions $1\to 5$, $2\to 5$, $3\to 5$, and $4\to 5$, although still they are recognized by PWGC. Here, ${{\Delta }^{*}}$ is still nonzero in correspondence to $1\to 5$, $2\to 5$, $3\to 5$, and $4\to 5$, while the corresponding elements of B are vanishing. According to our previous discussion, the interactions $1\to 5$, $2\to 5$, $3\to 5$, and $4\to 5$, detected by PWGC, should not be discarded: combining the results from CGC and PWGC, we obtain the correct causal pattern even in presence of strong synchronization.

HeLa [47] is a famous cell culture, isolated from a human uterine cervical carcinoma in 1951. HeLa cells have acquired cellular immortality, in that the normal mechanisms of programmed cell death after a certain number of divisions have somehow been switched off. We consider the HeLa cell gene expression data of [48]. Data corresponds to 94 genes and 48 time points, with an hour interval separating two successive readings (the HeLa cell cycle lasts 16 hours). The 94 genes were selected from the full data set described in [49], on the basis of the association with cell cycle regulation and tumor development. We apply linear PWGC and linear CGC (using only another conditioning variable and using both IB and PB selection strategies, described in section 3). We note that the CGC approach is unfeasible in this case due to the limited number of samples. In this case, due to the limited number of samples, we do not use statistical testing assess the significance of the retrieved links. Rather, we introduce a threshold for the influence and analyze the pattern as the threshold is varied. In figure 8, results are reported as a function of the number of links found by PWGC n_pairwise (which increases as the threshold is decreased); we plot (1) the number of links found by PGC n_partial, (2) the number of links found by PGC and not by PWGC n_novel, which are thus a signature of synergy, (3) the percentage of pairwise links that can be explained as direct or indirect causalities of the PGC pattern (thus being consistent with the partial causality pattern), found using the matrix B₁ to detect the indirect links, which correspond to circuits like the one described in figure 6(a) and (4) the number of causality links found by PWGC and not consistent with PWGC n_unexplained, corresponding to redundancy. The two curves refer to the two selection strategies for partial conditioning.

Zoom In Zoom Out Reset image size

The low number of samples here allowed us to use only one conditioning variable, and therefore, to analyze the circuits of only three variables; a closely related analysis [50] has been proposed to study how a gene modulates the interaction between two other genes. Conversely, the true underlying gene regulatory network being unknown, we cannot assess the performances of the algorithms in terms of correctly detected links.

We note that both n_novel and n_unexplained assume relatively large values; hence, both redundancy and synergy characterize this data set. The PB selection strategy yields slightly higher values of n_novel and n_unexplained, emerging as a better discriminator of synergy and redundancy than IB. A comparison with the fully conditioned approach is not possible in this case. On the other hand, as far as the search for synergetic effects is concerned, we find that the synergetic interactions found by PCGC with the two strategies are not coinciding; indeed, only 10% of all synergetic interactions are found by both strategies. This suggests that, when searching for suppressors, several sets of conditioning variables should be used in CGC to explore more possible alternative pathways, especially when there is not a priori information on the network structure.

In this paper, we considered the inference, from time series data, of the information flow between subsystems of a complex network, which is an important problem in medicine and biology. In particular, we analyzed the effects that synergy and redundancy induce on the Granger causal analysis of time series.

Concerning synergy, we have shown that the search for synergetic contributions in the information flow is equivalent to the search for suppressors, i.e., variables that improve the predictive validity of another variable. Pairwise analysis fails to verify this kind of variable. However, fully multivariate GC solves this problem: conditioning on suppressor variables leads to nonzero GC. In cases when the number of samples is low, we have shown that PCGC is a valuable option, provided that the selection strategy to choose the conditioning variables, succeeds in picking the suppressors. In this paper, we considered two different strategies: choosing the most informative variables for the candidate driver node, or choosing the nodes with the highest pairwise influence to the target. From the several examples analyzed here, we have shown that the first strategy is viable in the presence of redundancy, whereas when the interaction pattern has a tree-like structure, the latter is preferable. However, the issue of selecting variables for PCGC deserves further attention, as it corresponds to the search for suppressor variables and correspondingly of synergetic effects. We have also provided evidence, for the first time, that synergetic effects are present in an epileptic brain in the preictal condition (just prior to seizure).

We have then shown that CGC approaches do not work well in the presence of redundancy. To handle redundancy, we propose to split the pairwise links into two subsets: those that correspond to indirect connections of the multivariate GC, and those that do not. The links that are not explained, as indirect connections are ascribed to redundancy effects and they are merged to those from CGC to provide the full causality pattern in the system. We have applied this approach to a genetic data set from the HeLa culture, and found that the underlying gene regulatory networks are characterized by both redundancy and synergy, hence these approaches are also promising w.r.t. the reverse engineering of gene regulatory networks.

In conclusion, we observe that the problem of inferring reliable estimates of the causal interactions in real dynamical complex systems, when limited a priori information is available, remains a major theoretical challenge. Recently, the most important results in this direction are related to the use of data-driven approaches, such as GC and transfer entropy. In this work, we have shown that in the presence of redundancy and synergy, combining the results from the pairwise and conditioned approaches may lead to more effective analyses.