Joint Statistics of Strongly Correlated Neurons via Dimensional Reduction
JJoint Statistics of Strongly Correlated Neurons via DimensionalReduction
Ta¸skın Deniz & Stefan Rotter
Bernstein Center Freiburg & Faculty of Biology,University of Freiburg, Hansastraße 9a, 79104 Freiburg, Germany (Dated: November 20, 2018)
Abstract
The relative timing of action potentials in neurons recorded from local cortical networks oftenshows a non-trivial dependence, which is then quantified by cross-correlation functions. Theoreticalmodels emphasize that such spike train correlations are an inevitable consequence of two neuronsbeing part of the same network and sharing some synaptic input. For non-linear neuron mod-els, however, explicit correlation functions are difficult to compute analytically, and perturbativemethods work only for weak shared input. In order to treat strong correlations, we suggest here analternative non-perturbative method. Specifically, we study the case of two leaky integrate-and-fireneurons with strong shared input. Correlation functions derived from simulated spike trains fitour theoretical predictions very accurately. Using our method, we computed the non-linear corre-lation transfer as well as correlation functions that are asymmetric due to inhomogeneous intrinsicparameters or unequal input. a r X i v : . [ q - b i o . N C ] S e p . INTRODUCTION Electric activity generated by different neurons in the brain is often strongly correlated[1–4]. These correlations arise as a result of shared input, or input components that arethemselves correlated. Correlated activity can be a consequence of shared background fluc-tuations [5], but strong correlations might also indicate synchronous action potentials atthe input indicating temporal coding. However, a clear-cut dichotomy between decorrelatedand synchronized dynamics does not exist [6] [7, 8]. Rather, one should consider these twoextremes as two faces of the same coin. Recent high-precision measurements reported verylow average correlations suggesting a mechanism of active decorrelation in cortical networks[9–12]. At the same time it was observed by intracellular measurements that nearby neuronsreceive very similar input [2–4].Several studies of pair correlations in neural networks relate structure and dynamicsassuming a fluctuating dynamics about a fixed point that is characterized by asynchronous(A) population activity and irregular (I) spike trains [11, 13–15]. They employ essentiallylinear perturbation theory [16, 17] to compute correlation functions. Nevertheless, someof these works push the limits of existing methods. First of all, a qualitatively differentAI state was observed in simulations of spiking neural networks with stronger couplings[18]. Secondly, a partial extension of the theory to the strongly correlated regime was basedon numerically determined spike response functions [14]. Thirdly, pair correlation studieswere generalized to higher-order correlations in recurrent networks by accounting for certainnetwork connectivity motifs [19 ? ]. These studies exploit and extend existing methods, butthey also demonstrate the need for a new approach.The main assumption of perturbation theory is that the common drive of the two neuronsis weak. Yet, this criterion depends on the background state and strength of interactionsin a given network. We showed previously that low background rates, for example, maylead to a breakdown of perturbation theory even for low correlation coefficients c [20]. Thismakes non-perturbative methods indispensable for modeling and analysis of correlationsas low spike rates are typical in experiments [21]. All in all, a proper treatment of strongcorrelations must take the non-linear correlation transfer into account, which appears to playan important role in sensory processing [22]. However, a unified and transparent frameworkto calculate correlations of all strengths for neuron models of the integrate-and-fire type does2ot exist.The pitfall of previously suggested non-perturbative methods are their immense compu-tational costs due to the high dimensionality of the discretized problem [23]. This makescomputations practically impossible for a large range of parameters. For instance, the nu-merical effort of computing the pair correlation problem scales as N , where N is the numberof grid points used to approximate single neuron membrane potentials. Although, limitingthe grid size is possible [24], a too coarse voltage grid fails to properly reflect the statisticsof leaky-integrate-and-fire neurons with Poisson input. The precision issue gets even moresevere for correlations of higher order, which are needed to parametrize the joint statistics ofmultiple neurons. With our methods, in contrast, we observed that joint membrane potentialdistributions of even strongly correlated neurons can be reduced to a small set of principalvectors via singular value decomposition (SVD). This suggests that strong correlations canbe computed with high precision resorting to subspaces of relatively low dimension. In thiswork, specifically, we devise a SVD based method that allows to compute spike correlationfunctions of two leaky integrate-and-fire neurons with high accuracy.Similar problems were studied analytically for arbitrary input correlations of the stochas-tic dynamics of neural oscillators [25], and for level-crossings of correlated Gaussian processes[26]. Related numerical work considered strong input correlations for integrate-and-fire neu-rons receiving white noise input [27] or shot noise input with nontrivial temporal correlations[28, 29]. The problem of analytically calculating the stationary distributions conditional ona spike from the exit current at the threshold is also discussed in the case of colored noise[28]. A method to deal with very strong input correlations ( c ≈
1) in a specific input model(correlated Poisson processes) was suggested by [30]. Our study further suggests a noveltechnique, extending [24], to solve 2D jump equations for leaky integrate-and-fire neuronsby mapping it to a Markov chain. This method provides an accurate estimate of the steadystate joint distribution of membrane potentials.We test our method for large input correlations c and demonstrate its power for differenttypes of correlation asymmetries by comparing our semi-analytical approach to correlationsextracted from simulated neuronal spike trains. We look at the full range of input correla-tions and provide an example of a non-linear correlation transfer function. Our method canbe extended to more general integrate-and-fire models, higher-order input correlations andthird order output correlations. However, we have to defer a detailed analysis of such cases3o future work. +++ ab c FIG. 1. (a) Two LIF neurons receiving private and shared shared inputs, both represented byexcitatory and inhibitory spike trains. (b) Two neurons with shared white noise input, with firstand second moments matched to (a). (c) Schematic describing the threshold crossing and resetmechanism that is part of the membrane potential dynamics of LIF neurons.
II. MODELS & METHODSA. Two neurons with shared Poisson input
The leaky inteagrate-and-fire neuron model with postsynaptic potentials of finite ampli-tudes was studied previously in [24]. The stochastic equation that describes the membranepotential dynamics of one particular neuron is given as τ m ˙ V ( t ) = − V ( t ) + h ex τ (cid:88) i S exi ( t ) + h in τ m (cid:88) j S inj ( t ) . (1)where h ex and h in represent the amplitudes of individual EPSPs and IPSPs, respectively,and S exi ( t ) and S inj ( t ) are the spike trains of excitatory and inhibitory presynaptic neurons,with each of their spikes represented by a Dirac delta function S exi ( t ) = (cid:88) n δ ( t − t in ) . (2)4n analogous definition holds for inhibitory neurons. In both cases, if the membrane poten-tial reaches the firing threshold, V th , a spike is elicited and the voltage is reset to its restingvalue at 0.In order to study correlations between the spike trains of two neurons we look at twocoupled stochastic equations, describing the membrane potentials of two neurons that sharea certain fraction of their excitatory and inhibitory input spikes˙ V ( t ) = − V ( t ) τ m, + h ex ( (cid:88) k S ex ,k ( t ) + (cid:88) m S exm ( t )) + h in ( (cid:88) l S in ,l ( t ) + (cid:88) n S inn ( t )) (3)˙ V ( t ) = − V ( t ) τ m, + h ex ( (cid:88) k S ex ,k ( t ) + (cid:88) m S exm ( t )) + h in ( (cid:88) l S in ,l ( t ) + (cid:88) n S inn ( t )) (4)where shared excitatory and inhibitory input spike trains are denoted as (cid:80) m S exs,m ( t ) and (cid:80) n S ins,n ( t ), respectively. Comparing the parameters of the jump model to shared white noiseinput we implicitly specified the firing rates of 6 independent Poisson processes (Fig. 1)˙ V ( t ) = − V ( t ) τ m, + µ + σ √ τ m, ξ ( t ) + σ c, √ τ m, ξ c ( t ) (5)˙ V ( t ) = − V ( t ) τ m, + µ + σ √ τ m, ξ ( t ) + σ c, √ τ m, ξ c ( t ) . (6)For notational convenience, shared input rates r ex,c and r inh,c are computed from a sharedWiener process ξ c with zero mean. Setting h ex = h and h in = gh , rates for each independentexcitatory and inhibitory process are given as r ex = 11 + g (cid:16) σ τ h + gµhτ (cid:17) (7) r in = 1( g + g ) (cid:16) σ τ h − µhτ (cid:17) . (8)Here we note that some combinations of ( µ, σ ) do not correspond to any combination ofpositive Poisson rates, as they must satisfy σ ≥ hµ (9)guarantee r in ≥ B. Discretized Markov operators for the LIF model
In this section we summarize the discrete approximation to the dynamics of a LIF neuron,as developed in [24]. The coarse graining of the membrane potential is given by the map R → Z , V (cid:55)→ (cid:98) V ∆ V (cid:99) . (10)5 ernel Density Estimation via Di↵usion Taskin DenizBCF FreiburgJuly 22, 2016
Kernel Density Estimation via Di↵usion
Taskin DenizBCF FreiburgJuly 22, 2016
Kernel Density Estimation via Di↵usion
Taskin DenizBCF FreiburgJuly 22, 2016
Kernel Density Estimation via Di↵usion
Taskin DenizBCF FreiburgJuly 22, 2016 LR T FIG. 2. Schematic illustration of Singular Value Decomposition in 2 dimensions, M = L Σ R T , forpositive definite matrices (det( M ) > L and R are orthogonal matrices, and D is a real diagonalmatrix. The non-negative diagonal elements σ and σ of the matrix Σ are the so-called singularvalues of the matrix M . The probability density function P ( v ) of the membrane potential then becomes a vector p = ( p i ) satisfying p i = (cid:90) v i +1 v i P ( v ) dv, (11)with v i = i ∆ V . As we impose a cut-off lower boundary V − of the voltage scale, the dimensionof the discrete state space is given as N = V th − V − ∆ V .The temporal evolution of the membrane potential distribution is now described in termsof a Markov process. The Markov propagator can be expressed as a juxtaposition of threeoperators: a decay operator D describing leaky integration, a jump operator J accountingfor the action of synaptic inputs, and a threshold-and-reset operator T that implementsspike generation upon threshold crossing. The discrete decay operator D is derived from itscontinuous counterpart D as follows∆ V D (cid:126)p = (cid:90) v i +1 v i DP ( x ) dx = (cid:90) v i +1 v i e ∆ t/τ m (cid:124) (cid:123)(cid:122) (cid:125) q P ( e ∆ t/τ m x ) dx = (cid:100) q ( i +1) (cid:101)− (cid:88) (cid:98) qi (cid:99) p j . ∆ V + (cid:16) ∆ V (cid:100) qi (cid:101) − qv i (cid:17) p (cid:98) q i (cid:99) + (cid:16) qv i +1 − ∆ V (cid:98) q ( i + 1) (cid:99) (cid:17) p (cid:98) q ( i +1) (cid:99) . D ij = k = (cid:100) q ( i +1) (cid:101)− (cid:88) (cid:98) qi (cid:99) δ j,k + (cid:16) (cid:100) qi (cid:101) − qv i (cid:17) δ j, (cid:98) q i (cid:99) + (cid:16) qv i +1 − (cid:98) q ( i + 1) (cid:99) (cid:17) δ j, (cid:98) q ( i +1) (cid:99) . (12)The jump distribution, which underlies the jump operator J , can be derived from thecount distribution of excitatory and inhibitory synaptic events in a small time interval ∆ t .Assuming that they arrive with Poisson statistics (maximum entropy) at a rate of r ex and r in , respectively, we obtain P ( m, n ) = a m b n m ! n ! e − ( a + b ) (13)where a = r ex ∆ t and b = r in ∆ t are the corresponding expected event counts. Dummyindices m and n correspond to excitatory and inhibitory counts. The jump distribution ofthe membrane potential γ = ( m − gn ) h is given by P ( γ ) = ∞ (cid:88) m,n =0 P ( m, n ) δ γ, ( m − gn ) h . (14)The jump operator J is then derived as( J p ) i = (cid:88) γ P ( γ ) p i − γ ∆ V with a jump matrix given by J ij = (cid:88) γ P ( γ ) δ j,i − γ ∆ V . (15)The threshold-and-reset operator T takes all the states above threshold and inserts them atthe reset potential. This is simply given as T ij = I j> Vth ∆ V δ i, Vr ∆ V + I j ≤ Vth ∆ V δ i,j (16)Finally, the time evolution matrix of the corresponding Markov chain is given as a productof the three operators described above M = T J D . (17)The discrete stationary distribution is M P = P , (18)and the corresponding stationary rate is given as r Markov = ( τ ref + [ 1 h (cid:88) i ≥ i th ( J D P ) i ] − ) − . (19)7 . Spike train correlations The cross-covariance function of two stationary spike trains S a ( t ) = (cid:80) l δ ( t − t al ) ( a = 1 , C ( τ ) = (cid:104) S ( t + τ ) S ( t ) (cid:105) − (cid:104) S ( t + τ ) (cid:105)(cid:104) S ( t ) (cid:105) (20)where (cid:104) S a ( t ) (cid:105) = r a , with (cid:104) . (cid:105) indicating the ensemble average. For the model we studiedhere stationary rates were computed using Eq. 21. The cross-covariance function can beexpressed in terms of the conditional firing rate r | ( τ ), two stationary rates r and r , andthe amplitude r of a δ -function C ( τ ) = r δ ( τ ) + r ( r | ( τ ) − r ) . (21)Given the spiking neuron model considered here, Eq. 21 can be derived from the stationaryjoint membrane potential distribution P ( V , V ). Conditional on a spike at time t in thefirst neuron, one has to find the instantaneous distribution of the membrane potential of thesecond neuron P i | j ( V i ) = P ( V i | neuron j spikes at t ) . (22)The probability of observing a consecutive spike is given by the flux P flux ( V , θ ) with thenormalization in Eq. 23. The conditional flux is computed for the Markov approximationas in Eq. 47. In general we simply compute the exit rate at the threshold distributed over V by solving the following initial value problem P (0 , V ) = P | ( V ) = P flux ( V , θ ) (cid:16) (cid:90) θ ∞ P flux ( x, θ ) dx (cid:17) − (23)∆ t P ( t ) = M P ( t ) (24)where M is discrete time evolution matrix of the neuron model given by Eq. 17. M leadsto a (forward) time evolution in steps of ∆ t . The instantaneous conditional rate r | ( t ) inEq. 21 is computed with a formula similar to Eq. 19 r | ( t ) = r Markov ( t ) = ( τ ref + [ 1 h (cid:88) i ≥ i th ( J D P ( t )) i ] − ) − . (25)Finally, the covariance function C ( τ ) in Eq. 21 is derived by using r | ( t ) and two stationaryrates r and r and r . Note that the method described here can be generalized to higherorder correlations as well. 8 . Correlated jump distribution in 2D We now use a 2-dimensional state space describing the joint membrane potential evolutionof two neurons. Correlated and uncorrelated Poisson jumps push the 2D membrane potentialvector ( V , V ) into three independent directions, (1 , ,
1) and (1 , U , U ) to a new position ( V , V ) in state space drivenby 6 independent Poisson processes is obtained from a 2D convolution P ( V , V | U , U ) = (cid:90) P ( V − U − Z ) P ( V − U − Z ) P c ( Z ) dZ. (26)Inserting the expressions for the jump distributions in each direction, collecting all termsand using Eq. 13, we obtain P ( V , V ) = (cid:90) ∞ (cid:88) i,j =0 ∞ (cid:88) k,l =0 ∞ (cid:88) m,n =0 P ( i, j ) × P ( k, l ) P c ( m, n ) δ x − Z, ( i − gj ) h δ y − Z, ( k − gl ) h δ Z, ( m − gn ) h dZ. This expression is also valid for more general input statistic models that rely on a decompo-sition into statistically independent parts [30, 31]. Here we consider the shared Poisson inputmodel as described by Eq. 4. Integrating the expression with respect to Z and insertingthe mean event counts a i = r ex,i ∆ t and b i = r in,i ∆ t , the resulting expression is given in acompact form P ( V , V ) = ∞ (cid:88) i,j =0 ∞ (cid:88) k,l =0 ∞ (cid:88) m,n =0 e − ( a + b + a + b + a s + b s ) × a i b j i ! j ! a k b l k ! l ! a ms b ns m ! n ! δ V , ( i − m − g ( j − n )) h δ V , ( k − m − g ( l − n )) h . We can simplify this expression by choosing a regular grid according to h e = n ∆ V and h i = m ∆ V , for integers m and n . In practice, however, the resulting sum will be truncatedfor a given tolerance. We will use the matrix form of the discretized operators in the followingsubsection, which is equivalent to the above expression. E. Linear operators for correlated dynamics in 2D
Here we discuss the action of operators on state vectors of the discretized ( V , V ) space,assuming N bins in each dimension. We write the stationary voltage distribution in the9orm P ( X, Y ) ≡ N (cid:88) i,j =1 Ω ij X i ⊗ Y j . (27)where X and Y are two suitable orthogonal bases of R N . We will define a specific choicefor X and Y later in Section II F. As there is no crosstalk between the two neurons exceptby shared input leading to a correlated jump distribution, threshold and decay operators(tensors) are separable D D = D ⊗ D (28) T D = T ⊗ T . (29)Separability would also apply to the jump distribution for an input correlation coefficient of c = 0, corresponding to two neurons without shared input J D = J ⊗ J . (30)However, in the case of non-zero correlation, this relation holds only for a single path amongthe many connecting two points in state space . Therefore, every path must be taken intoaccount by considering the contribution of each operator to a movement in the oblique (1 , M D = T D J D D D (31)one can also find the stationary joint membrane potential distribution as the Perron-Frobenius eigenvector P of the propagator matrix M D M D P = P . (32)Regarding the correlated jump distribution there are two ways of constructing J operators.One possibility is described in Eq. 26. Alternatively, we construct linear Markov jumpoperators in 2D exploiting the commutativity of infinitesimal operators J D = e ( J p ⊗ I + I ⊗ J p + J c ⊗ J c ) . Here I is the identity matrix. Using the properties of the operator product ⊗ and commu-tativity of the individual factors, we can simplify this expression J D = ( e J p ⊗ I )( I ⊗ e J p )( e J c ⊗ J c ) .
10n order to expand the third term we define U and D operators as up and down transitionmatrices, where U corresponds to a 1-step up transition, and D corresponds to a 1-stepdown transition. This leads to U kij = δ i,j − k and D lij = δ i,j + l . (33)Hence, discrete approximations to infinitesimal generators of private components are givenas J i = a i U k + b i D l . (34)where matrix powers k and l are defined on h e = k ∆ V and h i = l ∆ V for simplicity.(Thisrestrictive assumption can be generalized easily by computing the continuous jump distribu-tion in Eq. 14 and then discretizing it, which leads to the same result.) On the other hand,we need to be careful with correlated spikes, which trigger jumps into the oblique direction(1 ,
1) with probability c mn = e − ( a c + b c ) m ! 1 n ! a mc b nc (35)for m excitatory and n inhibitory jumps. Expanding yields J D = (cid:88) m =0 (cid:88) n =0 c mn [ e J U km D ln ⊗ e J U km D ln ] . (36)As we noted before, the above construction is quite general and can be applied easily forgeneral amplitude distributions. We only need to consider a discrete amplitude distribution c mn in a given time bin of size ∆ t , as described above, see also [31, 32].A final remark on the method described in this section concerns the commutativity ofmatrices. This property leads to a numerically more economic expression J D = (cid:88) j ∈ J c ( j )[ e J O j ⊗ e J O j ] (37)where the set of integers J is defined as list of all jump numbers j = mk − nl . The coefficient c ( j ) = (cid:88) m,n e − ( a c + b c ) c mn δ j,mk − nl (38)is the probability of j jumps. The jump generator O is then defined in terms of matrixpowers O j = U j , j > I, j = 0 D j j < . (39)11 . Operator decomposition and SVD basis reduction The expansion method described above is straightforward, but rather cumbersome toimplement. We will now introduce a basis for the expansion of correlated jump operatorssuitable to reduce the cost of the computations involved, and helpful to increase the accuracyof a truncation. With our method, as compared to others, we have to solve linear equationsof lower dimensionality in order to get a better approximation for the correlation function.Some further arguments for selecting Singular Value Decomposition are discussed in theresults section. SVD of a matrix is given as M = L Σ R T (40)where the diagonal entries of the diagonal matrix Σ are the square roots of the non-zeroeigenvalues of MM T and M T M . Both matrices R and L are orthogonal with columnsconsisting of the eigenvectors of MM T and M T M , respectively. We show a 2D exampleSVD of Markov matrix in Fig. 2. For M replaced by a single-neuron time evolution matrix M ( M ), we define the matrix X ( Y ) by the selected orthogonal subspace of dimension K ( L ) of R ( R ), according to the largest K ( L ) singular values, resepectively. In order toproject J ( J )and J ( T ) to this subspace, we also extend the orthogonal basis X and Y to supra-threshold transitions˜ X = X I m , ˜ Y = Y I n , (41)where X is an M × K and Y is an N × L operator, respectively. I m and I n are identitymatrices, where m and n are the maximal supra-threshold jump numbers induced by bothprivate and shared inputs. The combined action of J and T for c = 0 is then expressed asreduced operators T J → X T T ˜ X ˜ X T J X, T J → Y T T ˜ Y ˜ Y T J Y, (42)which map M × M ( N × N ) to K × K ( L × L ) matrices. Below we use the same dimensionalreduction for correlated operators. In Eq. 37 the correlated jump operators are expressed as J D = (cid:88) j ∈ J c ( j ) A j ⊗ B j .
12 dimensional reduction is then achieved by using Eq. 42 in M D P = (cid:88) j ∈ J c ( j ) ( T A j D ⊗ T B j D ) P . In order to find P in Eq. 27 we need to solve Eq. 32 , which reads M D P = P . The projection operators X T ⊗ Y T satisfy ( X ⊗ Y )( X T ⊗ Y T ) P = P . Applying them tothe left hand side of this equation Eq. 32, we obtain( X T ⊗ Y T ) M D = ( X T ⊗ Y T ) M D ( X ⊗ Y )( X T ⊗ Y T ) P = (cid:88) j ∈ J c ( j ) (cid:16) ( X T T ˜ X ˜ X T A j XX T D X ) ⊗ ( Y T T ˜ Y ˜ Y T B j Y D Y ) (cid:17) Ω RHS = Ω . Ω ≡ ( X T ⊗ Y T ) P can be expressed in a simpler form as( X T ⊗ Y T ) P = ( X T ⊗ Y T )( (cid:88) ij Ω ij X i ⊗ Y j )= (cid:88) ij Ω ij X T X i ⊗ Y T Y j )= (cid:88) ij Ω ij e i ⊗ e j ≡ Ωfor monomials ( e i ) k = δ i,k . The reduced equation is then Q Ω = Ω (43)where Q is a tensor defined as Q = (cid:88) j ∈ J c ( j )( X T T ˜ X ˜ X T A j XX T D X ) ⊗ ( Y T T ˜ Y ˜ Y T B j Y D Y T Y ) . (44)The dimensionally reduced problem in Eq. 43 can then be solved in practice by reindexingtensor indices ( i, j, k, l ) (cid:55)→ ( I, K ) .
G. Conditional flux distribution
Using the decomposed 2D stationary distribution obtained by reduction, one can computethe flux distribution with the help of matrix operators. We compute the flux distribution13 b FIG. 3. Direct projections suffer from the imposed lower boundary and diverging dual eigenvectors.Therefore, we cannot increase the precision of our method using direct projections, as demonstratedabove (a) we show components log (Ω ij / max( S )) with N = 30 eigenvectors via dual spaceprojections. (b) Same as (a) with N = 80 eigenvectors. We observed that the example in (b) failsto converge as its maximum value is S , , because of numerical instabilities. using the 2D decay and jump operators with thresholding imposed only at one of the bound-aries. A scheme illustrating the situation is shown in Fig. 1c. Here we explain how to obtainthe flux distribution conditional on a spike in one neuron. In the general case of correlatedneurons, the action of the full operators J D and D D is given as˜ P J ,kl = J D D D P = (cid:88) k (cid:48) ,l (cid:48) (cid:88) j ∈ J c ( j )( A j D ⊗ B j D ) kl,k (cid:48) l (cid:48) ( X Tk (cid:48) Ω Y l (cid:48) ) , (45)where ˜ P J is a ( M + m ) × ( N + n ) matrix. The implicitly summed subspace indices, P ,kl = (cid:80) m,n X Tk (cid:48) ,m Ω mn Y n,l (cid:48) , are not shown and, A j , B j and c ( j ) are defined in Eq. 37. This equationcan be written in a concise form with implicit indices as˜ P J = (cid:88) j ∈ J c ( j ) A Tj D T X T Ω Y D B j = ˜ X T ˜Ω ˜ Y .
Again, for practical reasons, computations were reduced by projecting J D D D to extendedsubspaces ˜ X and ˜ Y similar to Eq. 42˜Ω = (cid:88) j ∈ J c ( j ) (cid:104) ( ˜ X T A j XX T D X ) ⊗ ( ˜ Y T B j Y D Y ) (cid:105) Ω . (46)In order to compute the probability of jumps we need to sum the probabilities for a jumpover the threshold V > θ or V > θ . The flux distribution is obtained as P flux ,k ∝ N + n (cid:88) l = N ˜ P J ,kl = N + n (cid:88) l = N ˜ X Tk ˜Ω ˜ Y l , P flux ,l ∝ M + m (cid:88) k = M ˜ P J ,kl = M + m (cid:88) k = M ˜ X Tk ˜Ω ˜ Y l (47)14efined for k < M and l < N . These expressions are normalized such that (cid:80) k P flux ,k = 1.The amplitude of the delta singularity at the reset potential is obtained as r = r ex,c N + n (cid:88) l = N M + m (cid:88) k = M ˜ X Tk ˜Ω ˜ Y l , (48)where r ex,c is the rate of excitatory shared spike trains. Once we have found the condi-tional flux distribution, we can solve the initial value problem defined in Eq. 23 in order todetermine the conditional rates r | ( t ) or r | ( t ). ac b FIG. 4. Reverse engineering of a known stationary distribution. (a) Projections of a knownstationary distribution (obtained by the full jump distribution as in 26) P = L Σ R T on the rightsingular vectors and (b) on the left singular vectors of the single-neuron time evolution matrix M .The result is shown here for symmetric neuron parameters L = R = W . (c) Reconstruction of astationary 2D joint membrane potential distribution. Singular vectors sorted by decreasing singularvalues and added one by one P ( K ) = (cid:80) Kk =1 σ k W k ⊗ W k , increasing the number of components fromtop left to bottom right. The convergence is relatively fast despite the rather high correlation of c = 0 .
7. Note that we used here a coarse grid ∆ V = 0 . O ( N ) operations. . Diffusion approximation vs. finite PSPs We compare the exit rate of the stochastic system with post-synaptic potentials of finiteamplitude with the analytic result obtained for the diffusion approximation [ ? ]. For smallenough PSPs the difference in rates of the two models is small r − sg = τ ref + τ m √ π (cid:90) Vth − µσVr − µσ e x [1 + erf( x )] dx (49) r − = τ ref + [ 1 h (cid:88) i ≥ i th ( J D p ) i ] − (50)We use the absolute difference between the two rates r error ( µ, σ ) = | r sg ( µ, σ ) − r Markov ( µ, σ ) | (51)to account for the accuracy of a specific space-time grid (∆ V, ∆ t ). I. Correlation coefficient and comparison to diffusion approximation
The correlation coefficient in 10a is computed with the formula C out ( c ) = r + (cid:82) r ( r | ( τ ) − r ) dt + (cid:82) r ( r | ( τ ) − r ) CV CV √ r r dt. (52)We used 49 for the stationary rate r , r is the amplitude of the δ -function in Eq. 48, andthe coefficient of variation, CV = σ ISI µ ISI , is computed with the equation CV = 2 πr (cid:90) Vth − µσVr − µσ e x dx (cid:90) y −∞ [1 + erf( x )] dy (53)as given in [13]. Thereby erf( x ) is the error function [33]. We computed correlation coeffi-cients of the diffusion approximation of the finite PSP system (meaning a 2D Fokker-Planckequation with the same c and σ s) in [20]. J. Direct numerical simulations and data smoothing
We used the neural simulation tool NEST [34] to perform numerical simulations of inputand output spike trains in the scenario described above. All analyses were based on dis-cretized voltage data obtained during simulations of 1 000 s duration using a time resolutionof ∆ t = 10 − s. 16mpirical voltage distributions were obtained by normalizing histograms appropriately.Further smoothing using a simple moving average was performed before comparing thesedistributions to the analytically obtained stationary distribution. We also performed thecomparison using cumulative distributions, as the implicit integration very efficiently reducesthe noise in the data. Two 2D distributions are compared via visual inspection of contourlines. We also directly compare spike train cross-correlation functions to assess efficiencyand accuracy of the method. K. Numerical evaluation of cross-correlation functions
We compute cross-correlation functions of spike trains from conditional PSTHs. One canexpress this as an integral over two variables τ = t − t and s = t + t with bin size ∆ C ( τ ) = 1∆ (cid:90) τ +∆ τ dτ (cid:48) u ( τ (cid:48) ) − l ( τ (cid:48) ) (cid:90) u ( τ (cid:48) ) l ( τ (cid:48) ) (cid:88) i,j δ ( τ (cid:48) − τ i ) δ ( s (cid:48) − s j ) ds (cid:48) (54)where we set u ( τ ) = T / √ | τ | , l ( τ ) = T / √ − | τ | . (55)with observation window T . L. Convergence and error bounds
The direct singular value decomposition of a 2D membrane potential distribution showsthat there are only few singular values that deviate significantly from 0 (Fig. 4). Thisbehavior does not depend strongly on the chosen discretization, but it does depend on theinput correlation coefficient c . Although singular vectors are not probability distributionsin their own respect, the singular vectors X i derived from the neuronal dynamics matrix(except the first few vectors) have the property (cid:88) k X ik ≈ (cid:88) ij P ij = (cid:88) i,j K,L (cid:88) k =0 ,l =0 X ik Ω kl Y jl = Σ S Σ (57)17s progressively small, where Σ k and Σ l are sums of k -th and l -th singular vectors of thefirst and the second matrix for k ≥ m and l ≥ n , respectivelyΣ err ( m, n ) = (cid:88) i,j K,L (cid:88) k = m,l = n X ik Ω kl Y jl = K,L (cid:88) k = m,l = n Σ k Ω kl Σ l . (58)This shows that the sum converges rather quickly. This error measure is related to projec-tions of 1D discrete stationary distribution P (satisfying M P = P ) to SVDs. All othereigenvectors of a Markovian matrix ( M P = λP for | λ | <
1) satisfy (cid:80) k ( P i ) k = 0. Wewant to avoid underestimating the total probability mass as a result of the truncated sumin Eq. 27. Hence, above we justify that the remainder of P projections after truncation canbe omitted up to a certain precision. On the other hand, in order to describe cumulativecontribution of singular vectors we look at the L distance of the omitted remainder (i.e. k ≥ m , l ≥ n ) E ( m, n ) = (cid:88) i,j | ∆ P ij ( m, n ) | = (cid:88) i,j (cid:12)(cid:12)(cid:12) K,L (cid:88) k = m,l = n X ik S kl Y jl (cid:12)(cid:12)(cid:12) (59)which describes how well the method converges self-consistently. Here we didn’t normalizethis equation for every term we added. Which means we just rely on fast convergence of P projections measured by Eq. 58, so first few error terms can be misleading. III. RESULTS
In order to treat strong correlations we devised a robust numerical method to study thejoint statistics of membrane potentials and spike trains of integrate-and-fire model neurons.The case study reported here covers the leaky integrate-and-fire (LIF) model with Pois-son input spike trains. However, our method can be easily generalized to non-linear leakfunctions [35], conductance based synaptic inputs [36] and more complex input correlationmodels [30, 31],. although we have to leave the details of such generalizations open. In thissection we will explain how to select a ‘good basis for expansion’, and we will give numericalexamples that demonstrate the power of the method.18 . SVD of joint probability distributions and choice of expansion basis
We started from a simple observation: The stationary joint membrane potential distri-bution for two neurons with independent input is given as P ( V , V ) = P ( V ) P ( V ) (60)where P ( V ) and P ( V ) are the stationary membrane potential distributions of two inde-pendent neurons, as described in Eq. 4. A similar relation for a discretized voltage grid canbe written as P = P ⊗ P . (61)For the case of shared input this simple relation is not valid any more. On the other hand,we observed that a value of the parameter c close to 0 will practically recruit only a smallnumber of additional principal components for any given precision, cf. Fig. 4. Here weperform a singular value decomposition of the full solution of Eq. 32 given in terms of thematrix P ij P = L Σ R T or P ij = N (cid:88) k =1 σ k L ik R kj (62)generalizing the case of independent neurons to also reconstruct the joint membrane potentialdistribution for neurons with shared input. As demonstrated in Fig. 4, convergence is ratherquick, even for moderate values of c .Another aim of our study was to gain some understanding about the influence of thespace-time grid. We observed that left and right singular vectors are of the form P kl = (cid:88) ij Ω ij ( X cik + a i δ k,r )( Y cjl + a j δ l,r ) (63)where X c and Y c reflect the quasi-continuous part of basis vectors and Ω is the couplingmatrix as defined in Eq. 27. Here we need to make sure that the emerging singularity at thereset bin is not causing any numerical problems. One needs to first consider a small timestep ∆ t and adapt the stepping in space ∆ x accordingly. A more thorough discussion of asuitable coarse graining strategy, however, is postponed to a later section of this paper.Here we suggest to use SVD as a method to achieve a dimensional reduction of the fullsystem. As it is a numerical method, its convergence and efficiency needs to be addressed.Generally, there are several different options to select a basis. Specifically, we use the19ight singular vectors of single-neuron Markov matrices. As demonstrated in Fig. 4, rightsingular vectors lead to an expansion that converges faster for coarse grids (e.g. V = 0 . V = 0 .
05 mV) the difference is less prominent (Figs. 5a and b),right singular vectors still converge slightly faster than left singular vectors (Fig. 5d). Rightsingular vectors of the single-neuron time evolution matrix yield an orthogonal coordinatesystem with very good properties. ac bd
FIG. 5. Comparison of using right or left singular vectors for a reconstruction of the joint mem-brane potential distribution. We observe that the right singular vectors have better convergenceproperties. (a) Mode coupling matrix Ω for a basis derived from right singular vectors. (c) Partialsum error (Eq. 59) for the basis of right singular vectors, corresponding to (a). (b) Mode couplingmatrix Ω for a basis derived from left singular vectors. (d) Difference of partial sum errors for leftsingular vectors corresponding to (b) and right singular vectors corresponding to (a). Red colorindicates positive sign, while blue color indicates negative sign of the error. The reconstructionwith right singular vectors converges slightly faster. Note that for the error measures consideredin (c) and (d) we didn’t take into account the bottom left 10 ×
10 entries of the matrix.
As reported previously [20], we may also use a direct analytical approach using the basis20
IG. 6. Convergence of the SVD-based approximation method using up to 50 singular vectorscorresponding to the largest singular values. (a) Mode coupling matrix Ω, defined by P = X T Ω Y (color represents log (Ω ij / max ij (Ω ij ))). (b) Reconstructed 2D membrane potential distributionbased on a coarse graining with 50 ×
50 grid points. (c) log of relative L error. The value givenat location ( i, j ) is the contribution to the reconstruction of P computed via summation of allvectors n > i , m > j (Eq. 59). (d) Error that arises from (cid:80) k X ik (cid:54) = 0 (Eq. 58). of the single-neuron Fokker-Planck operator, and its adjoint basis P ( X, Y ) = (cid:88) ij Ω ij X i ⊗ Y j (64)where X and Y are the left eigenvectors of the single neuron operators. However, the issueis that the discrete adjoint basis blows up at the lower boundary. The effect of this on ourapproximation is demonstrated in Fig. 3. In general, SVD eliminates a kernel of singularmatrices.In our treatment of the 2D Fokker-Planck equation, which is the infinitesimal limit of thetheory considered above, we used the basis and adjoint basis to project linear operators toa subspace. This has certain advantages as it satisfies constraints for marginal distributionsand preserves the Markov property to some extent. Positivity of the solution in the subspace21s not guaranteed, but time evolution is probability preserving ( (cid:80) i P i ( t ) = 1).First, SVD is computationally convenient, because it leads to using a real orthogonalbasis which resembles the eigenbasis of M . Second, numerical instabilities due to an ill-conditioned time evolution matrix M (some eigenvalues λ i ≈
0) are cured by SVD. Third,although the basis vectors implied by SVD have the disadvantage of not completely preserv-ing positivity, the deviation remains within tight bounds even for a relatively small numberof basis vectors.
B. Comparison to direct numerical simulations
We compare our SVD-based Markov theory and direct numerical simulations of spikingneurons both on the level of joint 2D membrane potential distributions and on the levelof spike train covariance functions, cf. Fig. 7. The empirical distributions derived from 2Dhistograms are slightly smoothed in order to compare them to the distributions derivedfrom the Markov theory on the level of contour lines. We also considered 2D cumulativedistribution functions, where the smoothing step can be omitted. Moreover, we computedoutput spike train covariance functions as described in methods section and compared themto the covariance functions obtained directly from the simulated spike trains.
C. Application 1: Non-linear correlation transfer
Two neurons that are driven by correlated input current will exhibit correlated outputspike trains. This transfer of correlation reflects an important property of neuronal dynamics,which is of particular relevance for understanding the contribution of neurons to networkdynamics. Recently, we were able to demonstrate, by exact analytical treatment, that thecorrelation transfer for leaky integrate-and-fire neurons is strongly non-linear [20]. Onlyfor weak input correlation it can be described by perturbative methods, and deviationsfrom linear response theory depend on the background firing rate. In the present work wedemonstrate the same non-linear correlation transfer, cf. Fig. 10. There we also demonstratehow the parameters of the spatial and temporal coarse graining affects the precision ofthe Markov approximation. Our main conclusion is that dimensional reduction via SVDsubspace projections makes it possible to achieve a superior precision with small bin sizes.22 bc d
FIG. 7. Effect of different parameters of neurons or input to neurons (here, σ asymmetry) on thejoint membrane potential distribution and spike cross-covariance function. Reconstruction of the2D joint membrane potential distribution using SVD (∆ V = 0 . t = 0 . Fine enough grids, however, could not be dealt with on typical computers without using thereduction suggested here.
D. Application 2: Asymmetric cross-covariance functions
Neurons in biological networks have widely distributed parameters, and this heterogeneitymay also influence information processing [38–40]. Moreover, robust asymmetries in spikecorrelations could lead to asymmetric synaptic efficacies, if they are subject to spike timingdependent plasticity [41, 42]. Our approach reveals a generic temporal asymmetry in cross-23
Time (ms) C o v a r i a n ce ( / s )
20 0 20
Time (ms) C o v a r i a n ce ( / s )
20 0 20
Time (ms) C o v a r i a n ce ( / s )
20 0 20
Time (ms) C o v a r i a n ce ( / s ) a bc d FIG. 8. Asymmetric cross-covariance functions in the strongly correlated regime ( c = 0 . µ asymmetry, µ (cid:54) = µ while all other paramters are the same, (b) σ asymmetry, σ (cid:54) = σ , (c) τ m asymmetry, τ (cid:54) = τ , which leads to µ (cid:54) = µ and σ (cid:54) = σ as privatespike train input rates are the same, (d) V th asymmetry, V th, (cid:54) = V th, . For specific parametervalues, see Table II. covariance functions, related to the heterogeneity of intrinsic neuron parameters and inputvariables. Such temporal asymmetry is more pronounced for larger values of c , especially inthe non-perturbative regime that we address in this work.We document here an application of our method to four types of asymmetries [20, 40]: µ asymmetry, σ asymmetry, τ m asymmetry and V th asymmetry. We quantified the asymmetryby a = χ /χ (specific values given in Table II), where χ is replaced by the respectiveparameter. Asymmetric correlations have been described previously, and they were bynumerical simulations and experimentally studied by [38, 40, 43, 44].*** 24 Voltage V o l t ag e P max P max Voltage V o l t ag e P max P max Voltage V o l t ag e P max P max Voltage V o l t ag e P max P max a bc d FIG. 9. Joint 2D membrane potential distribution of simulated neuron dynamics for c = 0 . w = 1 mV) is compared to the joint distributioncomputed with the SVD method (using a subspace of dimension 50 × s . Both types of non-equal neuronparameters can lead to similar distributions. Our method can deal with all such cases accurately.Results are presented here for different asymmetric parameters: (a) µ asymmetry, (b) σ asymmetry,(c) τ m asymmetry, which implies an asymmetry in µ and σ as well, as private spike train inputrates are the same, (d) V th asymmetry. For specific parameters, see Table II. IV. DISCUSSIONA. Relevance of the new approach presented here
Models of correlated neuronal activity describe the origin of correlations in spiking modelneurons, induced by the structure of the network and/or feedforward input. Such neuronmodels, however, are notoriously nonlinear. Nevertheless, most treatments rely on lineariza-tion and other simplifying assumptions, as nonlinear correlation transfer functions (i.e. rela-25 a bd e
Finer Space Grid
FIG. 10. Limits to the precision of cross-covariance functions and correlation transfer functions.(a) Correlation transfer function as a function of input correlation. We compare here analyticalresults (solid curves) described in a recent paper [20] with numerical results (orange symbols)obtained with the methods described in this paper. Non-perturbative correlation transfer functions C out ( C in ) in [20] for symmetric parameters and for high and low firing rates, respectively (blue: r b = 15 . = 0 .
5; green: r g = 1 .
13 Hz, CV = 0 . dC out dC in at C in = 0) are computed using perturbation theory [37]. Note thatwe added the obvious points C out (0) = 0 and C out (1) = 1 to the plot by hand. (b) Cross-covariancefunctions C ( τ ) (solid red curves, with unit Hz ) as a function of the lag τ . For non-infinitesimalPSPs there is a delta function at zero lag τ = 0 (blue stems, with unit Hz ), the amplitude of whichgrows as c increases. Figures from top left to bottom right correspond to different values of c . For(a) and (b) we chose ∆ V = 0 .
05 mV. Panels (c) and (d) are zoomed-in versions of the c = 0 . c = 0 .
95 (bottom) covariance functions (solid green curves) to demonstrate the effect ofthe grid (c) ∆ V = 0 .
05 mV vs. (d) ∆ V = 0 .
02 mV vs. (e) ∆ V = 0 .
01 mV. Further parameters aregiven in Table II . IG. 11. Here we demonstrate that fixing the space bin (here ∆ V = 0 .
02 mV), the choice of thetime bin affects the firing rate estimation [24]. The deviation of correlation coefficients for smallrates in Fig. 10a (orange triangles vs. dark green curve) is a result of a poor estimation of conditionalrates. The variance of the input ( σ ) is crucial in determining the appropriate temporal bin size,while the mean input ( µ ) is less effective. With increasing variance one observes an increasingfiring rate error. From these plots, we conclude that for a space bin ∆ V = 0 .
02 mV, a time binbetween ∆ t = 0 . t = 0 . tion of input and output correlations) are difficult to derive. Previous analytical approacheshave employed perturbation theory [16, 17] to study pairwise correlations under the assump-tion of weak input correlation [37, 45]. As a consequence of this technical convenience, we27till lack a systematic approach that allows us to deal with a broad range of correlations,and to gain an understanding of their implications for network dynamics. B. Extensions of the theory
All computations described in our paper can be applied to more general integrate-and-firemodels with anon-linear leak function Ψ( V ) τ m ˙ V = Ψ( V ) + J ex τ m (cid:88) k S exk ( t ) + J in τ m (cid:88) k S ink ( t ) . (65)We only need to rewrite the decay matrix as a discrete approximation of the differentialoperator D ( x ) = τ m ddt x − Ψ( x ).Other scenarios of interest are reflected by an altered amplitude distribution of the inputs.This is a natural consequence if individual synapses have different PSP amplitudes. It couldalso arise, however, if the population of input neurons has a non-trivial correlation structure.In particular, higher-order correlations have been treated in terms of specific amplitudedistributions [31, 46]. The method described in the present paper can be adapted to suchscenarios by simply using a modified definition of the jump matrix [30].Higher-order statistics on the output side is also compatible with our method, describingthe joint response behavior of three or more neurons that are driven by shared input. Third-order correlations can be computed in practice, because the projections work in the samefashion P lmn = (cid:88) ijk Ω ijk X il Y jm Z kn (66)now with a 3D jump operator given in the generic form J D = J ⊗ J ⊗ J . (67)This operator is again transformed with a basis derived from a SVD as J ⊗ J ⊗ J → XJ x X T ⊗ Y J Y T ⊗ ZJ Z T . (68)This procedure is computationally more demanding as we need to consider additional paths,although the scaling is not exponential. Under assumption of homogeneous shared input28same jump amplitudes driven by shared input in all directions) leads to an expressionsimilar to Eq. 37, J D = (cid:88) j ∈ J ψ ( j )[ e J O j ⊗ e J O j ⊗ e J O j ] . (69)Assuming joint stationarity of all three spike trains S ( t ), S ( t ) and S ( t ), we need to findthe joint third moment of the spike train statistics µ ( τ , τ ) = (cid:10) S ( τ ) S ( τ ) S (0) (cid:11) . (70)As shown above, second moments can be computed with our method (Fig. 9). In order toobtain the covariance function from the stationary 3D flux, the time evolution of the 2Dconditional flux at times ( τ , τ ) is needed. This is given as∆ t P | ( t ) = M P | ( t ) (71)which is computationally demanding as the numerical effort scales as O ( N ). However, thiscan be projected to the subspace with time dependent coupling matrix Ω as∆ t Ω( t ) = Q Ω( t ) . (72)This form has advantages over finite difference methods as e.g. suggested in [23]. Thecomputation of the third order moment defined in Eq. 70 requires a solution of Eq. 72at τ to find the second conditional distribution P | | ( τ ). Then we need to find the 1Dconditional distribution (e.g. for neuron 1) at t = τ by solving∆ t P | | ( t ) = M P | | ( t ) (73)similar to Eq. 23. This provides us with conditional rate r | | ( τ , τ ) and then the thirdmoment is given as µ ( τ , τ ) = r | | ( τ , τ ) r . (74)Details of this computation have to be deferred to future work, though. C. Boundary conditions and singularity
The joint membrane potential distribution has a singularity at the origin ( V , V ) =(0 ,
0) due to a coordinated reset caused by some shared input spikes. There is also a line29iscontinuity at V = 0 and V = 0, again due to the reset boundary condition. Thesesingularities are reflected in the right singular vectors X and Y . This is the exact reasonwhy we selected them as a basis to expand operators and joint distributions.We observed that a singularity (a δ -function) emerges when ∆ V is small and ∆ t is large,in relative terms. This is an issue even for the 1D discrete problem, and it is even moresevere for 2D problems as the amplitude of the singularity scales quadratically with ∆ V .This phenomenon occurs only if PSPs have a finite amplitude. As the PSP gets largerrelative to ∆ V , reset currents remain finite even in continuous time [24]. As a consequence,the limit to continuous variables must be taken with care, in particular for c > δ -singularity does not exist for the diffusion approximation [20]. However, the def-inition of the current at the origin again fails as the derivative is discontinuous in both V and V directions. The infinitesimal limit of the jump equation must be taken with care.There is no doubt that the jump equation is well-defined as the flux at the boundary is notlocal. However, the infinitesimal limit is problematic for correlated neurons ( c >
0) as theflux is not defined at the boundary of the 2D domain.
D. Precision, computational efficiency and grid selection
The selection of an appropriate grid in space and time is crucial for correlation compu-tations. The small residual offset between direct simulation and our new semi-analyticalcomputation (cf. Fig. 8 and 9), for example, can be considered as a discretization arti-fact. Although this issue would deserve a more systematic treatment, we report here someobservations that can guide grid selection:(i) For discrete solutions of the heat equation based on central difference scheme, con-vergence of 1D time evolution requires ∆ t (∆ V ) < τσ [47]. A similar rule also appliesin the 2D case considered here. In general, explicit discretization schemes of secondorder differential operators arising in the study of diffusion, require positivity and sta-bility conditions in the order of ∆ t = O ( τ (∆ V +∆ V ) σ ) [23]. In this work, we followed adiscretization scheme that approximately conserves probability, a Markovian approxi-mation [24]. However, we note once more that some grids may lead to violation of theMarkov property for too large ∆ t , as a result of boundary effects. This may createissues when the largest eigenvalue exceeds 1.30ii) To reflect small expected bin counts (especially for c ≈
1) adequately, one needs alarger ∆ t and a smaller ∆ V . This is in conflict with rule (i). Besides, we observethat smaller ∆ V for a fixed ∆ t actually leads to better firing rate approximation upto some point (Fig. 11).(iii) A finer grid requires more computational efforts to achieve a smooth correlation func-tion. The SVD reduction does not alter this behavior. Other dependencies and limitingfactors are indicated in Fig. 10. There are two constraining factors which are deter-mined by the selected precision of the approximation. One is the extent of the jumpdistribution, which affects the number of terms to be accumulated (size of the set J inEq. 37). For a fixed grid and a selected precision, this number increases with σ c . Thesecond constraint is the size of the SVD subspace. We know that as c gets closer to 1and C ( τ ) gets steeper we need to include more singular components.In Fig. 10 we illustrate how coarse graining affects the shape of the cross-covariancefunction C ( τ ). Although the precision of the approximation is limited by the subspaceprojections implied by SVD, the grid parameters ∆ t and ∆ V are the most important factorsto get the shape of the function right. However, for a fixed dimension of the SVD subspaceeven the finest grid would not be able to capture the singularity at zero time lag ( τ = 0).The grid effectively limits the precision of the approximation due to the reduced number ofdegrees of freedom. V. CONCLUSION
We developed a novel numerical method to compute the joint statistics and correlationfunctions for two LIF model neurons driven by shared input. Our approach can deal withthe full range of input correlations c , and the expansion converges fast. Also, our methodis widely generalizable and can deal with other scenarios that are biologically relevant.We observed in previous work [20] that low output firing rates generally require a non-perturbative treatment. If output rates are high, in contrast, and for high input firing rateswith small PSPs (diffusion regime), the approximation derived from linear response theory[37] is reasonably precise.We conclude that it is possible to compute correlation functions (in contrast to deriving31hem from simulations) for a wide range of models with finite PSP amplitudes, and also fora wide range of parallel spike train input models. Although there is currently no conclusivetheory for the selection of an appropriate spatio-temporal grid, we were able to come up withsome heuristics. The precision of even the first moment (firing rate) depends on the grid.Specifically if c is close to 1, the temporal component of the correlation function resemblesa delta function. In order to capture this phenomenon the grid must be fine enough.The innovation in our work is not only the formulation of correlation functions based on aMarkov chain approximation, but also a dimensional reduction. This helps us compute jointmembrane potential distributions. We showed that the number of components obtained bySVD needed to represent single neuron dynamical evolution matrices is small. This alsomeans that computations can presumably be generalized to higher-order correlations withonly moderately increased computational effort.Systematic benchmarking of our method has not yet been performed. However, we believeour method constitutes the only reasonable numerical approximation to the joint statisticsof strongly correlated neuronal dynamics with finite PSPs, apart from direct stochasticsimulations [32]. This approximation for reasonable grids in space and time was only viablewith SVD subspace projections. 32 PPENDIX: PARAMETERSTable 1.
Parameters for NEST simulations and semi-analytical computations
Neuron parameters: ( Fig. 3, Fig. 5, Fig. 6, Fig. 7)
Symbol Description Value V th voltage threshold 15 mV V r voltage reset 0 mV τ m membrane time constant 15 ms τ ref refractory period 1 ms h PSP 0.01 mV -0.1 mV∆ t time resolution 0.1 ms Input parameters µ mean input 10-15 mV σ , σ STD private input 2-10 mV σ c STD shared input 2-10 mV c input noise correlations 0-1 a asymmetry factor > able 2. Numerical results vs. NEST simulations : (Fig. 8, Fig. 9, Fig. 10)
Correlation asymmetry parameters
Symbol Description Value V th voltage threshold 15 mV V r voltage reset 0 mV τ m membrane time constant 15 ms τ ref refractory period 1 ms h PSP 0.1 mV µ mean input 12. mV σ STD total input 5. mV for (b) and (d) , 6. mV for (a) and (c) c input noise correlations 0.9 a σ asymmetry factor 1 / √ a µ asymmetry factor 10/13 a τ asymmetry factor 10/15 a V th asymmetry factor 13/15*** asymmetry factors: a = χ /χ . able 3 Neuron model parameters Symbol Description Unit S ex , S in spike trains 1/ms V membrane potential mV V th voltage threshold mV V r voltage reset mV t time ms τ m, , τ m, membrane time constant ms τ ref refractory period ms h ex = h , h in = gh PSP mV µ mean input mV σ , σ STD private input mV σ c STD shared input mV σ STD total input mV c input noise correlations 0-1 r e , r i ex & inh input rates Hz r , r output rates of 2 neruons Hz Table 4 Correlations and related notation
Symbol Description Unit C ( τ ) covariance function Hz r | ( τ ) conditional rate Hz P ( V ) probability distribution of V 1/mV P | ( V ) conditional probability distribution of V / mV P ( V , V ) joint probability distribution of ( V , V ) 1 / mV ∆ t discrete time evolution operator 1 M discrete time evolution matrix 1 able 5 Probability distributions and Markov approximation Symbol Description Unit P ( q, r ) probability for q excitatory, r inhibitory spike 1 P ( γ ) probability for jumps of length γ P ( V , V | V , , V , ) probability forjumps from ( V , , V , ) to ( V , V ) 1 P ( V , V ) probability for jumps from (0 ,
0) to ( V , V ) 1 T threshold matrix : N × ( N + n ) 1 J jump matrix : ( N + n ) × N D decay matrix: N × N M time evolution matrix : N × N T D threshold tensor : N × M × ( N + n ) × ( M + m ) 1 J D jump tensor : ( N + n ) × ( M + m ) × N × M D D decay tensor: N × M × N × M M D Time evolution tensor 2D : N × M × N × M U , D up or down operators in discrete space 1 J p , J p private V or V jump generators 1 J c , J c shared V or V jump generators 1 c mn coefficient of m excitatory and n inhibitory jumps 1 P stationary probability density in discrete space (mV) − P J,k probability for jumps k -th component of J P / mV P flux,k probability for jumps k -th component conditional exit flux 1 / mV∆ t time bin ms∆ V voltage bin mV a , b average count of private input 1 in ∆ t a , b average count of private input 2 in ∆ t a c , b c average count of shared input in ∆ t able 6 SVD reduction Symbol Description Unit L , R SVD left and right basis 1Σ singular value matrix 1 X , Y SVD subspaces 1˜ X , ˜ Y extended subspace 1 Q reduced operators defined on a selected subspace 1 M or N size of full grid i.e. matrix 1 m or n maximum number of jumps over the threshold 1 K or L size of SVD subspace 1 M + m or N + n size of full jump subspace 1 K + m or L + n size of reduced jump subspace 1 Table 7 Asymmetry parameters
Symbol Description Unit a σ asymmetry factor σ /σ a µ asymmetry factor µ /µ a τ asymmetry factor τ /τ a V th asymmetry factor V th, /V th, Frontiers in ComputationalNeuroscience , Vol. 8 (2014).[2] I. Lampl, I. Reichova, and D. Ferster, Neuron , 361 (1999).[3] M. Okun and I. Lampl, Nature Neuroscience , 535 (2008).[4] J. F. Poulet and C. C. Petersen, Nature , 881 (2008).[5] A. Arieli, A. Sterkin, A. Grinvald, and A. Aertsen, Science (New York, N.Y.) , 1868(1996), arXiv:arXiv:1011.1669v3.[6] B. Staude, S. Rotter, and S. Gr¨un, Neural Computation , 1973 (2008).[7] A. Kumar, S. Rotter, and A. Aertsen, Nat Rev Neurosci , 615 (2010).[8] B. Doiron, A. Litwin-Kumar, R. Rosenbaum, G. K. Ocker, and K. Josi´c, Nature Neuroscience , 383 (2016).[9] A. S. Ecker, P. Berens, G. a. Keliris, M. Bethge, N. K. Logothetis, and A. S. Tolias, Science(New York, N.Y.) , 584 (2010).[10] A. Renart, J. D. Rocha, P. Bartho, L. Hollender, N. Parga, A. Reyes, and K. D. Harris, ,587 (2010).[11] V. Pernice, B. Staude, S. Cardanobile, and S. Rotter, PLoS Computational Biology (2011),10.1371/journal.pcbi.1002059.[12] M. Helias, T. Tetzlaff, and M. Diesmann, PLoS Computational Biology (2014),10.1371/journal.pcbi.1003428, arXiv:1304.2149.[13] N. Brunel, Neurocomputing , 307 (2000).[14] V. Pernice, B. Staude, S. Cardanobile, and S. Rotter, Physical Review E - Statistical, Non-linear, and Soft Matter Physics , 1 (2012), arXiv:arXiv:1201.0288v2.[15] J. Trousdale, Y. Hu, E. Shea-Brown, and K. Josi´c, PLoS Computational Biology (2012),10.1371/journal.pcbi.1002408, arXiv:1110.4914.[16] N. Brunel and V. Hakim, Neural Computation , 1621 (1999), arXiv:9904278 [cond-mat].[17] B. Lindner and L. Schimansky-Geier, Physical Review Letters , 2934 (2001).[18] S. Ostojic, Nature neuroscience , 594 (2014).[19] S. Jovanovi´c, J. Hertz, and S. Rotter, Phys. Rev. E , 042802 (2015).[20] T. Deniz and S. Rotter, , 1 (2016), arXiv:1604.03619.[21] M. N. Shadlen and W. T. Newsome, The Journal of neuroscience : the official journal of theSociety for Neuroscience , 3870 (1998).[22] D. R. Lyamzin, S. J. Barnes, R. Donato, J. A. Garcia-Lazaro, T. Keck, and N. A. Lesica,Journal of Neuroscience , 8065 (2015).[23] R. Rosenbaum, F. Marpeau, J. Ma, A. Barua, and K. Josi´c, Journal of Mathematical Biology , 1 (2012), arXiv:1011.0669.[24] M. Helias, M. Deger, M. Diesmann, and S. Rotter, Frontiers in computational neuroscience , 29 (2010).[25] A. Abouzeid and B. Ermentrout, Physical Review E - Statistical, Nonlinear, and Soft MatterPhysics (2011), 10.1103/PhysRevE.84.061914, arXiv:arXiv:1101.1919v1.[26] T. Tchumatchenko, A. Malyshev, T. Geisel, M. Volgushev, and F. Wolf, Physical ReviewLetters , 2 (2010), arXiv:0810.2901.
27] R. D. Vilela and B. Lindner, Physical Review E - Statistical, Nonlinear, and Soft MatterPhysics , 1 (2009), arXiv:0912.2336.[28] T. Schwalger, F. Droste, and B. Lindner, Journal of Computational Neuroscience , 29(2015).[29] S. O. Voronenko, Stannat W, and Linder B, The Journal of Mathematical Neuroscience , 1(2015).[30] M. Schultze-Kraft, M. Diesmann, S. Gr¨un, and M. Helias, PLoS Computational Biology ,e1002904 (2013), arXiv:1207.7228.[31] A. Kuhn, A. Aertsen, and S. Rotter, Neural Computation , 67 (2003).[32] M. J. E. Richardson and R. Swarbrick, Phys. Rev. Lett. , 178102 (2010).[33] M. Abramowitz, Handbook of Mathematical Functions, With Formulas, Graphs, and Mathe-matical Tables, (Dover Publications, Incorporated, 1974).[34] M.-O. Gewaltig and M. Diesmann, Scholarpedia , 1430 (2007).[35] W. Gerstner and W. M. Kistler, Spiking neuron models : single neurons, populations, plasticity (Cambridge University Press, Cambridge, New York, Melbourne, 2002) autre tirage : 2006,2008 (4e).[36] A. Kuhn, A. Aertsen, and S. Rotter, Integration The Vlsi Journal , 2345 (2004).[37] J. de la Rocha, B. Doiron, E. Shea-Brown, K. Josi´c, and A. Reyes, Nature , 802 (2007).[38] K. Padmanabhan and N. N. Urban, Nature Neuroscience , 1276 (2010).[39] M. Y. Yim, J. Wolfart, A. Aertsen, and S. Rotter, Frontiers in Computational Neuroscience (2012), 10.3389/conf.fncom.2012.55.00030.[40] M. Y. Yim, A. Aertsen, and S. Rotter, Physical Review E - Statistical, Nonlinear, and SoftMatter Physics , 1 (2013), arXiv:1208.5350.[41] A. Morrison, M. Diesmann, and W. Gerstner, Biological Cybernetics , 459 (2008).[42] B. Babadi and L. F. Abbott, PLoS Computational Biology , e1002906 (2013).[43] S. Ostojic, N. Brunel, and V. Hakim, J Neurosci , 10234 (2009).[44] M. Y. Yim, A. Kumar, A. Aertsen, and S. Rotter, Journal of Computational Neuroscience ,293 (2014).[45] E. Shea-Brown, K. Josi´c, J. de la Rocha, and B. Doiron, Physical Review Letters , 108102(2008).[46] B. Staude, S. Rotter, and S. Gr¨un, in Analysis of Parallel Spike Trains , edited by S. Rotter nd S. Gr¨un (Springer, 2010) Chap. 12, pp. 253–283.[47] W. F. Ames, Numerical methods for partial differential equations (Academic press, 2014).(Academic press, 2014).