Inferring functional connectivity through graphical directed information

Abstract

Objective. Accurate inference of functional connectivity is critical for understanding brain function. Previous methods have limited ability distinguishing between direct and indirect connections because of inadequate scaling with dimensionality. This poor scaling performance reduces the number of nodes that can be included in conditioning. Our goal was to provide a technique that scales better and thereby enables minimization of indirect connections.Approach. Our major contribution is a powerful model-free framework, graphical directed information (GDI), that enables pairwise directed functional connections to be conditioned on the activity of substantially more nodes in a network, producing a more accurate graph of functional connectivity that reduces indirect connections. The key technology enabling this advancement is a recent advance in the estimation of mutual information (MI), which relies on multilayer perceptrons and exploiting an alternative representation of the Kullback-Leibler divergence definition of MI. Our second major contribution is the application of this technique to both discretely valued and continuously valued time series.Main results. GDI correctly inferred the circuitry of arbitrary Gaussian, nonlinear, and conductance-based networks. Furthermore, GDI inferred many of the connections of a model of a central pattern generator circuit inAplysia, while also reducing many indirect connections.Significance. GDI is a general and model-free technique that can be used on a variety of scales and data types to provide accurate direct connectivity graphs and addresses the critical issue of indirect connections in neural data analysis.

Keywords: aplysia; causality; directed information; functional connectivity; indirect connectivity; information theory; mutual information.

Figure 1.

Examples of motifs addressed by graphical methods such as GDI. (a) In the case of node 1 directly influencing node 2 and node 2 directly influencing node 3, an indirect proxy connection from node 1 to node 3 will be detected by GDI and eliminated by conditioning on the activity of node 2. (b) In the case of an extended path from node 1 to node 4, simple pairwise analyses will catch two indirect connections. In contrast, GDI conditions analyses on all other nodes’ activity and therefore the indirect connections will be eliminated. (c) In the case of a sink, i.e. where multiple nodes influence one node, simple pairwise analyses may not detect direct connections because of thresholding methods that are often used in analyzing graphs. Because three connections to node 4 exist, each connection will appear to have a small influence on node 4 that may be sub-threshold. GDI’s conditioning on other edges will increase the relative contribution of a given direct edge, pushing it above threshold and allowing it to be detected.

Figure 2.

Overall process for estimating GDI between two nodes that relies on an MI estimation technique [14]. (a) Estimating GDI from X to Y includes conditioning on other nodes. For this visualization we specify two other nodes V and W. (b) First take windows of length M + 1 of nodes’ corresponding time series with M = 1 for this example. Gray points in windows are ignored because I_G(X → Y) does not include values at time i except for Y_i, which is because GDI is focused on causal rather than synchronous relationships. (c) Use values from windows as samples that capture joint statistics, and permute values to produce samples that capture independence. Permuted samples consist of windows from each time series that are no longer simultaneous, i.e. windows are effectively drawn independently from each time series in random orders. (d) Because GDI can be defined as the difference between two MI terms and therefore two KL-divergence terms, samples are used to estimate KL-divergences between the original data distributions and the shuffled data distributions. The extent to which a classifier can differentiate between the original and the shuffled data is quantified as MI, and the difference between the two MI terms provides an estimate of GDI.

Figure 3.

Continuously-valued example demonstrating the ability of graph theoretic methods to correctly differentiate excitatory and inhibitory connections. (a) A simple three-node system was analyzed, where X inhibits Y and W excites Y, which can be summarized as Y(t) = W(t − 0.125) − X(t − 0.125). (b) Snapshot of the sinusoidal waveforms generated. The slight dip around every peak of X combined with the typical peak of W produces the delayed but emphasized peak of Y 0.125 s later. If X did not exhibit this dip, it would be impossible to identify X as inhibitory because its activity would be fully masked by the excitatory effect of W. (c) Using the peak of the normalized cross-correlation (CC) R_XY or the time lagged correlation (TLC) ρ_XY results in incorrect identification of the connection from X to Y as excitatory. However, both methods do result in correct identification of the connection from W to Y as excitatory. (d) In order to perform partial inference, which accounts for the activity of other nodes in pairwise analysis, the residuals resulting from regressing X and Y against W (pictured) as well as from regressing W and Y against X are obtained. (e) Time lagged partial correlation (TLPC), i.e. the correlation coefficient between the residuals at different time lags, correctly identifies both the inhibitory and excitatory connections by accounting for other nodes’ activity.

Figure 4.

Discretely-valued example demonstrating the ability of graph theoretic methods to correctly differentiate excitatory and inhibitory connections. (a) A simple three-node system was analyzed, where X inhibits Y and W excites Y, both with a delay of approximately 0.125 s. (b) Snapshot of the spikes generated. The intermittent dropout of X allows for enhanced excitation of Y by W. (c) Both the normalized cross-correlation (CC) R_XY and time lagged correlation (TLC) ρ_XY fail to identify the connection from X to Y as inhibitory via use of peak values. However, both methods do correctly identify the connection from W to Y as excitatory. (d) In order to perform partial inference, which accounts for the activity of other nodes in pairwise analysis, the residuals resulting from regressing X (d1) and Y (d2) against W as well as from regressing W and Y against X are obtained. (e) Time lagged partial correlation (TLPC), i.e. the correlation coefficient between the residuals at different time lags, correctly identifies both the inhibitory and excitatory connections by accounting for other nodes’ activity.

Figure 5.

Scaling plots. (a) A Gaussian simulation was implemented where ρ = 0.6. (b) A binary symmetric channel simulation was implemented where X_i−1~Bernoulli(0.3), Pr[Y_i = X_i−1] = 0.9, and I(X→Y) is known analytically. (a1,b1) Analysis of how GDI scales with number of samples N. Error bars indicate mean and variance of GDI estimates across 50 independent simulations with 10 bootstrap iterations used per estimate. d_z indicates number of dimensions being conditioned, which in the top row corresponds to independent unit Gaussian random variables and in the bottom row corresponds to discrete uniform random variables (0 or 1). d_z is therefore proportional to the number of nodes being conditioned, and accordingly indicates how GDI scales with network size. (a2,b2) Analysis of how GDI scales with number of dimensions being conditioned for given sample sizes N. (a3,b3) Analysis of how GDI estimates converge as number of bootstrap iterations increases. Each line is the cumulative average GDI estimate over bootstrap iterations for an independent simulation with d_z = 0 and N = 2500 samples.

Figure 6.

Analysis of signed GDI performance on Gaussian network. (a) The network structure incorporates four motifs: an extended path where there is a trail of connections from node 6 to 11, a sink where 3 other nodes influence node 2, a multipath connection from node 6 to 10, and an isolated node. Each connection indicates that the most recent past value of the causal node influences the current value of the node with a circle (inhibited) or triangle (excited). (b1) Considering the 10 highest edges, signed non-GDI, i.e. DI without conditioning on other nodes, detects a number of indirect connections and therefore misses a number of direct connections. (b2) Ten highest signed GDI values are the true direct connections and do not include any indirect connections, which is achieved by its conditioning ability. (c) Comparing all estimated DI (c1), estimated GDI (c2), and analytic GDI (c3) values. DI misidentifies many false connections (c1), while GDI correctly recovers the graph (c2). Similarity between estimated GDI (c2) and analytic GDI (c3) highlights estimator accuracy.

Figure 7.

Analysis of GDI performance on non-Gaussian and nonlinear network. (a) The network structure is the same as that of the prior figure, incorporating four motifs: an extended path where there is a trail of connections from nodes 6 to 11, a sink where 3 other nodes influence node 2, a multipath connection from nodes 6 to 10, and an isolated node. Each connection indicates that the most recent past value of the causal node influences the current value of the node with a circle (inhibited) or triangle (excited). (b1) Considering the ten greatest edges, non-GDI misidentifies two indirect connections and therefore misses two direct connections. (b2) Considering the ten greatest edges, GDI correctly identifies the true direct connectivity structure and does not identify any indirect connections, which is achieved by its conditioning ability. (c) Heatmap of true connectivity structure. (d) Comparing all values of estimated DI (d1), which mistakenly identifies indirect connections, and GDI (d2), which correctly identifies the direct connections and reduces the indirect connections to negligible values.

Figure 8.

Analysis of signed GDI performance on network of conductance-based spiking neurons with same edge structure as prior two examples, however weights have been modified. (a) Structure of direct connectivity between neurons. Linewidths proportional to synaptic strength. (b) Raster plot of simulated spike times. (c) Heatmap of direct connectivity between neurons. Color intensity proportional to synaptic strength. (d1) DI detects the direct connections but also includes many false positives. (d2) GDI analysis minimizes indirect connections by conditioning pairwise analyses on other simulated neurons’ activity, better characterizing the true connectivity structure.

Figure 9.

Connectivity inference on simulated CPG network. (a) Raster plot of spike times resulting from simulation of CPG implemented in SNNAP. (b) Connectivity matrix for CPG network where values indicate the nA shift in rheobase. Positive values indicate a reduction in rheobase (excitatory connection) whereas negative values indicate an increase in rheobase (inhibitory). (c1) Non-graphical DI inference of CPG network connectivity. Many true connections are detected, however a variety of false positives are included particularly in the approximately top-left quadrant. (c2) GDI inference of CPG network connectivity. Many of the values representing false positives in the DI analysis are correctly reduced while the stronger true connections generally seem to be appropriately preserved.