Modelling a multiplex brain network by local transfer entropy

This paper deals with the information transfer mechanisms underlying causal relations between brain regions under resting condition. fMRI images of a large set of healthy individuals from the 1000 Functional Connectomes Beijing Zang dataset have been considered and the causal information transfer among brain regions studied using Transfer Entropy concepts. Thus, we explored the influence of a set of states in two given regions at time t (At Bt.) over the state of one of them at a following time step (Bt+1) and could observe a series of time-dependent events corresponding to four kinds of interactions, or causal rules, pointing to (de)activation and turn off mechanisms and sharing some features with positive and negative functional connectivity. The functional architecture emerging from such rules was modelled by a directional multilayer network based upon four interaction matrices and a set of indexes describing the effects of the network structure in several dynamical processes. The statistical significance of the models produced by our approach was checked within the used database of homogeneous subjects and predicts a successful extension, in due course, to detect differences among clinical conditions and cognitive states.

1 Workflow of the analytical approach used in this paper. Causal Interaction among pairs of nodes by the TRANSFER ENTROPY formalism [A n , B n ] à B n+1 FILTERING and NAMING RULES -8 rules emerge by distribution among subjects (Table 1) -4 rules among those are endowed with function-related names : ActS, ActO, TfS,TfO (Table 2) MULTIPLEX NETWORK ANALYSIS -Network dynamics, Centrality and Modularity -Definition of a, b and g hierarchical levels (Fig. 8) -Defininition of 0 th , 1 st and 2 nd order subgraphs (Fig. 3) Tracing single links among nodes (Tables 1 , 2 in Supplementary Materials) Figure S1: Flowchart of the analytical approach followed in this paper.
-Yellow boxes = Preliminary steps (Section 2.1 of the paper): data collection and preprocessing and signal extraction by the defined 90 ROIs, .
-Green boxes = Characterizing the dyadic brain region interactions (Section 2.2 of the paper and Section 3 of Suppl. Mat., below): a) time-series discretization and characterization of the three functional states (1, 0 and -1); b) estimating the causal interaction by the Transfer Entropy (TE) formalism of a whole set of 27 state combinations (rules); c) finding a significant subset of function-related rules, namely ActS, TfS, ActO and TfO. -Blue boxes = Multiplex network analysis (section 2.3 of the paper): a) characterizing the dynamics of the mutual interactions among ROIs from (90x89=8010, self-interaction excluded) adjacency matrices; b) studying the adjacency matrices associated to the four above rules by four network indexes, two of them to estimate dynamical features (reciprocity and subgraphs), and the other two for the centrality and modularity analysis; c) exploring the single links among ROIs in order to define possible functional circuits; -Red Box = Comparative and critical final considerations (section 4.3.2 in this paper).

Statistical significance of dyadic interactions
Given the 90 brain regions [1] in Figure S2, in order to study the significant distribution among subjects of the 27 (3 3 ) rules from the dyadic interactions, for each subject a square matrix of 8010 local Transfer Entropy (TE) [2] values (excluding 90 self-interactions) was considered and the rules associated to random events were filtered out by a series of statistical tests. Thus, since only the time-sequences are randomized, the marginal frequency of events in each brain region is unchanged, and the local TE could be re-calculated for each brain region combination, for all subjects and rules. At increasing signal threshold the significant events increase [3] [4], and the noise over the information transfer among brain regions is reduced, as confirmed by the Mean Square Error (MSE) between the local TE values of each rule for each subject and the corresponding random model. For some rules an increasing trend in both the sum of local TE and MSE points to a significant improvement of the information transfer at increasing signal threshold. Those rules should emerge as different from other rules at the highest signal threshold (= 1 S.D.), which is confirmed by a one way ANOVA (shown in the Results section of the main text).
Since MSE only accounts for symmetrical deviations from the random model, to identify rules with an increasing trend of TE values, a t-student for independent samples (Bonferroni corrected for the number of subjects) was performed. The rationale is that a significant decrease with respect to the random model indicates disadvantaged rules, in contrast with rules showing a significant increases of TE values, which are advantaged and possibly endowed with a functional meaning. In conclusion, we take as significant the rules characterized by the following properties: -increasing trend of TE values in parallel with signal threshold; -increasing trend of MSE in parallel with z-score threshold; -significantly higher values of local TE when compared with other rules; -significant probability in a right-tail statistical comparison test. In Figure S3 a concise scheme of the whole procedure is reported.
[Index] Figure S3: Statistical significance of the rules. To evaluate the significance of rules the following steps were used: -in the preliminary step (top center) the time-series at hand were randomized preserving the marginal frequencies of events (1, -1 and 0); -in step 1 (bottom left ) for each subject the TE average values for each rule (27) and for each threshold (0.25, 0.50, 0.75, 1) were calculated; -in step 2 (bottom right) the average Mean Square Error (MSE) between original data and random model were calculated for all rules and thresholds; -in the intermediate step only rules having an increasing trend of values (calculated in step 1 and 2) as a function of threshold were selected, that is 12 rules out of 27; -in step 3 (bottom right) for each subject the original and randomized data at the highest threshold (1) were analyzed by a t test right-tail, in order to filter out rules having TE values higher than the random model. -In the very last step 8 rules out of 12 showed significant results in the largest fraction of subjects. [Index] 4 Subgraphs in a multiplex network Different names are possible for a 3-nodes subgraph [5] [6] [7] due to its symmetric features, e.g. ActO/ActS/ActO = ActS/ActO/ActO = ActO/ActO/ActS. In any case, we use the standardized name shown in the left panel of Figure S4 and Figure S5.
In a Cycle triad the information transfer starts and ends at the same node. Conversely, in a Flux triad the information transfer starts from a source node and ends to a sink node, passing through a middle node, and the process was called feedforward loop [8]. Thus, in a Flux subgraph we characterize: -a source node with only out-connections; -a sink node with only in-connections; -a middle node with both in-and out-connections. Introducing Flux subgraphs in a multiplex frame needs to distinguish links passing between: -source and middle nodes; -middle and sink nodes; -source and sink nodes.
Thus, a Flux subgraph is identified by the layers having a link corresponding to the previously defined cases and in the same order.
For example, the 2 nd order subgraph: Flux ActO/ActS/TfO is formed by: 1) a link between source and middle nodes falling in the ActO layer; 2) a link between middle and sink nodes belonging to the ActS layer, and 3) a link between source and sink nodes belonging to the TfO layer.  Using three layers the symmetry found in the 1 nd order subgraphs ( Figure S4) is lost and different dynamic processes are included in the isophorms. [Index] 5 Multireciprocity and synchronous or sequential mechanisms.
In order to characterize possible dynamical mechanisms underlying the network architecture, the logical functioning of the rules can be calculated between the interacting rules in their multireciprocal structure. Two kinds of dynamics can be assumed: synchronous and sequential. In the first case nodes can act in the same time step changing their state simultaneously. The only combination for such a case is the reciprocal turn off rules in their own layer: TfS/TfS and TfO/TfO. However, this combination is not-reciprocal in our results, and both Turn Off rules appear unidirectional. In the second case, node A could act on node B at time step n, and node B on node A at the next n+1 step, pointing to a temporal sequence changing the state of both A and B.
Results of the multireciprocity analysis show the following significant reciprocal interaction rules: TfS/TfO, ActS/TfO, ActO/TfS. As shown in the upper panel of Figure S6, the TfS/TfO rules cannot act at the same time: TfS need the node state of A=1 (or -1) and B=-1 (or 1), while the TfO the node state A=1 (of -1) and B=1 (or -1), then since they need different combination of nodes state to work, they cannot operate at the same time. On the other hand, in the turn off rules the node A change the state of the node B in 0, and vice versa, then the node B does not available to carry any reciprocal action. Consequently, TfS and TfO cannot describe a reciprocal dynamics, but operate independently reciprocally. In the latter cases, both ActS/TfO and ActO/TfS share same properties: they cannot operate synchronously, but a sequence dynamics is possible. As described in the bottom panel of Figure S6, if the node A acts the ActS (or ActO) rule on the node B, the state of node B change to 1 (or -1) state at the n+1 time step; then, the node B is now ready to operate the Turn Off rule to node A using the reciprocal connection. Such a dynamics can describe a feedback mechanism between nodes A and B, in which the first (de-)activates the second one, and the second turn off the first one sequencenly. However, the dynamical characterization for both synchronous or sequential mechanisms is not satisfied . Bottom panel Nodes operating with ActS tend to be reciprocal turned off by the TfO rule (on the left), as well as nodes reciprocally connected by ActO and TfS (on the right). Both cases do not satisfy the synchronization mechanism, but can be related by a sequence dynamics in which: first, node B (de-)activates node A, then node A turns off node B in the following step. Such a dynamics describes a possible feed-back mechanism between reciprocal nodes, whose schematic representation is shown. [Index]

Quantitative study of Centrality indexes
In the Centrality analysis [9] nodes can be sorted by specific features accounting for interactions within the whole system. The node degree is used to characterize node centrality and, in a directed and weighted network, the corresponding metrics. Thus, in-, out-and totalstrength indicate nodes endowed with a major in-Flux or out-Flux. In particular, we defined as central the nodes having a value higher than, at least, the third quartile of its distribution (>75%). In a multilayer description, besides the study of the node degree distribution in each layer, the distribution of central nodes across different layers was reckoned according to the following definitions: -Local central node: if present in not more than 25% of the total layers; -Intermediate multiplex central nodes: if present in the range between 25% and 75% of the total layers; -Multiplex central nodes: if present in more than 75% of total layers. In Figure S7  [Index]   Figure S8: Centrality Index of β and γ hierarchical levels. Higher hierarchical levels show similar pattern of central nodes as in the α-hierarchical level. The asterisk (*) indicates values higher than 85% of the distribution and double asterisk (**) higher than 95%. [Index]

Single links analysis
Tracing single connections points to describe the possible regulation pathways of the brain network. For the four layers of the α-hierarchical level the highest 40 links (corresponding to the 0.5% of the total) are identified and reported in Table S1 and Table S2 .
In the left side of Table 1, 26 ActS links (65% of the total) concern the same contralateral brain region: among these, eight are reciprocal and ten are not reciprocal with a trend in the left to right lateralization (7). The remaining links tend to fall in the same functional or anatomical module: -a) lateral frontal cortex (F3OP to F3T); -b) medial frontal cortex (ACIN and GR to F1M); -c) limbic (PHIP to HIP); -d) visual cortex (V1 to LING, O2 to O1, O3 and FUSI to O2, FUSI to O3); -e) temporal lobe (HES to T1, T1 to T1P); -f) non-specific pathways (HES to RO). Finally, most of the previous links (11/14) are homolateral (5/11 and 6/11 for right and left, respectively).
In summary, the ActS rule shows: 1) a connection to the same contralateral brain region mutual and non-mutual; 2) mutual bidirectional and non-mutual connections with non-homologous contralateral brain regions. Although some bidirectional connections were detected, the reciprocity analysis did not find any significant reciprocal connection. However, the ActS remains the only rule having a small trend in this direction, to be possibly ascribed to the cohesive nature of the rule. The non-homologous brain regions fall mainly in the same functional or structural network, e.g.: a) DMN with anterior cingulum to frontal medial/medial orbitalis, precuneus to posterior cingulum; b) cingulum cortex with anterior to middle cingulum; c) limbic with parahippocampus to hippocampus; d) vision with fusiform to calcarine, lingual and middle occipital gyrus. Thus, the above connections point to a dynamics falling within some specific functional modules of the brain regions.
The TfS links (Table S1, right side) enlighten the different pathways associated to central nodes: the basal ganglia, acting as out-strength nodes, take connections from some regions of cognitive and motor associated areas (PRE, SMA, RO, IN and T1), mainly towards the putamen (PUT). Thalamus and caudate point to other sparse regions: THA to F1M, PCIN and PQ; CAU to F2. Consequently, the same basal ganglia, acting as in-strength nodes, taking connection from: -limbic and DMN (OC, GR, F1M, F1MO, ACIN) to the caudate (CAU), -the lateral frontal (F2) and cingulate cortex (ACIN and MCIN) to the putamen (PUT). Other connections include: -temporo-frontal areas (F3O to T3), -limbic and DMN (PHIP to F1M and F1M to AMYG), -insula cortex (IN) to lateral frontal areas (F3O and F3OP), -occipital/visual cortex and DMN (Q to F1M, FUSI to PCIN, T3P to O3), -the interacion between a sparse set of brain regions (T3 to HIP, F1 to PCL).
Thus, the basal ganglia seem to be an intermediate station starting from the limbic and the DMN towards the salience and motor network.
As for the ActO rule, a particular connectivity pattern arises (Table S2, Table S1: ActS and TfS single links. 1) The ActS rule shows mutual interactions between interhemispheric (right and left) part of the same cerebral region, intrahemispheric interactions between different brain regions, and recursive pathway within the same functional module; 2) the TfS rule describes a basal ganglia (caudate and putamen) pathway toward frontal and DMN/limbic brain regions. Minimal fraction (0.5% = 40 connections) of the highest value weight links.
In bold central nodes of in-and out-strength. Asterisk (*) and double asterisk (**) indicate a strength of the distribution values higher than 85% and 95%, respectively.  1) The ActO rule shows a recursive interactions between the DMN/libic and fronto-parietal brain regions; 2) in TfO rule several pathways are detected: from a sparse set of cortical and sub-cortical brain areas to the fronto-parietal and DMN/limbic brain regions. Minimal fraction (0.5% = 40 connections) of the highest value weight links.

ActO
In bold central nodes of in-and out-strength. Asterisk (*) and double asterisk (**) indicate a strength of the distribution values higher than 85% and 95%, respectively.