Extending Transfer Entropy Improves Identification of Effective Connectivity in a Spiking Cortical Network Model

Transfer entropy (TE) is an information-theoretic measure which has received recent attention in neuroscience for its potential to identify effective connectivity between neurons. Calculating TE for large ensembles of spiking neurons is computationally intensive, and has caused most investigators to probe neural interactions at only a single time delay and at a message length of only a single time bin. This is problematic, as synaptic delays between cortical neurons, for example, range from one to tens of milliseconds. In addition, neurons produce bursts of spikes spanning multiple time bins. To address these issues, here we introduce a free software package that allows TE to be measured at multiple delays and message lengths. To assess performance, we applied these extensions of TE to a spiking cortical network model (Izhikevich, 2006) with known connectivity and a range of synaptic delays. For comparison, we also investigated single-delay TE, at a message length of one bin (D1TE), and cross-correlation (CC) methods. We found that D1TE could identify 36% of true connections when evaluated at a false positive rate of 1%. For extended versions of TE, this dramatically improved to 73% of true connections. In addition, the connections correctly identified by extended versions of TE accounted for 85% of the total synaptic weight in the network. Cross correlation methods generally performed more poorly than extended TE, but were useful when data length was short. A computational performance analysis demonstrated that the algorithm for extended TE, when used on currently available desktop computers, could extract effective connectivity from 1 hr recordings containing 200 neurons in ∼5 min. We conclude that extending TE to multiple delays and message lengths improves its ability to assess effective connectivity between spiking neurons. These extensions to TE soon could become practical tools for experimentalists who record hundreds of spiking neurons.

1. For all time bins t, count the number of distinct firing patterns at i t+1 , i t , and j t . This is done to estimate transition probabilities (e.g. p(i t+1 , i t , j t )). For first-order TE, there are 2 3 = 8 possible patterns because a neuron is either spiking or not in a given time bin.
2. Using the estimated transition probabilities, compute the TE from j ⇒ i (T ij ) by using Equation 1 below.
The first step is computationally expensive, so care must be taken to minimize the number of steps taken for each pair of neurons. One way to picture the algorithm is shown in Figure 1, where a copy of the time-series for neuron i is shifted backwards by one time bin (neuron i ) and stacked on top of the time-series for neurons i and j. This copying and shifting of i allows us to look at each column of the stacked series and simply count the distinct patterns observed 1 . If we visit every column, we will have counted up all firing patterns present in the data and can proceed to estimate the transition probabilities. At a high-level, the first-order algorithm looks something like Listing 1. For each pair of N neurons, we count up the firing patterns for all time bins t ∈ {1 . . . D}, where D is the duration of our data, and then use those counts to compute the TE from j ⇒ i. While this algorithm works fine for dense time-series with a large proportion of spikes to time bins, it is unnecessarily slow for sparse time series. When the proportion of spikes to time bins is low, the majority of firing patterns will be empty (i.e. no neurons are spiking). How can we use this to our advantage?
Listing 1: First-order transfer entropy algorithm for dense time-series Because most time bins will have no neurons spiking, we only need to visit the time bins where at least one of our series (i , i, and j) has a spike. The number unvisited time bins will then be the count of our "no spike" firing pattern. If we assume that our time-series are stored as a series of integers representing the time bins that a given neurons was spiking, the algorithm will look like Listing 2.
We use the current and move next functions to represent a pointer into each time series that can be moved forward. Instead of visiting every time bin, t now jumps around to just the non-empty bins by moving to the minimum current value of all time-series. For each visited bin, we count of the firing patterns as before, but we need to compute the number of unvisited time bins (empty firing pattern) before calculating the TE from j ⇒ i. For delayed TE, we would like to calculate the TE between i and j for multiple delays. This is easily accomplished with the algorithm above by shifting the time-series for neuron i (and also for i ) backwards by some delay d.

First-Order Algorithmic Complexity
In order to predict the calculation time for a given data set, we need to identify which aspects of the TE algorithm significantly affect performance. We can then derive a formula for estimating the amount of time it will take to run any TE calculation.
The number of neurons N will definitely affect performance, since this quantity is used in the two outer for loops in Listing 2 on lines 4 and 5. In fact, if we assume the body of the j loop (starting on line 6) will take c seconds each time to run, then calculating TE for one delay time will take roughly cN 2 seconds. In computer science terminology, this means the calculation time will grow quadratically with the number of neurons.
What other factors will affect performance? The number of time bins we have to visit will certainly have an impact. For dense time-series, this will be D time bins. Sparse time-series, however, we require that we go through and count how many bins were visited for all pairs of neurons. Unfortunately, doing this is almost as computationally expensive as calculating TE in the first place! If we know the average firing rate of our neurons F , though, we can estimate the number of visited time bins for any pair of neurons by calculating F D. When our neurons all have a fairly consistent firing rate, this estimate will work well. As D increases, then, we should expect the calculation time to (roughly) grow linearly due to F D being linear in D.
For first-order TE, estimating the transition probabilities (line 46 of Listing 2) for a pair of neurons does not take any longer for different data sizes. We can therefore bundle this and other machine-specific details into a constant C, making our final calculation time formula: where N is the number of neurons, D is the duration (number of time bins), F is the average firing rate, and C is a constant that depends on the actual machine and desired time units.

Higher-Order Algorithm
Our higher-order TE algorithm is a straightforward generalization of the first-order algorithm in Listing 2. We consider more time bins in either the receiving or sending neuron (or both) by stacking additional, shifted copies of their respective time-series (see Figure 2). As before, we then count up the distinct firing patterns by looking at the columns of our stacked series. Estimating the transition probabilities now requires that we take the additional time bins into account, resulting in Equation 3.
where k is the order of the receiving neuron i and l is the order of the sending neuron j. The pseudo-code in Listing 3 gives a bird's eye view of the higher-order algorithm. The quantity R represents the total order of the calculation, which is k + l + 1. For first-order TE, R = 3 and we only had three time-series to deal with. For higher-order TE, ij series will be of length R and the number of distinct firing patterns will be 2 R . so. This is because the actual run time of one iteration of the loop on line 47 is so small doubling it a handful of times makes little difference until R gets large.
Putting it all together, we get the final TE calculation time formula in Equation 4. F DR and 2 R are added together because the 2 R loop is run after the F DR loop (instead of inside it). C is still a constant that represents machine-specific factors and the desired time units.