Joint statistics of strongly correlated neurons via dimensionality reduction

The leaky integrate-and-fire neuron model with postsynaptic potentials of finite amplitudes was studied previously in [25]. The stochastic equation that describes the membrane potential dynamics of one particular neuron is given as

$$\begin{align} \tau_m\dot{V}(t)=-V(t)+h_{{\rm ex}}\tau_m S^{{\rm ex}}(t) +h_{{\rm in}} \tau_m S^{{\rm in}}(t) \nonumber \end{align} \tag{ 1 }$$

where $h_{{\rm ex}}$ and $h_{{\rm in}}$ represent the amplitudes of individual EPSPs and IPSPs, respectively, and $S^{{\rm ex}}(t)$ and $S^{{\rm in}}(t)$ are the spike trains of excitatory and inhibitory presynaptic neurons, with each of their spikes represented by a Dirac delta function

$$\begin{align} S^{{\rm ex}}(t)=\sum_i \sum_{n}\delta(t-t^{i}_n). \nonumber \end{align} \tag{ 2 }$$

An analogous definition holds for inhibitory neurons. In both cases, if the membrane potential reaches the firing threshold, $V_{{\rm th}}$, a spike is elicited and the voltage is reset to its resting value at 0 (figure 1(c)).

Zoom In Zoom Out Reset image size

In order to study correlations between the spike trains of two neurons we look at two coupled stochastic equations, describing the membrane potentials of two neurons that share a certain fraction of their excitatory and inhibitory input spikes

$$\begin{align} \dot{V_1}(t)=-\frac{V_1(t)}{\tau_{m,1}}+h_{{\rm ex},1}(S^{{\rm ex}}_{1}(t)+S_{c}^{{\rm ex}}(t))+ h_{{\rm in},1}(S^{{\rm in}}_{1}(t) + S_{c}^{{\rm in}}(t)) \nonumber \end{align} \tag{ 3 }$$

$$\begin{align} \dot{V}_2(t)=-\frac{V_2(t)}{\tau_{m,2}}+h_{{\rm ex}.2}(S^{{\rm ex}}_{2}(t)+S_{c}^{{\rm ex}}(t))+ h_{{\rm in},2}(S^{{\rm in}}_{2}(t) + S_{c}^{{\rm in}}(t)) \nonumber \end{align} \tag{ 4 }$$

where shared excitatory and inhibitory input spike trains are denoted as $S^{{\rm ex}}_{c}(t)$ and $S^{{\rm in}}_{c}(t)$, respectively. Comparing the parameters of the jump model to shared white noise input we implicitly specified the firing rates of 6 independent Poisson processes (figure 1)

$$\begin{align} \dot{V_1}(t)=-\frac{V_1(t)}{\tau_{m,1}}+\mu_1+\frac{\sigma_1}{\sqrt{\tau_{m,1}}} [\sqrt{1-c}\:\xi_1\pm\sqrt{c}\:\xi_c]) \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \dot{V}_2(t)=-\frac{V_2(t)}{\tau_{m,2}}+\mu_2+\frac{\sigma_2}{\sqrt{\tau_{m,2}}} [\sqrt{1-c}\:\xi_2\pm\sqrt{c}\:\xi_c]). \label{eq_dif_model} \nonumber \end{align} \tag{ 6 }$$

We present here one way to parametrize the white noise input with a shared component: $\frac{\sigma_a}{\sqrt{\tau_{m, a}}} [\sqrt{1-c}\:\xi_a\pm\sqrt{c}\:\xi_c])$ for a  =  {1, 2}, where the parameter c is the input correlation coefficient. The diffusion limit above was studied in [21] employing a spectral expansion of 2D Fokker–Planck equation. Later we compare the results achieved with two different models (Poisson versus diffusion) in figure 11. We parametrize Poisson input firing rates (as in equations (16) and (17)) by solving the corresponding white noise input parameter equations (as in equations (9) and (10)). Additionally, for notational convenience, shared input rates $r_{{\rm ex}, c}$ and $r_{{\rm in}, c}$ are computed from a shared Wiener process $\xi_c$ with zero mean. Setting $h_{{\rm ex}, a}=h_a$ and $h_{{\rm in}, a}=gh_a$, Campbell's theorem for shot noise [32] states that mean and variance of the corresponding diffusion approximation are given by

$$\begin{align} \mu_{a}=\tau_{m,a}h_a(r_{{\rm ex},a}-gr_{{\rm in},a}) \nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \sigma^2_{a} (1-c)=\tau_{m,a}h_a^2(r_{{\rm ex},a}+g^2 r_{{\rm in},a}) \nonumber \end{align} \tag{ 8 }$$

where $a ={1, 2}$ is the neuron index. (This can be applied to the discrete time case we deal with here.) Additionally, the shared spike train rates are found from

$$\begin{align} \label{mu_r} 0=\tau_{m,1}h_1(r_{{\rm ex},c}-gr_{{\rm in},c}) \nonumber \end{align} \tag{ 9 }$$

$$\begin{align} \sigma^2_{1} c=\tau_{m,1}h_1^2(r_{{\rm ex},c}+g^2 r_{{\rm in},c}) \label{sigma_r} \nonumber \end{align} \tag{ 10 }$$

as well as

$$\begin{align} 0=\tau_{m,2}h_2(r_{{\rm ex},c}-gr_{{\rm in},c}) \nonumber \end{align} \tag{ 11 }$$

$$\begin{align} \sigma^2_{2} c=\tau_{m,2}h_2^2(r_{{\rm ex},c}+g^2 r_{{\rm in},c}) \nonumber \end{align} \tag{ 12 }$$

which demands that parameters satisfy

$$\begin{align} \label{the_constraint_phi} \frac{\sigma^2_{1}}{\tau_{m,1}h_1^2}=\frac{\sigma^2_{2}}{\tau_{m,2}h_2^2}=\phi. \nonumber \end{align} \tag{ 13 }$$

By solving equations (9) and (10), the firing rates for each independent excitatory and inhibitory process are given as

$$\begin{align} r_{{\rm ex},a}=\frac{1}{1+g}\Big{(}\frac{\sigma_a^2(1-c)}{\tau_{m,a} h_a^2}+\frac{g\mu_a}{h_a\tau_{m,a}}\Big{)} \nonumber \end{align} \tag{ 14 }$$

$$\begin{align} r_{{\rm in},a}=\frac{1}{g(1+g)}\Big{(}\frac{\sigma_a^2(1-c)}{\tau_{m,a} h_a^2}-\frac{\mu_a}{h_a\tau_{m,a}}\Big{)} \nonumber \end{align} \tag{ 15 }$$

considering equation (13) and defining $\varphi_a=\frac{\mu_a}{h_a\tau_{m, a}}$, this becomes

$$\begin{align} \label{r_ex} r_{{\rm ex},a}=\frac{1}{1+g}\Big{(}(1-c)\phi+g\varphi_a\Big{)}\nonumber \end{align} \tag{ 16 }$$

$$\begin{align} r_{{\rm in},a}=\frac{1}{g(1+g)}\Big{(}(1-c)\phi-\varphi_a\Big{)} \label{r_in} \nonumber \end{align} \tag{ 17 }$$

and similarly

$$\begin{align} \label{r_c_ex} r_{{\rm ex},c}=\frac{c\phi}{1+g}\nonumber \end{align} \tag{ 18 }$$

$$\begin{align} r_{{\rm in},c}=\frac{c\phi}{g(1+g)}. \label{r_c_in} \nonumber \end{align} \tag{ 19 }$$

Here we note that some combinations of $(\mu, \sigma)$ do not correspond to a combination of positive Poisson rates, as they must satisfy

$$\begin{align} \sigma^2\geqslant h\mu \nonumber \end{align} \tag{ 20 }$$

to guarantee $r_{{\rm in}}\geqslant 0$ for the inhibitory rate.

In this section we summarize the discrete approximation to the dynamics of a LIF neuron, as developed in [25]. The coarse graining of the membrane potential is given by the map

$$\begin{align} \mathbb{R} \rightarrow \mathbb{Z}, \quad V \mapsto \lfloor \frac{V}{\Delta V}\rfloor. \nonumber \end{align} \tag{ 21 }$$

The probability density function P(v) of the membrane potential then becomes a vector $\boldsymbol{p}$ with components satisfying

$$\begin{align} \label{p_little} p_i = \frac{1}{\Delta V}\int^{v_{i+1}}_{v_i} P(V)\,{\rm d}V, \nonumber \end{align} \tag{ 22 }$$

with the voltage vector component

$$\begin{align} \label{v_little} v_i=i \, \Delta V. \nonumber \end{align} \tag{ 23 }$$

As we impose a cut-off lower boundary V₋ of the voltage scale, the dimension of the discrete state space is given as $N=\frac{V_{{\rm th}}-V_-}{\Delta V}$.

The temporal evolution of the membrane potential distribution is now described in terms of a Markov process. The Markov propagator can be expressed as a juxtaposition of three operators: a decay operator $\mathcal{D}$ describing leaky integration, a jump operator $\mathcal{J}$ accounting for the action of synaptic inputs, and a threshold-and-reset operator $\mathcal{T}$ that implements spike generation upon threshold crossing. These matrices are computed in discrete time with time bin $\Delta t$. The discrete decay operator $\mathcal{D}$ acting on the probability density vector $\boldsymbol{p}$ is derived from its continuous counterpart D, using equations (22) and (23) (explained graphically in figure 3), as follows

$$\begin{align} \label{D_integral} (\mathcal{D}\boldsymbol{p})_i\Delta V= \int^{v_{i+1}}_{v_i} DP(V)\,{\rm d}V = \int^{v_{i+1}}_{v_i} \underbrace{{\rm e}^{\Delta t/\tau_m}}_{\eta}P({\rm e}^{\Delta t/\tau_m}V)\,{\rm d}V \nonumber \end{align} \tag{ 24 }$$

and with change of variables $x=\eta V$ this leads to

$$\begin{align*} \int^{\eta v_{i+1}}_{\eta v_i} P(x)\,{\rm d}x&= \sum_{j=\lceil \eta i\rceil}^{\lceil\eta (i+1) \rceil-1} p_j.\Delta V+ \Big{(} \Delta V \lceil \eta i \rceil-\eta v_i \Big{)} p_{\lfloor \eta i\rfloor}- \Big{(}\Delta V \lceil \eta (i+1)\rceil -\eta v_{i+1}\Big{)} p_{\lfloor \eta (i+1)\rfloor}\nonumber \\ &=\Delta V \Big{[}\sum_{j=\lceil \eta i\rceil}^{\lceil\eta (i+1) \rceil-1} p_j+ \Big{(} \lceil \eta i \rceil-\eta i \Big{)} p_{\lfloor \eta i\rfloor}- \Big{(} \lceil \eta (i+1)\rceil-\eta (i+1)\Big{)} p_{\lfloor \eta (i+1)\rfloor}\Big{]} \nonumber \end{align*}$$

where drift operator D is defined as $DP(x) =\underbrace{{\rm e}^{\Delta t/\tau_m}}_\eta P({\rm e}^{\Delta t/\tau_m}x)$. This definition leads to the following decay matrix [25]

$$\begin{align} \label{D_matrix} \mathcal{D}_{ij}= \sum_{k=\lceil \eta i\rceil}^{\lceil \eta(i+1) \rceil-1} \delta_{j,k}+ \Big{(} \lceil \eta i \rceil-\eta v_i \Big{)} \delta_{j,{\lfloor \eta_i\rfloor}}- \Big{(} \lceil \eta (i+1)\rceil-\eta (i+1)\Big{)} \delta_{j,{\lfloor \eta(i+1)\rfloor}}. \nonumber \end{align} \tag{ 25 }$$

The jump distribution, which underlies the jump operator $\mathcal{J}$, can be derived from the count distribution of excitatory and inhibitory synaptic events in a small time interval $\Delta t$. Assuming that they arrive with Poisson statistics (maximum entropy) at a rate of $r_{{\rm ex}}$ and $r_{{\rm in}}$, respectively, we obtain

$$\begin{align} \label{eq_P_q_r} q(m,n)= {\rm e}^{-(a+b)}\frac{a^m \: b^n}{m! \: n!} \nonumber \end{align} \tag{ 26 }$$

where $a= r_{{\rm ex}}\Delta t$ and $b= r_{{\rm in}}\Delta t$ are the corresponding expected event counts. Dummy indices m and n correspond to excitatory and inhibitory counts. The jump distribution of the membrane potential is given by

$$\begin{align} \label{eq_jump_dist} Q(\gamma)=\sum_{m,n=0}^{\infty}q(m,n)\delta\big{(}\gamma-(m-gn)h\big{)} \nonumber \end{align} \tag{ 27 }$$

where the list of points {γ} corresponding to excitatory and inhibitory jumps mentioned in equation (26) is the set of all possible transitions $\{(m, n) \}\rightarrow \{\gamma_{\kappa}\}$ with $\kappa(m, n)=\lfloor \frac{(m-gn)h}{\Delta V} \rfloor$. The jump operator $\mathcal{J}$, that acts on the probability density vector $\boldsymbol{p}$, is then derived as

$$\begin{align*} (\mathcal{J}\boldsymbol{p})_{i}\Delta V=\int^{v_{i+1}}_{v_i}{\rm d}V\int_{V_-}^{{V_+}} {\rm d}\gamma Q(\gamma) P(V-\gamma) =\int^{v_{i+1}}_{v_i}{\rm d}V\sum_{j}\int_{v_i-j}^{v_{i-j+1}} {\rm d}\gamma Q(\gamma) P(V-\gamma) \nonumber \end{align*}$$

where V₋ and V₊ are lower and upper bounds of the voltage space. Here we denote the shift in the discrete index as $\kappa=\lfloor \frac{\gamma}{\Delta V} \rfloor$ with a jump matrix given by

$$\begin{align} \mathcal{J}_{ij}=\sum_{m,n} q(m,n)\delta_{j,i-\kappa(m,n)} \rightarrow \sum_{\kappa} c(\kappa) \delta_{j,i-\kappa} = c(i-j). \nonumber \end{align} \tag{ 28 }$$

where κ is implicitly a function of (m, n). We note that there may be multiple pairs (m, n) that lead to same κ, hence

$$\begin{align} c(\kappa)=\sum_{\{m(\kappa),n(\kappa)\}} q(m,n). \nonumber \end{align} \tag{ 29 }$$

For homogenous jump distribution we study in this work, the discretization of the jump operator is simply a map of continuous convolution to discretized convolution. Hence in principle $\mathcal{J}_{ij}$ is a cyclic matrix, although this property is violated at V₋ and V₊ boundaries where one needs to truncate the discretization. The approximation is precise as P(V) vanishes fast with application of the decay operator in vicinity of V  <  V₋ and similarly around V  >  V₊ by combined effect of reset and jump operators. Note that the truncation point $V_+>(V_{{\rm th}}+h_{{\rm ex}})$, as for finite $\Delta t$ multiple excitatory spikes are possible. Moreover, by selecting the precision of the truncated jump distribution we determine jump numbers n_jump (for details look at [25]), that is typically increasing with excitatory rate $r_{{\rm ex}}$. (This implies that the size of $\mathcal{J}$ is then $(N+n_{\rm jump})\times N$.) The truncation error bound doesn't lead to large n_jump because the jump distribution decays exponentially fast for large number of spikes (see equation (26)).

The threshold-and-reset operator $\mathcal{T}$ takes all the states above threshold and inserts them at the reset potential. This is simply given as

$$\begin{align} \mathcal{T}_{ij}={\bf 1}_{j> \frac{V_{{\rm th}}}{\Delta V}} \delta _{i,\frac{V_{\rm r}}{\Delta V}}+{\bf 1}_{j\leqslant \frac{V_{{\rm th}}}{\Delta V}} \delta _{i,j} \label{eq_threshold_matrix} \nonumber \end{align} \tag{ 30 }$$

where ${\bf 1}$ is a vector with ${\bf 1}_i=1$ and the size of the threshold matrix $\mathcal{T}$ is $N\times(N+n_{\rm jump})$. We note that the matrix multiplication $\mathcal{T} \mathcal{J}$ results in an $N\times N$ matrix. Finally, the time evolution matrix of the corresponding Markov chain is given as a product of the three operators described above

$$\begin{align} \mathcal{M}=\mathcal{T}\mathcal{J}\mathcal{D} \label{eq_M1d} \nonumber \end{align} \tag{ 31 }$$

where $\mathcal{M}$ is a matrix with dimension $N\times N$. The discrete stationary distribution is

$$\begin{align} \mathcal{M}\boldsymbol{p}_0=\boldsymbol{p}_0, \nonumber \end{align} \tag{ 32 }$$

and the corresponding stationary rate is given as

$$\begin{align} \label{eq_rate_markov} r_{{\rm Markov}}=\Big{(}\tau_{{\rm ref}}+\Big{[}\frac{\Delta V}{\Delta t} \sum_{i \geqslant i_{{\rm th}}} (\mathcal{J}\mathcal{D}\boldsymbol{p}_0)_i\Big{]}^{-1}\Big{)}^{-1}. \nonumber \end{align} \tag{ 33 }$$

For a given $\Delta t$, the precision of the numerical approximation depends on the order N as discussed in [25] and as shown here in figures 11 and 12. We discuss effect of the grid further in section 4.

The cross-covariance function of two stationary spike trains $S_a(t) = \sum\nolimits_l \delta(t-t^a_l)$ (a  =  1, 2) is defined as

$$\begin{align} C_{12} (\tau)= \langle S_1(t+\tau) S_2(t) \rangle -\langle S_1(t+\tau)\rangle \langle S_2(t) \rangle \nonumber \end{align} \tag{ 34 }$$

where $\langle S_a(t)\rangle = r_a$, with $\langle \cdot\rangle$ indicating the ensemble average. For the model we studied here stationary rates were computed using equation (35). Similarly the cross-covariance function can be expressed in terms of the conditional firing rate $r_{1\vert 2}(\tau)$, two stationary rates r₁ and r₂, and the amplitude r₀ of a δ-function

$$\begin{align} \label{eq_corr_rate} C_{12} (\tau)= r_0\delta(\tau)+r_2 (r_{1\vert 2}(\tau)-r_1). \nonumber \end{align} \tag{ 35 }$$

Given the spiking neuron model considered here, equation (35) can be derived from the stationary joint membrane potential distribution $P_{0}(V_1, V_2)$. Conditional on a spike at time t₀ in the first neuron, one has to find the instantaneous distribution of the membrane potential of the first neuron

$$\begin{align} P_{1\vert 2}(V_1;t_0){\rm d}V\equiv {\rm Prob}\big{(}{\rm{~1st~neuron~in~}} [V_1,V_1+{\rm d}V_1] \: \big{\vert} \: {\rm{~2nd~neuron~spikes~at}} t_0\big{)} \nonumber \end{align} \tag{ 36 }$$

and the same for both neurons changing roles. The probability of observing a consecutive spike is given by the (hypothetical) continuous flux density $F(V_1, V_{{\rm th}, 2};t_0)$ with the normalization in equation (37). The conditional flux is computed for the Markov approximation later in equation (65). We present here only its relation to correlation function equation (35). In general we simply compute the exit rate at the threshold distributed over V₁ by solving the following initial value problem

$$\begin{align} \label{eq_IVP} P_{1\vert 2}(V_1;t_0)=F(V_1,V_{{\rm th},2};t_0)\Big{(}\int_{\infty}^{V_{{\rm th},1}} F(x,V_{{\rm th},2};t_0) \,{\rm d}x\Big{)}^{-1} \nonumber \end{align} \tag{ 37 }$$

$$\begin{align} P_{1\vert 2}(V_1;t_0+\Delta t)=M_1 P_{1\vert 2}(V_1;t_0). \nonumber \end{align} \tag{ 38 }$$

where M₁ is a (rather hypothetical) discrete time ($\Delta t$) evolution operator of the neuron model in equations (92) and (93). The set of coupled jump equations, corresponding to equations (92) and (93), may be solvable only in discrete space. The discretization of M₁ is given by equation (31) as $\mathcal{M}_1$. $\mathcal{M}_1$ leads to a (forward) time evolution in steps of $\Delta t$ of $\boldsymbol{p}$. We note that if the renewal process (e.g. output spike trains) is stationary and in equilibrium we can select t₀  =  0 without loss of generality. Yet applying the space discretization the linear time evolution is then

$$\begin{align} \label{eq_IVP_desc} \boldsymbol{p}_{1\vert 2}(0)=\frac{1}{\Delta V}\int_{v_{i}}^{v_{i+1}}{\rm d}xP_{1\vert 2}(V;0) \nonumber \end{align} \tag{ 39 }$$

$$\begin{align} \boldsymbol{p}_{1\vert 2}(t_n)=\frac{1}{\Delta V}\int_{v_{i}}^{v_{i+1}}{\rm d}x M_1^{n} P_{1\vert 2}(x,0)=\mathcal{M}_1^{n}\boldsymbol{p}_{1\vert 2}(0) \nonumber \end{align} \tag{ 40 }$$

The instantaneous conditional rate $r_{1\vert 2}(t)$ in equation (35) is computed with a formula similar to equation (33)

$$\begin{align} \label{eq_rate_markov_t} r_{1\vert 2}(t_n)=r_{{\rm Markov}}(t_n)=\Big{(}\tau_{{\rm ref}}+\Big{[}\frac{\Delta V}{\Delta t} \sum_{i \geqslant i_{{\rm th}}} (\mathcal{J}\mathcal{D}\boldsymbol{p}_{1\vert 2}(t_n))_i\Big{]}^{-1}\Big{)}^{-1}. \nonumber \end{align} \tag{ 41 }$$

Finally, the covariance function $C_{21} (\tau)$ in equation (35) is derived by using $r_{1\vert 2}(t)$ and two stationary rates r₁ and r₂ and r₀. Note that the method described here can be generalized to higher order correlations as well.

We also demonstrate here that with the parametrization in section 2.1 the desired correlation coefficients c is achieved. For convenience, we define C[.,.] and A[.] as integrals ($\int_{-\infty}^{\infty}d \tau...$) of the covariance and auto-covariance functions respectively. The correlation coefficient of two input terms in equations (92) and (92) is then

$$\begin{align*} C_{\rm Poisson}(c)&=\frac{C\Big{[}h_{{\rm ex},1}(S^{{\rm ex}}_1+S^{{\rm ex}}_c)+h_{{\rm in},1}(S^{{\rm in}}_1+S^{{\rm in}}_c), h_{{\rm ex},2}(S^{{\rm ex}}_2+S^{{\rm ex}}_c)+h_{{\rm in},2}(S^{{\rm in}}_2+S^{{\rm in}}_c)\Big{]}} {\sqrt{A\Big{[}(h_{{\rm ex},1}(S^{{\rm ex}}_1+S^{{\rm ex}}_c)+h_{{\rm in},1}(S^{{\rm in}}_1+S^{{\rm in}}_c)\Big{]} A\Big{[} h_{{\rm ex},2}(S^{{\rm ex}}_2+S^{{\rm ex}}_c)+h_{{\rm in},2}(S^{{\rm in}}_2+S^{{\rm in}}_c)\Big{]}}}\nonumber \\ &=\frac{h_{1}h_{2} (r^{{\rm ex}}_c+g^2r^{{\rm in}}_c)} {\sqrt{\Big{(}h_{1}^2[(r^{{\rm ex}}_1+r^{{\rm ex}}_c)+g^{2}(r^{{\rm in}}_1+r^{{\rm in}}_c)]\Big{)} \Big{(} h_{2}^2[(r^{{\rm ex}}_2+r^{{\rm ex}}_c)+g^{2}(r^{{\rm in}}_2+r^{{\rm in}}_c)]\Big{)}}} \nonumber \\ &=\frac{h_{1}h_{2} c\phi} {\sqrt{\Big{(} h_{1}^2[(r^{{\rm ex}}_1+g^{2}r^{{\rm in}}_1)+(r^{{\rm ex}}_c+g^{2}r^{{\rm in}}_c)]\Big{)} \Big{(} h_{2}^2[(r^{{\rm ex}}_2+g^{2}r^{{\rm in}}_2)+(r^{{\rm ex}}_c+g^{2}r^{{\rm in}}_c)]\Big{)}}} \nonumber \\ &=\frac{c\phi} {\sqrt{\Big{(} (1-c)\phi+c\phi\Big{)} \Big{(} (1-c)\phi+c\phi\Big{)}}} \nonumber \\ &=\frac{c\phi} {\sqrt{\phi \phi}} =c \nonumber \end{align*}$$

where we plugged in $(r^{{\rm ex}}_a+g^{2}r^{{\rm in}}_a)=(1-c)\phi$ and $(r^{{\rm ex}}_c+g^{2}r^{{\rm in}}_c)=c\phi$, which are derived from formulas equations (9)–(19).

We now use a 2-dimensional state space describing the joint membrane potential evolution of two neurons. Correlated and uncorrelated Poisson jumps push the 2D membrane potential vector $(V_1, V_2)$ into three independent directions, (1, 0), (0, 1) and (1, 1), allowing jumps in the positive (excitatory) and negative (inhibitory) direction, respectively. Hence, the jump distribution from an initial position $(V_1-\gamma_1, V_2-\gamma_2)$ to a new position $(V_1, V_2)$ in state space driven by 6 independent Poisson processes is obtained from a 2D convolution

$$\begin{align} Q(\gamma_1,\gamma_2)=\int Q_{1}(\gamma_1-Z)Q_{2}(\gamma_2-Z)Q_{c}(Z)\,{\rm d}Z. \nonumber \end{align} \tag{ 42 }$$

The equation above is valid for homogenous jump distribution we study here, where we simply omitted the variables indicating the initial point $(V_1-\gamma_1, V_2-\gamma_2)$ for the sake of compact notation. This can easily be generalized, for jump distributions that depend on V. Inserting the expressions for the jump distributions in each direction, collecting all terms and using equation (26), we obtain

$$\begin{align} Q(\gamma_1,\gamma_2)=&\int \,{\rm d}Z \sum_{i,j=0}^{\infty}\sum_{k,l=0}^{\infty}\sum_{m,n=0}^{\infty} q_{1}(i,j) q_{2}(k,l) q_{c}(m,n) \delta\big{(}\gamma_1-Z-(i-gj)h\big{)} \nonumber \\ &\hspace{3cm}\times \delta\big{(}\gamma_2-Z-(k-gl)h\big{)}\delta\big{(}Z-(m-gn)h\big{)}.\nonumber \end{align} \tag{ 43 }$$

This expression is also valid for more general input statistic models that rely on a decomposition into statistically independent parts [31, 33]. The corresponding amplitude distribution q(m, n) can cover cases of non-Poissonian spike trains with nontrivial higher-order statistics, e.g. MIP or SIP processes [33]. Here we consider the shared Poisson input model as described by equations (92) and (93). Integrating the expression with respect to Z and inserting the mean event counts $a_i= r_{{\rm ex}, i}\Delta t$ and $b_i= r_{{\rm in}, i}\Delta t$, the resulting expression is given in a compact form

$$\begin{align} \nonumber Q(\gamma_1,\gamma_2)=&\sum_{i,j=0}^{\infty}\sum_{k,l=0}^{\infty}\sum_{m,n=0}^{\infty} {\rm e}^{-(a_1+b_1+a_2+b_2+a_s+b_s)} \frac{a_1^i \: b_1^{\,j}}{i! \: j!} \frac{a_2^k \: b_2^l}{k! \: l!} \frac{a_s^m \: b_s^n}{m! \: n!} \: \nonumber \\&\hspace{3cm}\times\delta\big{(}\gamma_1-(i-m-g(\,j-n))h\big{)}\:\delta\big{(}\gamma_2-(k-m-g(l-n))h\big{)}. \label{P_jump_2d} \nonumber \end{align} \tag{ 44 }$$

We can simplify this expression by choosing a regular grid according to $h_e= n \Delta V$ and $h_i= m \Delta V$, for integers m and n. In practice, however, the resulting sum will be truncated for a given tolerance. We will use the matrix form of the discretized operators in the following section, which is equivalent to the above expression.

Here we discuss the action of operators on state vectors of the discretized $(V_1, V_2)$ space, assuming N bins in each dimension. We write the stationary voltage distribution in the form

$$\begin{align} \label{eq_p_0} \boldsymbol{p}_0(X,Y) \equiv \sum_{i,j=1}^N \Omega_{ij} \: X_i\otimes Y_j. \nonumber \end{align} \tag{ 45 }$$

where X and Y are two suitable orthogonal bases of $\mathbb{R}^N$. We will define a specific choice for X and Y later in section 2.6. As there is no crosstalk between the two neurons except by shared input leading to a correlated jump distribution, threshold and decay operators (tensors) are separable

$$\begin{align} \mathcal{D}_{\rm 2D}=\mathcal{D}_1\otimes \mathcal{D}_2 \nonumber \end{align} \tag{ 46 }$$

$$\begin{align} \mathcal{T}_{\rm 2D}=\mathcal{T}_1\otimes \mathcal{T}_2. \nonumber \end{align} \tag{ 47 }$$

Separability would also apply to the jump distribution for an input correlation coefficient of c  =  0, corresponding to two neurons without shared input

$$\begin{align} \mathcal{J}_{\rm 2D}=\mathcal{J}_1\otimes \mathcal{J}_2. \nonumber \end{align} \tag{ 48 }$$

Regardless of the method, once we have computed the correct Markov matrix for 2D via

$$\begin{align} \mathcal{M}_{\rm 2D}=\mathcal{T}_{\rm 2D} \: \mathcal{J}_{\rm 2D} \: \mathcal{D}_{\rm 2D} \label{eq_M2d} \nonumber \end{align} \tag{ 49 }$$

one can also find the stationary joint membrane potential distribution as the Perron–Frobenius eigenvector $\boldsymbol{p}_0$ of the propagator matrix $\mathcal{M}_{\rm 2D}$

$$\begin{align} \mathcal{M}_{\rm 2D} \boldsymbol{p}_0=\boldsymbol{p}_0. \label{eq_full_markov} \nonumber \end{align} \tag{ 50 }$$

Regarding the correlated jump distribution there are two ways of constructing $\mathcal{J}$ operators. One possibility is described in equation (44). Alternatively, we construct linear Markov jump operators in 2D exploiting the commutativity of infinitesimal operators

$$\begin{align*} \mathcal{J}_{\rm 2D}={\rm e}^{(J_{1} \otimes I +I\otimes J_{2}+J_{c} \otimes J_{c})}. \nonumber \end{align*}$$

Here I is the identity matrix. Using the properties of the operator product $\otimes$ and commutativity of the individual factors, we can simplify this expression

$$\begin{align*} \mathcal{J}_{\rm 2D}=({\rm e}^{J_{1}} \otimes I)(I \otimes {\rm e}^{J_{2}}) ({\rm e}^{J_{c} \otimes J_{c}}). \nonumber \end{align*}$$

In order to expand the third term we define U and D operators as up and down transition matrices, where U corresponds to a 1-step up transition, and D corresponds to a 1-step down transition. This leads to

$$\begin{align} U^{k}_{ij}=\delta_{i,j-k} \quad{\rm{and}}\quad D^{l}_{ij}=\delta_{i,j+l}. \nonumber \end{align} \tag{ 51 }$$

Hence, discrete approximations to infinitesimal generators of private components are given as

$$\begin{align} J_i=a_iU^{k}+b_iD^{l} \nonumber \end{align} \tag{ 52 }$$

where index i  =  {1, 2, c} and matrix powers k and l are defined on $h_e= k \Delta V$ and $h_i= l \Delta V$ for simplicity. (This restrictive assumption can be generalized easily by computing the continuous jump distribution in equation (27) 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

$$\begin{align} q_{c}(m,n)={\rm e}^{-(a_c+b_c)}\frac{1}{m!}\frac{1}{n!}a_c^mb_c^{n} \nonumber \end{align} \tag{ 53 }$$

for m excitatory and n inhibitory jumps. Expanding yields

$$\begin{align} \label{eq_J_2D_1} \mathcal{J}_{\rm 2D}=\sum_{m=0} \sum_{n=0}q_{c}(m,n) [{\rm e}^{J_1}U^{km} D^{ln}\otimes {\rm e}^{J_2}U^{km} D^{ln}]. \nonumber \end{align} \tag{ 54 }$$

As we noted before, the above construction is quite general and can be applied easily for general amplitude distributions. We only need to consider a discrete amplitude distribution q_c(m, n) in a given time bin of size $\Delta t$, as described above, see also [33, 34].

A final remark on the method described in this section concerns the commutativity of matrices. This property leads to a numerically more efficient expression

$$\begin{align} \label{eq_J_2D} \mathcal{J}_{\rm 2D}=\sum_{j \: \in \: \mathbb{J}} c(\,j) [{\rm e}^{J_1}O^{\,j}\otimes {\rm e}^{J_2}O^{\,j}] \nonumber \end{align} \tag{ 55 }$$

where the set of integers $\mathbb{J}$ is defined as list of all jump numbers j  =  mk  −  nl. The coefficient

$$\begin{align} c(\,j)= \sum_{m,n} q_{c}(m,n)\delta_{j,mk-nl} \nonumber \end{align} \tag{ 56 }$$

is the probability of j jumps. The jump generator O is then defined in terms of matrix powers

$$\begin{align} O^{\,j}=\left\{\begin{array}{@{}ll@{}} U^{\,j}, & j>0 \nonumber \\ I,& j=0 \nonumber \\ D^{\,j}& j<0 \nonumber \\ \end{array}\right.. \nonumber \end{align} \tag{ 57 }$$

We emphasize that in case heterogenous PSPs the grid can be tuned to satisfy e.g. $\frac{h_{{\rm ex}, 1}}{{\rm d}V_1}=\frac{h_{{\rm ex}, 2}}{{\rm d}V_2}$ so that the shared spikes always generate a diagonal (j, j) transition in the index space. In general, in case of two distinct jump distributions driven by shared spikes in V₁ and V₂ directions, we can rescale the space w.r.t. jump distributions so that equation (55) holds. Additionally, we note that we assumed $h= k {\rm d}x$ and $gh= l {\rm d}x$ for convenience. We can relax this assumption by using full 1D jump matrices (in homogenous case a cyclic matrix) instead of matrix exponential ${\rm e}^{J}$ and 1D jump distribution instead of c(j) in equation (55). We confirm that for a fine grids the corresponding difference is marginal.

The expansion method described above is straightforward, but rather cumbersome to implement. We will now introduce a basis for the expansion of correlated jump operators suitable to reduce the cost of the computations involved, and helpful to increase the accuracy of a truncation. With our method, as compared to others, we have to solve linear equations of lower dimensionality in order to get a better approximation for the correlation function. Some further arguments for selecting Singular Value Decomposition are discussed in the results section. SVD of a matrix is given as

$$\begin{align} \mathcal{M}=\mathcal{L}\Sigma \mathcal{R}^T \nonumber \end{align} \tag{ 58 }$$

where the diagonal entries of the diagonal matrix $\Sigma$ are the square roots of the non-zero eigenvalues of $\mathcal{M}\mathcal{M}^T$ and $\mathcal{M}^T \mathcal{M}$. Both matrices $\mathcal{R}$ and $\mathcal{L}$ are orthogonal with columns consisting of the eigenvectors of $\mathcal{M}\mathcal{M}^T$ and $\mathcal{M}^T \mathcal{M}$, respectively. We show a 2D example SVD of Markov matrix in figure 2. For $\mathcal{M}$ replaced by a single-neuron time evolution matrix $\mathcal{M}_1$ ($\mathcal{M}_2$), we define the matrix X (Y) by the selected orthogonal subspace of dimension K (L) of $\mathcal{R}_1$ ($\mathcal{R}_2$), according to the largest K (L) singular values, resepectively. In order to project $\mathcal{J}_1$ ($\mathcal{J}_2$)and $\mathcal{J}_2$ ($\mathcal{T}_2$) to this subspace, we also extend the orthogonal basis X and Y to supra-threshold transitions

$$\begin{align} \tilde{X}= \left(\begin{array}{@{}ll@{}} X & 0 \nonumber \\ 0 & I_{m} \end{array}\right), \qquad \tilde{Y}= \left(\begin{array}{@{}ll@{}} Y & 0 \nonumber \\ 0 & I_{n} \end{array}\right), \nonumber \end{align} \tag{ 59 }$$

where X is an $M\times K$ and Y is an $N\times L$ operator, respectively. I_m and I_n are identity matrices, where m and n are the maximal supra-threshold jump numbers induced by both private and shared inputs (here we drop the index 'jump' in n_jump and m_jump). For a detailed discussion of maximum jump numbers, denoted here as m and n, check [25]. The combined action of $\mathcal{J}$ and $\mathcal{T}$ for c  =  0 is then expressed as reduced operators

$$\begin{align} \label{eq_reduced_operators} \mathcal{T}_1 J_1 \rightarrow X^T \mathcal{T}_1\tilde{X} \tilde{X}^T \mathcal{J}_1 X, \qquad \mathcal{T}_2 J_2 \rightarrow Y^T \mathcal{T}_2 \tilde{Y} \tilde{Y}^T \mathcal{J}_2 Y, \nonumber \end{align} \tag{ 60 }$$

which map $M\times M$ ($N\times N$) to $K \times K$ ($L \times L$) matrices. Below we use the same dimensional reduction for correlated operators. In equation (55) the correlated jump operators are expressed as

$$\begin{align*} \mathcal{J}_{\rm 2D}= \sum_{j \: \in \: \mathbb{J}}c(\,j) \: A_j \otimes B_j. \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

A dimensional reduction is then achieved by using equation (60) in

$$\begin{align*} \mathcal{M}_{\rm 2D}\boldsymbol{p}_0= \sum_{j \: \in \: \mathbb{J}} c(\,j) \: (\mathcal{T}_1 A_j \mathcal{D}_1\otimes \mathcal{T}_2 B_j\mathcal{D}_2) \boldsymbol{p}_0. \nonumber \end{align*}$$

In order to find $\boldsymbol{p}_0$ in equation (45) we need to solve equation (50), which reads

$$\begin{align*} \mathcal{M}_{\rm 2D}\boldsymbol{p}_0=\boldsymbol{p}_0. \nonumber \end{align*}$$

The projection operators $X^T\otimes Y^T$ satisfy $(X \otimes Y) (X^T \otimes Y^T)\boldsymbol{p}_0=\boldsymbol{p}_0$. Applying them to the left hand side of this equation equation (50), we obtain

$$\begin{align*} (X^T \otimes Y^T)\mathcal{M}_{\rm 2D}&=(X^T \otimes Y^T)\mathcal{M}_{\rm 2D}(X \otimes Y) (X^T \otimes Y^T)\boldsymbol{p}_0\nonumber \\ &=\sum_{j \: \in \: \mathbb{J}} c(\,j) \:\Big{(}(X^{T}\mathcal{T}_1\tilde{X} \tilde{X}^T A_j X X^T\mathcal{D}_1 X) \otimes (Y^T \mathcal{T}_2 \tilde{Y} \tilde{Y}^T B_j Y \mathcal{D}_2 Y)\Big{)} \Omega\nonumber \\ &\mathop{=}\limits^{\rm RHS}\Omega. \nonumber \end{align*}$$

$\Omega\equiv (X^T \otimes Y^T)\boldsymbol{p}_0$ can be expressed in a simpler form as

$$\begin{align*} (X^T\otimes Y^T) \boldsymbol{p}_0&=(X^T\otimes Y^T) (\sum_{ij} \Omega_{ij}X_i \otimes Y_j)\nonumber \\ &=\sum_{ij}\Omega_{ij} X^TX_i \otimes Y^T Y_j)\nonumber \\ &=\sum_{ij}\Omega_{ij} e_i \otimes e_j\equiv \Omega \nonumber \end{align*}$$

for monomials $(e_i)_k=\delta_{i, k}$. The reduced equation is then

$$\begin{align} \label{eq_reduced_eigen} \mathcal{Q} \Omega=\Omega \nonumber \end{align} \tag{ 61 }$$

where $\mathcal{Q}$ is a tensor defined as

$$\begin{align} \label{eq_Q} \mathcal{Q}=\sum_{j \: \in \: \mathbb{J}} c(\,j)(X^{T}\mathcal{T}_1\tilde{X} \tilde{X}^T A_j X X^T\mathcal{D}_1 X) \otimes (Y^T \mathcal{T}_2 \tilde{Y} \tilde{Y}^T B_j Y \mathcal{D}_2 Y^T Y). \nonumber \end{align} \tag{ 62 }$$

The dimensionally reduced problem in equation (61) can then be solved in practice by reindexing tensor indices $(i, j, k, l) \mapsto (I, K)$.

Using the decomposed 2D stationary distribution obtained by reduction, one can compute the flux distribution with the help of matrix operators. We compute the flux distribution using the 2D decay and jump operators with thresholding imposed only at one of the boundaries. A scheme illustrating the situation is shown in figure 1(c). Here we explain how to obtain the flux distribution conditional on a spike in one neuron. In the general case of correlated neurons, the action of the full operators $\mathcal{J}_{\rm 2D}$ and $\mathcal{D}_{\rm 2D}$ is given as

$$\begin{align} \tilde{\boldsymbol{p}}_{{\rm J},kl}=\mathcal{J}_{\rm 2D}\mathcal{D}_{\rm 2D}\boldsymbol{p}_0=\sum_{k^{\prime},l^{\prime}} \sum_{j \: \in \: \mathbb{J}} c(\,j) (A_j \mathcal{D}_1\otimes B_j \mathcal{D}_2)_{kl,k^{\prime}l^{\prime}} (X^T_{k^{\prime}} \Omega Y_{l^{\prime}}), \nonumber \end{align} \tag{ 63 }$$

where $\tilde{\boldsymbol{p}_{{\rm J}}}$ is a $(M+m)\times(N+n)$ matrix. The implicitly summed subspace indices, $\boldsymbol{p}_{0, kl}=\sum\nolimits_{m, n}X^T_{k^{\prime}, m} \Omega_{mn} Y_{n, l^{\prime}}$ with dummy indices m and n, are not shown and, A_j, B_j and c(j) are defined in equation (55). This equation can be written in a concise form with implicit indices as

$$\begin{align*} \tilde{\boldsymbol{p}}_{{\rm J}}= \sum_{j \: \in \: \mathbb{J}} c(\,j) A_{j}^{T} \mathcal{D}^T_{1}X^T\Omega Y \mathcal{D}_{2} B_{j}=\tilde{X}^{T}\tilde{\Omega} \tilde{Y}. \nonumber \end{align*}$$

Again, for practical reasons, computations were reduced by projecting $\mathcal{J}_{\rm 2D}\mathcal{D}_{\rm 2D}$ to extended subspaces $\tilde{X}$ and $\tilde{Y}$ similar to equation (60)

$$\begin{align} \tilde{\Omega} =\sum_{j \: \in \: \mathbb{J}} c(\,j) \:\Big{[}(\tilde{X}^T A_j X X^T\mathcal{D}_1 X) \otimes (\tilde{Y}^T B_j Y \mathcal{D}_2 Y)\Big{]} \Omega. \nonumber \end{align} \tag{ 64 }$$

In order to compute the probability of jumps we need to sum the probabilities for a jump over the threshold $V_1>V_{{\rm th}, 1}$ or $V_2>V_{{\rm th}, 2}$. The flux distribution is obtained as

$$\begin{align} \label{eq_flux} \boldsymbol{f}_{1\vert 2,k} \propto \sum_{l=N}^{N+n} \tilde{\boldsymbol{p}}_{{\rm J},kl}=\sum_{l=N}^{N+n} \tilde{X}^{T}_k\tilde{\Omega}\tilde{Y}_l, \qquad \boldsymbol{f}_{2\vert 1,l}\propto \sum_{k=M}^{M+m} \tilde{\boldsymbol{p}}_{{\rm J},kl}= \sum_{k=M}^{M+m}\tilde{X}^{T}_k\tilde{\Omega}\tilde{Y}_l \nonumber \end{align} \tag{ 65 }$$

defined for k  <  M and l  <  N. Here we have $\tilde{X}^{T}_k\tilde{\Omega}\tilde{Y}_l=\sum\nolimits_{ij}\tilde{X}_{k, j}\tilde{\Omega}_{ij}\tilde{Y}_{j, l}$ for the sake of notation. These expressions are normalized such that $\sum\nolimits_{k} \boldsymbol{f}_{1\vert 2, k} \Delta x=1$ and $\sum\nolimits_{l}\boldsymbol{f}_{2\vert 1, l} \Delta y=1$. We denote the procedure above symbolically as

$$\begin{align} \label{flux_op} \boldsymbol{f}_{1\vert 2}=\mathcal{F}_{1\vert 2}(\boldsymbol{p}_0) \nonumber \end{align} \tag{ 66 }$$

$$\begin{align} \boldsymbol{f}_{2\vert 1}=\mathcal{F}_{2\vert 1}(\boldsymbol{p}_0). \nonumber \end{align} \tag{ 67 }$$

The amplitude of the delta singularity at the reset potential is obtained as

$$\begin{align} \label{eq_delta_amp} r_0=\frac{(\Delta V)^2}{\Delta t}\sum_{l=N}^{N+n} \sum_{k=M}^{M+m} \tilde{X}^{T}_k\tilde{\Omega}\tilde{Y}_l, \nonumber \end{align} \tag{ 68 }$$

where $\Delta V$ and $\Delta t$ are space and time bins consecutively. Once we have found the conditional flux distribution, we can solve the initial value problem defined in equation (37) in order to determine the conditional rates $r_{1\vert 2}(t)$ or $r_{2\vert 1}(t)$.

We compare the exit rate of the stochastic system with post-synaptic potentials of finite amplitude with the analytic result obtained for the diffusion approximation [16]. For small enough PSPs the difference in rates of the two models is small

$$\begin{align} \label{eq_rates} r_{sg}^{-1}=\tau_{{\rm ref}}+\tau_{m}\sqrt{\pi} \int_{\frac{V_{\rm r}-\mu}{\sigma}}^{\frac{V_{{\rm th}}-\mu}{\sigma}} {\rm e}^{x^2}[1+{\rm erf}(x)] \,{\rm d}x\nonumber \end{align} \tag{ 69 }$$

$$\begin{align} r_{{\rm Markov}}=\Big{(}\tau_{{\rm ref}}+\Big{[}\frac{\Delta V}{\Delta t} \sum_{i \geqslant i_{{\rm th}}} (\mathcal{J}\mathcal{D}\boldsymbol{p}_0)_i\Big{]}^{-1}\Big{)}^{-1}. \nonumber \end{align} \tag{ 70 }$$

We use the absolute difference between the two rates

$$\begin{align} r_{{\rm error}}(\mu,\sigma)= \vert r_{sg}(\mu,\sigma)-r_{{\rm Markov}}(\mu,\sigma)\vert \nonumber \end{align} \tag{ 71 }$$

to account for the accuracy of a specific space-time grid $(\Delta V, \Delta t)$.

The correlation coefficient in figure 11(a) is computed with the formula

$$\begin{align} \label{eq_corr_coeff} C_{\rm out}(c)=\frac{r_0+\int r_2 (r_{1\vert 2}(\tau)-r_1)\,{\rm d}\tau +\int r_1 (r_{2\vert 1}(\tau)-r_2)\,{\rm d}\tau}{CV_1 CV_2\sqrt{r_1r_2}}. \nonumber \end{align} \tag{ 72 }$$

We used (69) for the stationary rate r, r₀ is the amplitude of the δ-function in equation (68), and the coefficient of variation, $CV^2=\frac{\sigma^2_{ISI}}{\mu^2_{ISI}}$, is computed with the equation

$$\begin{align} CV^{2}=2\pi r^2 \int_{\frac{V_{\rm r}-\mu}{\sigma}}^{\frac{V_{{\rm th}}-\mu}{\sigma}} {\rm e}^{x^2}\,{\rm d}x \int^{y}_{-\infty}[1+{\rm erf}(x)]^2 \,{\rm d}y \nonumber \end{align} \tag{ 73 }$$

as given in [13]. The result of linear response theory for output spike train correlations is given in [35] as

$$\begin{align} \label{eq_Cout_pert} C_{{\rm out},{\rm pert}}^{(1)}=\frac{c\sigma_1 \sigma_2 \frac{dr_1}{{\rm d}\mu}\frac{dr_2}{{\rm d}\mu}}{{\rm CV}_1 {\rm CV}_2\sqrt{r_1r_2}} =\frac{c \: r_1^{3/2} r_2^{3/2} [{\rm e}^{x_t^2} {\rm erf}(x_t)-{\rm e}^{x_r^2} {\rm erf}(x_r)][{\rm e}^{y_t^2} {\rm erf}(y_t)-{\rm e}^{y_r^2} {\rm erf}(y_r)]}{{\rm CV}_1 {\rm CV}_2}. \nonumber \end{align} \tag{ 74 }$$

Thereby ${\rm erf}(x)$ is the error function [36]. We computed correlation coefficients of the diffusion approximation of the infinitesimal PSP system (meaning a 2D Fokker–Planck equation with the same c and σs) in [21] (results are shown in figure 11(a)).

We used the neural simulation tool NEST [37] to perform numerical simulations of input and output spike trains in the scenario described above. All analyses were based on discretized voltage data obtained during simulations of 1000 s duration using a time resolution of $\Delta t=10^{-4}$ s.

Empirical voltage distributions were obtained by normalizing histograms appropriately. Further smoothing using a simple moving average was performed before comparing these distributions to the analytically obtained stationary distribution. We also performed the comparison using cumulative distributions, as the implicit integration very efficiently reduces the noise in the data. Two 2D distributions are compared via visual inspection of contour lines. We also directly compare spike train cross-correlation functions to assess efficiency and accuracy of the method.

We compute cross-correlation functions of spike trains from conditional PSTHs. One can express this as an integral over two variables $\tau=t_1-t_2$ and $s=t_1+t_2$ with bin size $\Delta$

$$\begin{align} C(\tau)=\frac{1}{\Delta}\int_{\tau}^{\tau+\Delta} \frac{{\rm d}\tau^{\prime}}{u(\tau^{\prime})-l(\tau^{\prime})} \int_{l(\tau^{\prime})}^{u(\tau^{\prime})} \sum_{i,j} \delta(\tau^{\prime}-\tau_i)\delta(s^{\prime}-s_j)\,{\rm d}s^{\prime} \nonumber \end{align} \tag{ 75 }$$

where we set

$$\begin{align} u(\tau)=T/\sqrt{2}+\vert \tau\vert, \qquad l(\tau)=T/\sqrt{2}-\vert \tau\vert. \nonumber \end{align} \tag{ 76 }$$

with observation window T.

The direct singular value decomposition of a 2D membrane potential distribution shows that there are only few singular values that deviate significantly from 0 (figure 4). This behavior does not depend strongly on the chosen discretization, but it does depend on the input correlation coefficient c. Although singular vectors are not probability distributions in their own respect, the singular vectors X_i derived from the neuronal dynamics matrix (except the first few vectors) have the property

$$\begin{align} \label{Sum_to_0} \sum_k X_{ik} \approx 0 \nonumber \end{align} \tag{ 77 }$$

provided the discretization is fine enough. This behavior is demonstrated in figure 7. The contribution of each sum to the overall normalization

$$\begin{align} \sum_{ij}\boldsymbol{p}_{0,ij}=\sum_{i,j} \sum_{k=0,l=0}^{K,L} X_{ik} \Omega_{kl} Y_{jl}=\Sigma_1S\Sigma_2 \nonumber \end{align} \tag{ 78 }$$

is progressively small, where $\Sigma_{1k}$ and $\Sigma_{2l}$ are sums of kth and lth singular vectors of the first and the second matrix for $k\geqslant m$ and $l \geqslant n$, respectively

$$\begin{align} \label{eq_err_sum} \Sigma_{err}(m,n)=\sum_{i,j} \sum_{k=m, l=n}^{K,L} X_{ik} \Omega_{kl} Y_{jl}=\sum_{k=m, l=n}^{K,L} \Sigma_{1k}\Omega_{kl}\Sigma_{2l}. \nonumber \end{align} \tag{ 79 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This shows that the sum converges rather quickly. This error measure is related to projections of 1D discrete stationary distribution $\boldsymbol{p}_0$ (satisfying $\mathcal{M}_1\boldsymbol{p}_0=\boldsymbol{p}_0$) to SVDs. All other eigenvectors of a Markovian matrix ($\mathcal{M}_1\boldsymbol{p}_0=\lambda \boldsymbol{p}_0$ for $\vert \lambda\vert <1$) satisfy $\sum\nolimits_k (P_i)_k=0$. We want to avoid underestimating the total probability mass as a result of the truncated sum in equation (45). Hence, above we justify that the remainder of $\boldsymbol{p}_0$ projections after truncation can be omitted up to a certain precision. On the other hand, in order to describe cumulative contribution of singular vectors we look at the L₁ distance of the omitted remainder (i.e. $k\geqslant m$, $l \geqslant n$)

$$\begin{align} \label{eq_err_sum_2} E(m,n)=\sum_{i,j} \vert \Delta \boldsymbol{p}_{0,ij}(m,n)\vert =\sum_{i,j} \Bigl\vert \sum_{k=m, l=n}^{K,L} X_{ik} S_{kl} Y_{jl} \Bigr\vert \nonumber \end{align} \tag{ 80 }$$

which describes how well the method converges self-consistently. Here we didn't normalize this equation for every term we added. Which means we just rely on fast convergence of $\boldsymbol{p}_0$ projections measured by equation (79), so first few error terms can be misleading.

We started from a simple observation: The stationary joint membrane potential distribution for two neurons with independent input is given as

$$\begin{align} P(V_1,V_2)=P_1(V_1)P_2(V_2) \nonumber \end{align} \tag{ 81 }$$

where $P_1(V_1)$ and $P_2(V_2)$ are the stationary membrane potential distributions of two independent neurons, as described in equations (92) and (93). A similar relation for a discretized voltage grid can be written as

$$\begin{align} \boldsymbol{p}=\boldsymbol{p}_{1} \otimes \boldsymbol{p}_2. \nonumber \end{align} \tag{ 82 }$$

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 small number of additional principal components for any given precision, see figure 4. Here we perform a singular value decomposition of the full solution of equation (50) given in terms of the matrix $\boldsymbol{p}_{ij}$

$$\begin{align} \boldsymbol{p}=\mathcal{L}\Sigma \mathcal{R}^{T} \qquad{\rm{or}}\qquad \boldsymbol{p}_{ij} = \sum_{k=1}^N \sigma_k L_{ik} R_{kj} \nonumber \end{align} \tag{ 83 }$$

generalizing the case of independent neurons to also reconstruct the joint membrane potential distribution for neurons with shared input. As demonstrated in figure 4, convergence is rather quick, even for moderate values of c.

Another aim of our study was to gain some understanding about the influence of the space-time grid. We observed that left and right singular vectors are of the form

$$\begin{align} \boldsymbol{p}_{kl}=\sum_{ij} \Omega_{ij} (X^{c}_{ik}+d_i \delta_{k,r})(Y^{c}_{jl}+d_j \delta_{l,r}) \nonumber \end{align} \tag{ 84 }$$

where X^c and Y^c reflect the quasi-continuous part of basis vectors and $\Omega$ is the coupling matrix as defined in equation (45). Here we need to make sure that the emerging singularity at the reset bin is not causing any numerical problems. One needs to first consider a small time step $\Delta t$ and adapt the stepping in space $\Delta V$ accordingly. A more thorough discussion of a suitable 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 full system. 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 the right singular vectors of single-neuron Markov matrices. As demonstrated in figure 4, right singular vectors lead to an expansion that converges faster for coarse grids (e.g. V  =  0.5 mV). Although for finer grids (e.g. V  =  0.05 mV) the difference is less prominent (figures 5(a) and (b)), right singular vectors still converge slightly faster than left singular vectors (figure 5(d)). Right singular vectors of the single-neuron time evolution matrix yield an orthogonal coordinate system with very good properties.

Zoom In Zoom Out Reset image size

As reported previously [21], we may also use a direct analytical approach using the basis of the single-neuron Fokker–Planck operator, and its adjoint basis

$$\begin{align} \boldsymbol{p}_0(X,Y)=\sum_{ij}\Omega_{ij} X_i \otimes Y_j \nonumber \end{align} \tag{ 85 }$$

where X and Y are the left eigenvectors of the single neuron operators. However, the issue is that the discrete adjoint basis blows up at the lower boundary. The effect of this on our approximation is demonstrated in figure 6. In general, SVD eliminates a kernel of singular matrices.

Zoom In Zoom Out Reset image size

In our treatment of the 2D Fokker–Planck equation in [21], which is the infinitesimal limit of the theory considered above, we used the basis and adjoint basis to project linear operators to a subspace. This has certain advantages as it satisfies constraints for marginal distributions and preserves the Markov property to some extent. Positivity of the solution in the subspace is not guaranteed, but time evolution is probability preserving ($\sum\nolimits_i P_i(t)=1$).

First, SVD is computationally convenient, because it leads to using a real orthogonal basis which resembles the eigenbasis of $\mathcal{M}$. Singular vectors (e.g. X_i) quickly converge to the span of the non-orthogonal basis $\mathcal{M}$, as equation (77) holds for sufficiently big i (shown at figure 7(d)). Second, numerical instabilities due to an ill-conditioned time evolution matrix $\mathcal{M}$ (some eigenvalues $\lambda_i \approx 0$) are cured by SVD. Third, although the basis vectors implied by SVD have the disadvantage of not completely preserving positivity, the deviation remains within tight bounds even for a relatively small number of basis vectors.

Zoom In Zoom Out Reset image size

We compare our SVD-based Markov theory and direct numerical simulations of spiking neurons both on the level of joint 2D membrane potential distributions and on the level of spike train covariance functions, see figure 8. The empirical distributions derived from 2D histograms are slightly smoothed in order to compare them to the distributions derived from the Markov theory on the level of contour lines. We also considered 2D cumulative distribution functions, where the smoothing step can be omitted. Moreover, we computed output spike train covariance functions as described in methods section and compared them to the covariance functions obtained directly from the simulated spike trains.

Zoom In Zoom Out Reset image size

Two neurons that are driven by correlated input current will exhibit correlated output spike trains. This transfer of correlation reflects an important property of neuronal dynamics, which is of particular relevance for understanding the contribution of neurons to network dynamics. Recently, we were able to demonstrate, by exact analytical treatment, that the correlation transfer for leaky integrate-and-fire neurons is strongly non-linear [21]. Only for weak input correlation it can be described by perturbative methods, and deviations from linear response theory depend on the background firing rate. In the present work we demonstrate the same non-linear correlation transfer, see figure 11. Low output firing rates generally require a non-perturbative treatment, while the approximation derived from linear response theory [38] is reasonably precise if firing rates are high (as we demonstrated in [21] figure 5(a)). We considered firing rates between 1 and 25 Hz, and values for ${\rm CV}^2$ between 0.5 and 1, consistent with what is reported in neocortical neurons in vivo.

There we also demonstrate how the parameters of the spatial and temporal coarse graining affects the precision of the Markov approximation. Our main conclusion is that dimensional reduction via SVD subspace projections makes it possible to achieve a superior precision with small bin sizes. Fine enough grids, however, could not be dealt with on typical computers without using the reduction suggested here. Hence our approach in comparison to [24] has practical advantage not only in partially solving the technical problem of singularity for e.g. c  >  0.7, but also in providing computational economy. This is practical regarding the capacity of an everyday computer which can hardly deal with arrays of order $10^4 \times 10^4$ with floating point precision.

The study presented here is complementary to our work on LIF neurons with shared white noise [21]. Although we have some parameter restrictions here as our approach is numerical and parameters affect the computational cost, we can safely extrapolate all results in [21], to our current work with finite jump amplitudes and Poisson statistics. Our approach can deal with almost the full range of input correlations $-1\leqslant c<1$, and the expansion converges fast (figures 5(b) and (d)). Also, our method is widely generalizable [24].

In this work we focus on strongly correlated neuronal dynamics. We calculated several examples of joint membrane potentials and correlation functions in strongly correlated regime (presented in figures 8–11). All in all, the regime $c\rightarrow 1$ is indeed relevant. Strong correlations of membrane potentials were observed in nearby neurons of cortical networks [2–4], compatible with the high degree of shared input suggested from neuroanatomical studies. On the contrary, there are studies that conclude that neurons decorrelate [9] by an active mechanism [10, 12]. Here we show that decorrelation in the strongly correlated regime can also arise from parameter heterogeneities (figure 10(a)). Such parameters can be external variables such as input current mean and variance, or intrinsic parameters such as membrane time constant, size of PSPs and threshold.

Neurons in biological networks have widely distributed parameters, and this heterogeneity may also influence information processing [41–43]. Moreover, robust asymmetries in spike correlations could lead to asymmetric synaptic efficacies, if they are subject to spike timing dependent plasticity [44, 45]. Our approach reveals a generic temporal asymmetry in cross-covariance functions, related to the heterogeneity of intrinsic neuron parameters and input variables. Such temporal asymmetry is more pronounced for larger values of c, especially in the non-perturbative regime that we address in this work.

We document here an application of our method to four types of asymmetries [21, 43]: μ asymmetry, σ asymmetry, $\tau_m$ asymmetry and $V_{{\rm th}}$ asymmetry. We quantified the asymmetry by $a=\chi_1/\chi_2$ (specific values given in table A2), where χ is replaced by the respective parameter. Asymmetric correlations have been described previously, and they were by numerical simulations and experimentally studied by [41, 43, 46, 47].

In the strongly correlated regime the correlation transfer function is non-linear [31, 38] and the dynamics is quite sensitive to heterogeneities of the input and of the parameters assigned to the two neurons [41–43]. Recent experiments demonstrated that asymmetric correlation functions arise in neocortical neurons as well [41–43]. Correlation asymmetries could make an important contribution to structure formation in networks through Hebbian learning on short time scales in the range of the membrane time constant of neurons [44, 45]. All in all it is important to demonstrate that our method predicts joint membrane potential distributions (figure 10), correlation functions (figure 8 and (figure 9)) and correlation coefficients (figure 11) to a good extent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed that our results in figure 9 shows a slight mismatch with the empirical correlation functions around $\tau=0$. There we present only the time-resolved part of $C(\tau)$, ignoring the delta-function $r_0\delta(\tau)$. It is evident from examples in figure 11 that the integral of the correlation function around $\tau=0$ does not depend strongly on the grid, as long as a fine enough grid is selected. However, depending on $\Delta t$ and $\Delta V$, while amplitude r₀ can be over-estimated, $r_{1\vert 2}(\epsilon)$ can be underestimated for a small $\vert \epsilon\vert >0$, or vice verse. The main factor that induces such a mismatch around $\tau=0$ is $\Delta t$, as multiple jumps over the threshold $(V_{{\rm th}}, V_{{\rm th}})$ to $(V_1>V_{{\rm th}}, V_1>V_{{\rm th}})$ simply determine r₀, an over-estimation there would lead to an under-estimation of $r_{1\vert 2}(\epsilon)$. On the other hand, we have selected (the numerical limit $\tau \rightarrow 0$) $\epsilon=\Delta t/2$ (e.g. in figure 9) in order to avoid synchronous spikes in the empirical distribution. However, a larger bin $\epsilon=\Delta t$ would indeed underestimate the limiting value (results not shown here).

The precision of the approximation depends on several other factors, one is the number of SVDs. When $c\rightarrow 1$ for symmetric parameters (in practice c  >  0.95 is hard to deal with), it is evident that the membrane potential becomes singular and the expansion requires more and more terms. A second factor is the selected precision of jump distribution. This translates to number of terms we sum in equation (55). Another factor that may change the approximation is given by the asymptotic rates. In [25] the probability density obtained by the three different methods: Direct simulation, numerical solution of the Markov system, and Fokker–Planck theory. The agreement of all three variables in the relevant range below the threshold is excellent. However, the equilibrium firing rate shows (for 2 neurons, denoted here by r₁ and r₂) a deviation in the range of 0.5 Hz. Although in that case empirical and theoretical correlation functions could have a similar shape, the difference in asymptotic rates (e.g. $\lim\nolimits_{\tau\rightarrow \infty} r_{1\vert 2}(\tau)=r_1$ and r₂ can lead to a shift in the central peak e.g. $C_{12}(\epsilon)=r_{1\vert 2}(\epsilon)r_2$.

We point out that overall there are four parameters in jump equations (92) and (92), which can result in asymmetric 2-dimensional joint distributions: h_a, $\tau_{m, a}$, $\varphi_a=\frac{\mu_a}{h_a\tau_{m, a}}$ and $V_{{\rm th}, a}$ (as we can fix V_r without loss of generality). Such asymmetries were studied due to parameter heterogeneity of $\mu_a$, $\sigma_a$, $\tau_{m, a}$ and $V_{{\rm th}, a}$ in figures 9 and 10. In our parametrization $\sigma_a$ heterogeneity is a result of h_a heterogeneity when $\tau_m$ is constant. However, the asymmetry in the diffusion approximation can be reduced to dimensionless variables $x_{{\rm th}, a}=\frac{V_{{\rm th}, a}-\mu_a\tau_{m, a}}{\sigma_a\sqrt{\tau_{m, a}}}$ and $x_{r, a}=\frac{V_{r, a}-\mu_a\tau_{m, a}}{\sigma_a\sqrt{\tau_{m, a}}}$, and the membrane time constant $\tau_{m, a}$ [21].

The Poisson model with large input rates r_a and small PSPs h (relative to $V_{{\rm th}}$) fits well to diffusion approximation [25, 32]. As long as Poisson input doesn't deviate significantly from gaussianity the forth parameter is redundant. For example proportional increase in h can compansate increase in $V_{{\rm th}}$ such that $N_h=\frac{V_{{\rm th}}}{h}$ is constant. Nevertheless, our treatment is general e.g. we can deal with input models that lead to non-Gaussian jump distributions [31, 33]. Such non-Gaussianity can be result of input with small firing rates (in our study, practically small with respect to $1/\Delta t$), large PSPs [25] or higher order correlations [31, 33].

Models of correlated neuronal activity describe the origin of correlations in spiking model neurons, induced by the structure of the network and/or feedforward input. Such neuron models, however, are notoriously nonlinear. Nevertheless, most treatments rely on linearization and other simplifying assumptions, as nonlinear correlation transfer functions (i.e. relation of input and output correlations) are difficult to derive. Previous analytical approaches have employed perturbation theory [16, 17] to study pairwise correlations under the assumption of weak input correlation [35, 38]. As a consequence of this technical convenience, we still 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.

All computations described in our paper can be applied to more general integrate-and-fire models with a non-linear leak function $\Psi(V)$

$$\begin{align} \tau_m\dot{V}=\Psi(V)+h_{{\rm ex}}\tau_m S^{{\rm ex}}(t) +h_{{\rm in}} \tau_m S^{{\rm in}}(t). \nonumber \end{align} \tag{ 86 }$$

We only need to rewrite the decay matrix as a discrete approximation of the differential operator $D(x) = \tau_m \frac{\rm d}{{\rm d}t} x - \Psi(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 could also 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 amplitude distributions [33, 48]. The method described in the present paper can be adapted to such scenarios by simply using a modified definition of the jump matrix [31]. This requires generalization of jumps matrices ${\rm e}^{J_1}\rightarrow \mathcal{J}_1$ and jump distribution of shared spikes $c(\,j) \rightarrow\psi(\,j)$.

Higher-order statistics on the output side is also compatible with our method, describing the 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 same fashion

$$\begin{align} \boldsymbol{p}_{lmn}=\sum_{ijk}\Omega_{ijk} X_{il} Y_{jm} Z_{kn} \nonumber \end{align} \tag{ 87 }$$

now with a 3D jump operator given in the generic form

$$\begin{align} \mathcal{J}_{\rm 3D}=\mathcal{J}_1 \otimes \mathcal{J}_2 \otimes \mathcal{J}_3. \nonumber \end{align} \tag{ 88 }$$

This operator is again transformed with a basis derived from a SVD as

$$\begin{align} \mathcal{J}_1 \otimes \mathcal{J}_2 \otimes \mathcal{J}_3 \rightarrow X \mathcal{J}_x X^{T} \otimes Y \mathcal{J}_2 Y^{T} \otimes Z \mathcal{J}_3 Z^{T}. \nonumber \end{align} \tag{ 89 }$$

This procedure is computationally more demanding as we need to consider an additional dimension, although the scaling is not exponential. Under assumption of homogeneous shared input (same jump amplitudes driven by shared input in all directions) leads to an expression similar to equation (55),

$$\begin{align} \label{eq_J_3D} {\mathcal{J}}_{\rm 3D}=\sum_{j \: \in \: {\mathcal{J}}} \psi(\,j)[{\rm e}^{\mathcal{J}_1}O^{\,j}\otimes {\rm e}^{\mathcal{J}_2}O^{\,j} \otimes {\rm e}^{\mathcal{J}_3}O^{\,j}]. \nonumber \end{align} \tag{ 90 }$$

Assuming joint stationarity of all three spike trains S₁(t), S₂(t) and S₃(t), we need to find the joint third moment of the spike train statistics

$$\begin{align} \label{eq_third_moment} \mu_{123} (\tau_1,\tau_2)= \bigl\langle S_1(\tau_1) S_2(\tau_2) S_3(0) \bigr\rangle. \nonumber \end{align} \tag{ 91 }$$

In general a simplified approach to spiking neural network dynamics is the ultimate goal of theoretical studies [11, 13–18, 35, 38]. In case of a large scale interacting network of spiking neurons a brut-force approach is doomed to fail. We argue that our approach can be generalized to interacting neurons as well, by modifying the time-evolution tensor $\mathcal{M}$ and its counterpart $\mathcal{Q}$. This can be achieved by an additional multiplicative term $\mathcal{K}_{\rm 2D}(\boldsymbol{p})$ representing the effect of presynaptic spike trains. Yet, it is neither trivial to construct, nor easy to treat such a non-linear interaction term. Merely immense dimensionality of the system limits possible computations.

As far as the theme of this paper goes, we believe a weak input case (or a.k.a perturbative regime) isn't of interest to us. Weak input, leading to a linear modification, can be integrated to the approach implemented above or by some well-known treatments [16, 17]. The interesting point here is understanding dynamics of spiking neural networks beyond perturbative regime e.g. stochastic dynamics on the verge of synchronization. For instance, computations of correlations is of extreme importance in a strongly coupled network regime. Yet a general method thereof is not introduced so far. We believe a modification of our method could provide some further insight into balanced networks with strong input [18].

All in all, we want to demonstrate here some technical issues of such an approach. As an example, we can look at a case with excitatory interactions in non-perturbative regime. An oversimplified scenario is that the recurrent input can shift the distribution like any other Poisson input with shared and private components as in figure 1. This simple model is

$$\begin{align} \dot{V_1}(t)=-\frac{V_1(t)}{\tau_{m,1}}+h_{{\rm ex}}(S^{{\rm ex}}_{1}(t)+S_{c}^{{\rm ex}}(t))+ h_{{\rm in}}(S^{{\rm in}}_{1}(t) + S_{c}^{{\rm in}}(t))+ j_{{\rm ex}} S_2(t) \label{eq_model_1} \nonumber \end{align} \tag{ 92 }$$

$$\begin{align} \dot{V_2(t)}=-\frac{V_2(t)}{\tau_{m,2}}+h_{{\rm ex}}(S^{{\rm ex}}_{2}(t)+S_{c}^{{\rm ex}}(t))+ h_{{\rm in}}(S^{{\rm in}}_{2}(t) + S_{c}^{{\rm in}}(t)) \label{eq_model_2} \nonumber \end{align} \tag{ 93 }$$

where we consider only connections ($2\rightarrow 1$) and S₂ represents the output spike trains of the neuron 2. Here we simply assume that the output spike trains are Poisson, which means we will omit the time resolved correlation function as in figure 9 or 5, just for sake of argument. Assume that κ is the discrete jump index driven by excitatory connections with PSP $j_{{\rm ex}}$ and U is an up-transition matrix. The instantaneous counts of the incoming spikes relative to V₁ are given by $\mathcal{F}_{1\vert 2}(\boldsymbol{p})$, where $\mathcal{F}_{\cdot\vert \cdot}(\cdot)$ is simply the flux operator equation (66). In general the time evolution operator has the form

$$\begin{align} \mathcal{M}_{\rm 2D}(\boldsymbol{p})=\mathcal{T}_{\rm 2D} \mathcal{K}_{\rm 2D}(\boldsymbol{p})\mathcal{J}_{\rm 2D} \mathcal{D}_{\rm 2D}. \nonumber \end{align} \tag{ 94 }$$

In which case $\mathcal{K}_{\rm 2D}(\boldsymbol{p})$ represents interactions and the simple example above reads

$$\begin{align} \mathcal{K}_{\rm 2D}(\boldsymbol{p})=\Big{[}\Big{(}I-{\rm diag}(\mathcal{F}_{1\vert 2}(\boldsymbol{p})) +U^{\kappa}{\rm diag}(\mathcal{F}_{1\vert 2}(\boldsymbol{p}))\Big{)}\otimes I]\Big{]}. \nonumber \end{align} \tag{ 95 }$$

where ${\rm diag}(\boldsymbol{x})$ is a diagonal matrix with components $x_i\delta_{ij}$. Note that $\mathcal{K}_{\rm 2D}(\boldsymbol{p})$ can have a general form with multiple synaptic interactions and multiple transitions in finite time $\Delta t$ etc. The non-linear equation with stationarity condition (assuming that the system is Markovian) then reads

$$\begin{align} \mathcal{M}_{\rm 2D}(\boldsymbol{p})\boldsymbol{p}=\boldsymbol{p} \nonumber \end{align} \tag{ 96 }$$

notice that if $\kappa=0$, this reduces to uncoupled equation we studied in this work. Dealing with this system requires a further study, as we can't easily handle it with our current techniques.

Here an ad hoc trick can be implemented e.g. first studying Gaussian fixed point of the network dynamics via self consistent equations [16] and using the background parameters in order to set-up SVDs and so on. Second, finding an initial background voltage distribution $\boldsymbol{p}_{(0)}$ using

$$\begin{align} \mathcal{M}_{\rm 2D}\boldsymbol{p}_{(0)}=\boldsymbol{p}_{(0)} \nonumber \end{align} \tag{ 97 }$$

Third, implementing an iteration starting with $K_{\rm 2D}(\boldsymbol{p}_{(0)})$ and improving the distribution via

$$\begin{align} \mathcal{M}_{\rm 2D}(\boldsymbol{p}_{(0)})\boldsymbol{p}_{(1)}=\boldsymbol{p}_{(1)} \nonumber \end{align} \tag{ 98 }$$

and so on for $\boldsymbol{p}_{(n)}\rightarrow \boldsymbol{p}_{(n+1)}$. This is advantageous because we can use a fixed SVD basis. Yet, it is rather a computationally vast procedure with no guaranteed success. No doubt, this case suggests further investigation, especially regarding the output spike train correlation. Such non-Markovian input may be of key importance to understand possible deviations from Gaussianity at the boundary of a transition from weakly coupled regime to strongly coupled regime [14, 18], which shows non-trivial temporal properties of correlation functions.

The joint membrane potential distribution has a singularity at the origin $(V_{1}, V_{2})=(0, 0)$ due to a coordinated reset caused by some shared input spikes. There is also a line discontinuity at V₁  =  0 and V₂  =  0, again due to the reset boundary condition. These singularities are reflected in the right singular vectors X and Y. This is the exact reason why we selected them as a basis to expand operators and joint distributions.

We observed that a singularity (a δ-function) emerges when $\Delta V$ is small and $\Delta t$ is large, in relative terms. This is an issue even for the 1D discrete problem [25], and it is even more severe for 2D problems. The amplitude of the singularity is constant and hence $\boldsymbol{p}_{00}$ scales with $(\Delta V){}^{-2}$. This phenomenon occurs only if PSPs have a finite amplitude. As the PSP gets larger relative to $\Delta V$, reset currents remain finite even in continuous time [25]. As a consequence, the limit to continuous variables must be taken with care, in particular for c  >  0.

The δ-singularity does not exist for the diffusion approximation [21]. However, the definition of the current at the origin again fails as the derivative is discontinuous in both V₁ and V₂ directions. There is no doubt that the jump equation is well-defined [49] as the flux at the boundary is not local. However, the infinitesimal limit is problematic for correlated neurons (c  >  0) since the flux is not defined at the boundary of the 2D domain $(V_{{\rm th}}, V_{{\rm th}})$.

To check the accuracy of our theory, we compare it to a numerically obtained solution for the equilibrium distribution, correlation function and firing rate. In 1D and 2D alike, we model the membrane potential as a Markov process the state of which is completely determined by the membrane potential distribution $\boldsymbol{p}(v_1, v_2, t)$. The discrete time process evolves between adjacent time steps $\Delta t$. As stated in [25] where authors showed a 1D discrete model, these transitions are composed of the sequence of decay ($\mathcal{D}$), voltage jumps due to synaptic input ($\mathcal{J}$) and the thresholding operation ($\mathcal{T}$). Likewise we construct 2D operators in a similar fashion and decompose the correlated jump distribution by using SVDs. We discretize the voltage space and determine the transition matrix $\mathcal{(TJD)}_{\rm 2D}$ between discrete voltage bins. Due to the Perron–Frobenius theorem [50], the column-stochastic transition matrix (satisfying $\sum\nolimits_j\mathcal{(TJD)}_{ij}=1$ $\forall i$) has an eigenvalue $\lambda=1$ that corresponds to the equilibrium state $\boldsymbol{p}_0(v_1, v_2)$. We determine the corresponding eigenvector numerically which gives the equilibrium distribution and allows us to calculate the correlation function. For details see section 'models and methods'.

The selection of an appropriate grid in space and time is crucial for correlation computations. The small residual offset between direct simulation and our new semi-analytical computation (see figures 9 and 10), for example, can be considered as a discretization artifact. Although this issue would deserve a more systematic treatment, we report here some observations that can guide grid selection:

• (i)  
  For discrete solutions of the heat equation based on central difference scheme, convergence of 1D time evolution requires $\frac{\Delta t}{(\Delta V){}^2} < \frac{\tau}{\sigma^2}$ [51]. A similar rule also applies in the 2D case considered here. In general, explicit discretization schemes of second order differential operators arising in the study of diffusion, require positivity and stability conditions in the order of $\Delta t = O(\frac{\tau (\Delta V_1+\Delta V_2){}^2}{\sigma^2})$ [24]. In this work, we followed a discretization scheme that approximately conserves probability, a Markovian approximation [25]. However, we note once more that some grids may lead to violation of the Markov property for too large $\Delta t$, as a result of boundary effects. This may create issues when the largest eigenvalue exceeds 1.
• (ii)  
  To reflect small expected bin counts (especially for $c \approx 1$) adequately, one needs a larger $\Delta t$ and a smaller $\Delta V$. This is in conflict with rule (i). Besides, we observe that smaller $\Delta V$ for a fixed $\Delta t$ actually leads to better firing rate approximation up to some point (figure 12). The variance of the input ($\sigma^2$) is crucial in determining the appropriate temporal bin size, while the mean input (μ) is less effective. With increasing variance one observes an increasing firing rate error (figure 12).
• (iii)  
  A finer grid requires more computational efforts to achieve a smooth correlation function. The SVD reduction does not alter this behavior. Other dependencies and limiting factors are indicated in figure 11. There are two constraining factors which are determined by the selected precision of the approximation. One is the extent of the jump distribution, which affects the number of terms to be accumulated (size of the set $\mathbb{J}$ in equation (55)). For a fixed grid and a selected precision, this number increases with $\sigma_c$. The second constraint is the size of the SVD subspace. We know that as c gets closer to 1 and $C(\tau)$ gets steeper we need to include more singular components.

Zoom In Zoom Out Reset image size

In figure 11 we illustrate how coarse graining affects the shape of the cross-covariance function $C(\tau)$. Although the precision of the approximation is limited by the subspace projections implied by SVD, the grid parameters $\Delta t$ and $\Delta V$ are the most important factors to get the shape of the function right. However, for a fixed dimension of the SVD subspace even the finest grid would not be able to capture the singularity at zero time lag ($\tau=0$). The grid effectively limits the precision of the approximation due to the reduced number of degrees of freedom.