causal strength of complex networks [58, 59]. However, the effects of excitations and inhibitions could not be differentiated in its original form. Based on the GC, we pro-pose a computational algorithm (presented in Section 3.2.3) under the assumption of best linear predictor (BLP) for analyzing neuronal networks by estimating the synaptic weights among them. The idea of the mathematical assumption BLP is that the weighted voltage-fluctuation of the pre-synaptic neurons should be the best linear explanation for the voltage-fluctuation of the post-synaptic neuron among the network. Using this inter-pretation, the GC analysis can be extended to measure both excitatory and inhibitory effects between neurons without too much extra computational complexity. The appro-priateness of the BLP assumption was examined by three sorts of simulated networks:
those with linear, almost linear, and nonlinear network structures. To illustrate the ap-plication of the proposed method, real spike trains from the anterior cingulate cortex (ACC) and the striatum (STR) were analyzed.
It is worth noting that spike trains are non-equally spaced data and are regarded as being from a point process. Filtering is usually required for converting them to equally spaced time series for further GC analyses [57]. This study adopted the Gaussian kernel filtering or binning (depending on the situation) to convert spike trains into time series data for the following three main reasons: (1) it reduces the complexity of analysis, and considers also the effect of temporal summation of action potentials, (2) under suitable preprocessing, even short, sparse spike trains can be converted, so that the standard autoregression modeling can be applied [75], (3) most important of all, spike trains can be filtered to form close approximations to the firing rates or the voltage-fluctuations of the underlying neurons [42].
The rest of this chapter is organized as follows. In Section 2, we briefly introduce the so-called Granger causality index, and then extend it to measure both excitatory and inhibitory effects between network nodes by using the BLP assumption. Section 3 presents three network models to ensure the appropriateness of the proposed algorithm.
In Section 4, we apply the algorithm to analyze real spike train data. Section 5 provides some discussion about the results obtained from Section 3–4, shortcomings of the method, and related future works.
3.2.1 An introduction to the GC
Let X = (x1, x2, . . . , xn) be a stationary n-dimensional time series process with zero mean. The p-th order linear autoregressive model for X is given by
x1t =
p
X
r=1
a1,1r x1t−r+ · · · +
p
X
r=1
a1,nr xnt−r+ 1t
x2t =
p
X
r=1
a2,1r x1t−r+ · · · +
p
X
r=1
a2,nr xnt−r+ 2t ...
xnt =
p
X
r=1
an,1r x1t−r+ · · · +
p
X
r=1
an,nr xnt−r+ nt,
(3.1)
where ai,jr is the projection coefficient from the i-th time series onto the j-th time series at time lag r, representing the coupling strength from node j to node i in the network. The residuals 1, 2, . . . , n are zero-mean uncorrelated white noises with covariance matrix Σ. The diagonal entries {Σii = V ar(i), i = 1, . . . , n} measure the accuracy of the autoregressive prediction to each node based on the information from time stamps t − 1 to t − p.
To see whether the information contained in time series xj is useful in explaining the state of time series xi, namely, the importance of node j to node i, we can exclude the time series variable xj from (4.7) to obtain a reduced1 (n − 1)-dimensional autoregressive model with residual series ηi,j of xi and corresponding prediction error Γjii = V ar(ηi,j).
Here Γjii measures the accuracy of the prediction of xi based only on the previous values in time series {x1, . . . , xj−1, xj+1, . . . , xn}. If Σii in (4.7) is significantly less than Γjii in the reduced model in some suitable statistical sense, then we say that xj Granger-cause xi. This causality can be quantified by the GC index from xj to xi formulated as:
Fj→i= ln Γjii
Σii. (3.2)
It is clear that Fj→i = 0 when Γjii = Σii, i.e., xj has no causal influence on xi, and Fj→i > 0 when xj Granger-cause xi. Notice that Fj→i is defined only for j 6= i and is always nonnegative, i.e., Σii is bounded above by Γjii, since the full model defined in (4.7) should fit the data better than the reduced model. Finally, we note that the GC index (4.8) is significant if the corresponding coefficients ai,jr are jointly significantly different from zero. This can be assessed via an F -test on the null hypothesis that ai,jr are zero [29, 60]. The projection coefficients ai,jr and prediction errors Σii can be obtained
1The i−th equation of the reduced model reads xit=
p
X
r=1
bi,1r x1t−r+· · ·+
p
X
r=1
bi,j−1r xj−1t−r+
p
X
r=1
bi,j+1r xj+1t−r+
· · · +
p
X
r=1
bi,nr xnt−r+ ηi,jt , where the b’s are the corresponding projection coefficients.
by solving the Yule-Walker equation [39] and an efficient model order can be determined using the Akaike Information Criterion (AIC):
AIC(p) = 2 log(det(Σ)) + 2n2p
T , (3.3)
where T is the total length of the time series. More information about the GC can be found in the author’s previous work [61] and also [4, 13, 20].
3.2.2 Synaptic weights estimation
Now, consider the multivariate zero-mean time series, X = (x1, x2, . . . , xn), consisting of the trajectories of membrane voltage of n distinct neurons. Suppose that the i-th neuron is triggered by some other k neurons in the network, say {i1, i2, . . . , ik}-th neurons, with synaptic weights {αi1, αi2, . . . , αki}. For convenience, we assume that 1 ≤ k := k(i) ≤ n−1.
The case k = 0 means that the i-th neuron is not triggered by others, thus is relatively easy to deal with. The weights are assumed to be nonzero, if some αsi is zero, then we can just remove the corresponding index is from the trigger set. Positive and negative weights represent excitatory and inhibitory influences on the i-th neuron, respectively.
In general, the trigger set Ii := {i1, i2, . . . , ik} and the corresponding weights αi :=
{αi1, αi2, . . . , αik} in the network can not be identified and estimated easily due to the un-derlying complex dynamics. However, under the assumption of the best linear predictor (BLP) (Definition 1 below), Ii and αi can be approximated effectively. The results are described in the following proposition.
Definition 1.
Let x and y be two stationary time series with zero-means. Then we say that y forms the best linear predictor (BLP) of x among a variable set G if σ2(x|¯x, ¯y) < σ2(x|¯x, ¯z), ∀z ∈ G, where σ2(x|¯x, ¯y) := minp,{fr},{dr}E{xt−Pp
r=1[fryt−r + drxt−r]}2. Proposition 1.
In the situation described above, if further the weighted trajectory ui := αi1xi1 + αi2xi2+
· · · + αikxik, made by the trigger set and the corresponding weights, forms the BLP to the trajectory of the i-th neuron (namely, xi) among the whole network; then based on the GC framework, Ii can be completely identified and the estimate of αi can be obtained as
ˆ
αi := {Pp
r=1ai,ir 1,Pp
r=1ai,ir 2, . . . ,Pp
r=1ai,ir k} up to a scale factor.
Proof:
Let ui := α1ixi1 + · · · + αikxik form the BLP of xi, then there exist a positive integer p and projection coefficients {fri, r = 1, 2, . . . , p}, {dir, r = 1, 2, . . . , p} such that xit = Pp
r=1[friuit−r+ dirxit−r] + it, where i is a stationary white noise possessing the smallest
variance among the whole network. Replacing ui with the weighted trajectory, we obtain
xit =
p
X
r=1
[friuit−r+ dirxit−r] + it
=
p
X
r=1
[fri(αi1xit−r1 + · · · + αikxit−rk ) + dirxit−r] + it
=
p
X
r=1
[αi1frixit−r1 + · · · + αikfrixit−rk + dirxit−r] + it, (3.4)
which represents the underlying but unknown network structure of {xi, xi1, xi2, . . . , xik}.
On the other hand, fitting to data the same equation as the i-th equation in (4.7), we have the following empirical regression (compared to the theoretical regression (3.4))
xit =
p
X
r=1
[ai,1r x1t−r + · · · + ai,nr xnt−r] + ˜it. (3.5)
Let ¯Ii := {1, 2, . . . , i − 1, i + 1, · · · , n} − Ii be the complement of the trigger set Ii. We note that ˜it ≡ it if ai,sr = 0, ∀r = 1, . . . , p and s ∈ ¯Ii; otherwise ˜i and i are totally different but with V ar( ˜it) = V ar(it) since (3.4) has the smallest residual variance among the whole network and (3.5) has more degree of freedom (coefficients) than (3.4).
If the trajectories of the ¯Ii-th neurons are stochastically independent of both the i-th and Ii-th neurons, then we have ai,sr = 0, for r = 1, . . . , p and s ∈ ¯Ii. Comparing (3.5) with (3.4), we then have
p
X
r=1
ai,ir 1 = αi1
p
X
r=1
fri , . . . ,
p
X
r=1
ai,ir k = αik
p
X
r=1
fri. (3.6)
Since the projection coefficients in (3.5) can be obtained by solving the Yule-Walker equation or simply by the least-squares method, (3.6) immediately leads to
αi1 p
X
r=1
ai,ir 1
= pαi2 X
r=1
ai,ir 2
= · · · = pαik X
r=1
ai,ir k ,
(3.7)
provided
p
X
r=1
fri 6= 0. For αis and
p
X
r=1
ai,ir s to have the same sign for s = 1, 2, . . . , k, we
can, without loss of generality, assume that
p
X
r=1
fri > 0. If
p
X
r=1
fri < 0, then −ui is used to replace ui.
If some of the trajectories of the ¯Ii-th neurons are linearly dependent of the i-th or Ii-th neurons, then ai,sr 6= 0, for some r ∈ {1, · · · , p}, s ∈ ¯Ii and the projection
coefficients in (3.5) are thus affected, resulting in a biased estimation of (3.7). However, this predicament can be solved by virtue of the assumption of BLP and the concept of GC. Since the i in (3.4) has the smallest variance among the whole network, taking out any element of {xs : s ∈ ¯Ii} from the regression (3.5) does not increase the variance of ˜i. According to the concept in (4.8), we can correct the model coefficients by ruling out all the useless information of the ¯Ii-th neurons.
We end this subsection by the following remarks.
Remark 1.
The idea behind the BLP mathematical assumption is that the weighted voltage-trajectory of the trigger neurons should be the best linear explanation for the voltage-trajectory of the target neuron among the whole network. Based on this interpretation, the GC index can be extended to measure both excitatory and inhibitory effects in virtue of the esti-mated synaptic weights.
Remark 2.
The synaptic weight αis and the summation of the projection coefficients
p
X
r=1
ai,ir s are
forced to have the same sign for s = 1, 2, . . . , k, because positive and negative
p
X
r=1
ai,ir s refer to positive and negative correlations between xi and xis, respectively.
Remark 3.
The case k = 0 means that the i-th neuron is not triggered by other neurons in the network, therefore Fj→i= 0, ∀j and there is no synaptic weights to be estimated.
Remark 4.
For readers dealing with sparse networks, L1 regularization (or LASSO) would be an use-ful technique for fitting to data a sparse regression to avoid overfitting and the problem of multiple testing [1,49]. In this scenario, a computational much more efficient approach would be first running LASSO to learn the network structure and then using GC to get the causal strength.
Remark 5.
Some arguments in the proof such as the stochastic independence and the estimations of projection coefficients assume that the law of large numbers (LLN) holds. In the case of small samples or limited data, estimation errors come into exsitence thus some statistical inferences in the proof may fail. However, the proof holds for most large-sample cases.
3.2.3 The algorithm
Here, we present a step by step algorithm for computing the proposed index (named neuron synaptic index, NSI) from multiple spike train data.
Step 1 : Properly smooth the spike train data by kernel filtering (Gaussian kernels are commonly used) to acquire the approximate membrane voltage trajectories to the underlying voltage evolution of the neurons in the network.
Step 2 : Subtract the mean value from each trajectory to form zero-mean time series and then fit the vector autoregressive model in (4.7). Appropriate model order can be obtained beforehand by using AIC in (4.9).
Step 3 : Compute all the GC indices by (4.8) for all pairs of neurons, i.e., i, j = 1, 2, . . . , n with i 6= j and also perform F -tests to ensure statistical significance.
Step 4 : For each node i = 1, 2, . . . , n, refine the autoregressive model by ruling out the information about the ¯Ii-th neurons, i.e., the neurons with insignificant Fj→i, to correct the projection coefficients.
Step 5 : For each node, compute the synaptic weights of the trigger neurons by simply summing the projection coefficients up to the model order p as shown in (3.7).
Step 6 : For each node, take the weighted trajectory as a new explanatory time series and then compute the GC index from this weighted time series to that of the node.
Step 7 : Finally, the NSI is then defined to be the l1-normalized 2 estimated weights obtained in Step 5 multiplying the GC index obtained in Step 6 (see (B.5) in Appendix).