A growing number of biological experiments are producing datasets consisting of large numbers of interacting variables, from genomics to neural networks to eco-systems, giving rise to the nascent field of systems biology. The eminent challenge of this discipline is how to simplify the analysis of these high dimensional dynamical systems while retaining their relevant features. In particular, despite the rich, complex dynamics of the brain, and the highly interconnected and non-linear nature of its information processing capabilities, the bulk of the literature on brain imaging involves slight variations on the main theme: the identification of the degree of correlation between the activation of a local brain area and external markers. A variety of attempts have been made trying to go beyond this paradigm, including ICA, Volterra kernels and supervised classification techniques [1–4]; however, the applicability of these multi-variate methods has been restricted to networks with only a few hundred vertices.

Recently, a different approach was introduced [5, 6], based on the analysis of the structure of pair-wise correlations between voxels in fMRI, embedding it in a graph structure whose nodes are the voxels and edges between them are defined by the covariance exceeding a threshold. With this relatively simple approach, we were able to demonstrate that the resulting networks have universal statistical topological properties, like scale-free connectivity and small-worldness [7], shared with other biological networks. Moreover, the spatial distribution of these networks is determined by the specific task the brain is engaged in, as opposed to capturing a task-independent network with little functional sensitivity. Similar ideas have been used to analyze and discriminate a variety of functional and dysfunctional brain states, and indeed have confirmed our initial results [8–10].

Another approach to capture the dynamics of complex systems is the causality analysis pioneered by Granger. The essence of this approach is to identify possible causal relationship between the variables of a system by analyzing how much the time course of one variable contributes to that of another one, based on an auto-regressive model. The method has been applied with success to neural data in the context of a few electro-physiological recordings, and to small numbers of brain areas represented by aggregates of several voxels [11, 12]. One aspect of Granger's method and its derivatives that makes it superior to a simple covariance analysis is that it takes care of the transitivity and symmetry of the covariance: if a variable *A* determines the dynamical behavior of variable *B*, and *B* in its turn drives *C*, the covariance will find symmetric "links" *A* ↔ *B*, *B* ↔ *C* and *A* ↔ *C*; similarly, if *A* drives *B* and *C* but the latter two do not influence each other, the covariance will still find B and C to be linked. Multivariate methods based on Granger's causality are able to break the symmetry and the transitivity, and therefore provide a more accurate description of the internal dynamics of the system. These methods have nevertheless several limitations arising from the need of a very large number of samples, as the auto-regressive model requires a matrix whose dimensions are (*M* × *T*)^{2}, where *M* is the number of brain areas, and *T* the number of time units looking into the past [13, 14].

In the present work, we try to circumvent both the limitations of covariance-based analysis and auto-regression-based causality analysis. We extend our previous findings by attempting to capture a larger signature of the dynamics of the brain by including directionality in the edges of the network, based on the concept of delayed covariance.

In our previous work, we defined a functional network by considering all functional voxels {*v*_{
i
}} as possible nodes; their covariance determines whether a binary functional link (or edge) exists between them: *c*_{
ij
}= ⟨(*v*_{
i
}(*t*) - $\overline{v}$_{
i
})(*v*_{
j
}(*t*) - $\overline{v}$_{
j
})${\sigma}_{i}^{-1}{\sigma}_{j}^{-1}$, where $\overline{v}$_{
i
}= ⟨(*v*_{
i
}(*t*)⟩_{
t
}and ${\sigma}_{i}^{2}$ = ⟨(*v*_{
i
}(*t*) - $\overline{v}$_{
i
})^{2}⟩, such that if the correlation between *i* and *j* exceeds a threshold, a functional link is considered, and none otherwise: **if** *c*_{
ij
}> *C*_{
T
}**then** *d*_{
ij
}= 1, **else** *d*_{
ij
}= 0. We extend here this approach by considering the delayed or lagged covariance: *c*_{
ij
}(*τ*) = ⟨(*v*_{
i
}(*t* - *τ*)- $\overline{v}$_{
i
})(*v*_{
j
}(*t*) - $\overline{v}$_{
j
})${\sigma}_{i}^{-1}{\sigma}_{j}^{-1}$. We reason as follows: if there is a significant peak of the covariance between *i* and *j* at zero lag, then there is a potential binary symmetric link between them, as before. However, if the significant peak is not at zero lag, then we will consider that the preceding voxel, and only it, has a directed link pointing to the succeeding one. That is, **if** *c*_{
ij
}(*τ* = 0) > *C*_{
T
}**then** *d*_{
ij
}= *d*_{
ji
}= 1; **else if** *c*_{
ij
}(*τ* > 0) > *C*_{
T
}**then** *d*_{
ij
}= 1 **and** *d*_{
ji
}= 0; **else** *d*_{
ij
}= 0.

In other words, two voxels whose activity is highly correlated and simultaneous are considered to be symmetrically linked; a voxel that is highly correlated with the future of another one will be considered as a "source", and the latter as a "sink". This approach clearly can break the symmetry of the covariance, but as described cannot deal with the problem of the transitivity described above. Taking into consideration the relatively poor temporal resolution of fMRI, we reasoned that for the time being we could only in earnest tackle one of the confounding sources of undirected links, namely the explanation of a zero-lag covariance (i.e. undirected link) between two voxels by the presence of a common source, of which they are both targets. That is, after identifying all sources and sinks, every potential undirected link that can be explained by a common source is removed. We also considered possible reductions of triangulations of directed links (as in *A* → *B*, *B* → *C*, *A* → *C*), but we did not find this approach to have a significant impact, presumably due to the low temporal resolution of the signal (data not shown, see the Methods section).