A kinetic theory approach to capturing interneuronal correlation: the feed-forward case

Abstract

We present an approach for using kinetic theory to capture first and second order statistics of neuronal activity. We coarse grain neuronal networks into populations of neurons and calculate the population average firing rate and output cross-correlation in response to time varying correlated input. We derive coupling equations for the populations based on first and second order statistics of the network connectivity. This coupling scheme is based on the hypothesis that second order statistics of the network connectivity are sufficient to determine second order statistics of neuronal activity. We implement a kinetic theory representation of a simple feed-forward network and demonstrate that the kinetic theory model captures key aspects of the emergence and propagation of correlations in the network, as long as the correlations do not become too strong. By analyzing the correlated activity of feed-forward networks with a variety of connectivity patterns, we provide evidence supporting our hypothesis of the sufficiency of second order connectivity statistics.

Appendices

Appendices

A Method for solving the kinetic theory equations

The first step in solving the kinetic theory equations (2) is to rewrite them in conservative form, i.e., divergence form $\frac{\partial {\rho}}{\partial {t}} = - \nabla \cdot {\mathbf{J}}$, where J is a flux density. The terms describing the advection due to the leak current are already the divergence of a flux density. We need to rewrite the integrals corresponding to the voltage jumps in response to excitatory input. Recall that F _A (x) is the complementary cumulative distribution function and f _A (x) is the probability density of the random jump size A, so that $\frac{\partial F_A}{\partial x}=-f_A(x)$. Applying this identity and the fundamental theorem of calculus, we obtain the following equalities. We rewrite the terms due to independent input as

$$\begin{array}{lll} &\int_{v_{\mathrm{reset}}}^{v_1}f_A(v_1-\theta_1)\rho(\theta_1,v_2,t){\text d}\theta_1-\rho(v_1,v_2,t) \\\\ &{\kern12pt} =-\frac{\partial}{\partial v_1}\int_{v_{\mathrm{reset}}}^{v_1}F_A(v_1-\theta_1)\rho(\theta_1,v_2,t){\text d}\theta_1, \notag \\\\ &\int_{v_{\mathrm{reset}}}^{v_2}f_A(v_2-\theta_1)\rho(v_1,\theta_2,t){\text d}\theta_2-\rho(v_1,v_2,t) \\\\ &{\kern12pt} =-\frac{\partial}{\partial v_2}\int_{v_{\mathrm{reset}}}^{v_2}F_A(v_2-\theta_2)\rho(v_1,\theta_2,t){\text d}\theta_2, \end{array} $$

(15)

The integrals involving F _A are probability flux densities analogous to those defined in Nykamp and Tranchina (2000).

The jumps due to synchronous input are a little more complicated. There are many ways to write them in conservative form. Since the evolution of ρ is symmetric in v ₁ and v ₂, we kept the equation symmetric by dividing the diagonal jumps into a symmetric pair of horizontal and vertical jumps.

$$\begin{array}{lll} &&{\kern-7pt}\int_{v_{\mathrm{reset}}}^{v_1}\int_{v_{\mathrm{reset}}}^{v_2}f_A(v_1-\theta_1)f_A(v_2-\theta_2)\rho(\theta_1,\theta_2,t){\text d}\theta_1{\text d}\theta_2 \\\\ &&{\kern21pt} -\rho(v_1,v_2,t) \\\\ &&{\kern19pt}=-\frac{1}{2}\frac{\partial}{\partial v_1}\biggl(\int_{v_{\mathrm{reset}}}^{v_1}F_A(v_1-\theta_1)\rho(\theta_1,v_2,t){\text d}\theta_1 \\\\ &&{\kern70pt} +\int_{v_{\mathrm{reset}}}^{v_1}\int_{v_{\mathrm{reset}}}^{v_2}F_A(v_1-\theta_1)f_A(v_2-\theta_2) \\\\ &&{\kern70pt} \times\rho(\theta_1,\theta_2,t){\text d}\theta_2{\text d}\theta_1\biggr) \\\\ &&{\kern21pt} -\frac{1}{2}\frac{\partial}{\partial v_2}\biggl(\int_{v_{\mathrm{reset}}}^{v_2}F_A(v_2-\theta_2)\rho(v_1,\theta_2,t){\text d}\theta_2 \\\\ &&{\kern61pt} +\int_{v_{\mathrm{reset}}}^{v_2}\int_{v_{\mathrm{reset}}}^{v_1}F_A(v_2-\theta_2)f_A(v_1-\theta_1) \\\\ &&{\kern62pt} \times\rho(\theta_1,\theta_2,t){\text d}\theta_1{\text d}\theta_2\biggr). \end{array} $$

(16)

Although these integrals do not correspond to physical horizontal and vertical jumping of the voltage, we view their formulation simply as an intermediate step to developing a numerical scheme to solve the original Eq. (2).

As a result, we obtain the following integro-differential equation:

$$\begin{array}{lll} \frac{\partial \rho}{\partial t}(v_1,v_2,t)&\!=\!&-\!\nabla\cdot{\mathbf{J}}(v_1,v_2,t) \!+\!\delta(v_1\!-\!v_{\mathrm{reset}})J_{\mathrm{reset},1}(v_2,t) \\\\ && +\ \delta(v_2\!-\!v_{\mathrm{reset}})J_{\mathrm{reset},2}(v_1,t) \\\\ && +\ \delta(v_1-v_{\mathrm{reset}})\delta(v_2-v_{\mathrm{reset}})J_{\mathrm{reset},3}(t), \label{eq:kt_conservative} \end{array} $$

(17)

where

$$\begin{array}{rll} {\mathbf{J}}(v_1,\!v_2,\!t)&\!=\!&{\mathbf{J}}_{\mathrm{leak}}({\kern-.5pt}v_1,\!v_2,\!t{\kern-.5pt})\!+\!{\mathbf{J}}_{\mathrm{ind}}(v_1,v_2,t)\!+\!{\mathbf{J}}_{\mathrm{syn}}({\kern-.5pt}v_1,\!v_2,\!t{\kern-.5pt}), \\\\ {\mathbf{J}}_{\mathrm{leak}}&\!=\!&(J_{\mathrm{leak}}^1,J_{\mathrm{leak}}^2),{\mathbf{J}}_{\mathrm{ind}}\!=\!(J_{\mathrm{ind}}^1,J_{\mathrm{ind}}^2), \\\\ {\mathbf{J}}_{\mathrm{syn}}&\!=\!&(J_{\mathrm{syn}}^1,J_{\mathrm{syn}}^2), \\\\ J_{\mathrm{leak}}^1(v_1,v_2,t) &\!=\!&-\frac{v_1-Er}{\tau}\rho(v_1,v_2,t), \\\\ J_{\mathrm{ind}}^1(v_1,\!v_2,\!t)&\!=\!&\nu_{\mathrm{ind}}(t)\int_{v_\mathrm{reset}}^{v_1}F_A(v_1-\theta_1)\rho(\theta_1,v_2,t){\text d}\theta_1, \\\\ J_{\mathrm{syn}}^1(v_1,\!v_2,\!t)&\!=\!&\frac{1}{2}\nu_{\mathrm{syn}}(t)\biggl(\!\int_{v_{\mathrm{reset}}}^{v_1}\!F_A(v_1-\theta_1)\rho(\theta_1,v_2,t){\text d}\theta_1 \\\\ &&{\kern35pt}+\!\int_{v_{\mathrm{reset}}}^{v_1}\!\int_{v_{\mathrm{reset}}}^{v_2}\!F_A(v_1\!-\!\theta_1)\!f_A(v_2\!-\!\theta_2) \\\\ &&{\kern35pt}\times\!\rho(\theta_1,\theta_2,t){\text d}\theta_2{\text d}\theta_1\!\biggr), \end{array} $$

(18)

All second components are analogous by symmetry. The reset terms are defined in Eq. (3).

We do not have local conservative equality between the flux density across the threshold $(J_{\mathrm{ind}}^1(v_{\mathrm{th}},v_2,t)+J_{\mathrm{syn}}^1(v_{\mathrm{th}},v_2,t))$ and the reset J _reset,1(v ₂,t). Due to our definition of the flux density $J_{\mathrm{syn}}^1(v_{\mathrm{th}},v_2,t)$, it includes the effect of synchronous inputs where the voltage of neuron 2 jumps from v ₂ to higher voltages (and hence appears in the reset term J _reset,1(v ₂,t) at those higher voltage or even in J _reset,3(t) if neuron 2 also crossed threshold). Nonetheless, one can verify that the system is globally conservative because the total flux across threshold equals the total reset:

$$\begin{array}{lll} &\int_{v_{\mathrm{reset}}}^{v_{\mathrm{th}}}(J_{\mathrm{ind}}^1(v_{\mathrm{th}},v_2,t)+J_{\mathrm{syn}}^1(v_{\mathrm{th}},v_2,t)){\text d}v_2 \\\\ &{\kern21pt} +\int_{v_{\mathrm{reset}}}^{v_{\mathrm{th}}}(J_{\mathrm{ind}}^2(v_1,v_{\mathrm{th}},t)+J_{\mathrm{syn}}^2(v_1,v_{\mathrm{th}},t)){\text d}v_1 \\\\ &{\kern21pt} =\int_{v_{\mathrm{reset}}}^{v_{\mathrm{th}}}J_{\mathrm{reset},1}(v_2,t){\text d}v_2+\int_{v_{\mathrm{reset}}}^{v_{\mathrm{th}}}J_{\mathrm{reset},2}(v_1,t){\text d}v_1 \\\\ &{\kern21pt} +J_{\mathrm{reset},3}(t) \end{array} $$

(19)

To numerically solve the equations, we used the finite volume approach. We discretized the domain [0,v _th] ×[0,v _th] in (v ₁,v ₂) space into squares of size Δv ×Δv with Δv = 0.0125. We let v _1,j  = v _2,j  = (j − 1/2)Δv so that the point (v _1,j , v _2,k ) is the center of square (j,k). Each point ρ _j,k(t) represented the integral of ρ(v ₁,v ₂,t) over the square (j,k). The divergence $\nabla \cdot {\mathbf{J}}$ at the center of each square was approximated by differences in the integral of the flux along the boundary of the square, and we approximated the resulting integrals along each line segment with the midpoint rule. We linearly interpolated ρ _j,k to estimate its value along the boundary for J _leak. To estimate the value of the integrals for J _ind and J _syn, we used Simpson’s rule.

We reduced the dimension of the system of equations by exploiting the symmetry ρ _j,k = ρ _k,j. We discretized in time using the trapezoid method with Δt = 0.5 ms. At each time step, we solved the system of linear equations using GMRES (Saad and Schultz 1986). In this way, we used a numerical method that is second order accurate in both time and space.

Since we set $v_{\mathrm{reset}} \!<\! E_\text r$, the advection flux J _leak is non-zero along locations of voltage reset, preventing the J _reset,1 and J _reset,2 delta-function source terms from forming a delta-function in ρ(v ₁,v ₂,t). However, the double delta-function source due to the reset J _reset,3 will form a delta-function component of ρ. The flux J _reset,3 represents a finite probability per unit time that both neurons reset to (v _reset,v _reset). At that point, the voltage pair will advect along the line $\{(v_1,v_2)|v_1\!=\!v_2, 0 \!<\! v_1 \!<\! E_\text r\}$ until either neuron receives an input. Hence, there will be finite probability that the voltages of a neuron pair reside along this line. The delta-function in ρ along this line will prevent second order convergence of the above numerical method if we apply it directly to Eq. (17). To ensure convergence, we divide ρ into two components, a smooth component ρ _s and a delta function component along the diagonal with weight function ρ _δ :

$$ \rho(v_1,v_2,t)=\rho_s(v_1,v_2,t)+\delta(v_1-v_2)\rho_{\delta}(v_1) $$

(20)

Plugging this expression into the original evolution Eq. (17) and then grouping the delta-function terms will result in coupled evolution equations for ρ _s and ρ _δ . Since ρ _s and ρ _δ are smooth, we can solve the resulting equations using our numerical scheme and achieve second order accuracy.

B Calculating connectivity statistics for different network classes

We calculate the expected number of inputs W ¹ and the fraction of shared inputs β for networks with arbitrary outgoing degree distributions and incoming connections determined randomly. We examine the connectivity from a presynaptic population 1 with N ₁ neurons onto a postsynaptic population 2 with N ₂ neurons.

Let $\hat{W}_{ij}=1$ if there is a connection from neuron j in population 1 onto neuron i in population 2; otherwise, let $\hat{W}_{ij}=0$. Let $d_j = \sum_{i=1}^{N_2} \hat{W}_{ij}$ be outgoing degree of neuron j in population 1, i.e., the number of connections from neuron j onto all neurons in population 2. Let the function f(k) be the outgoing degree distribution so that

$$ \Pr(d_j = k) = f(k) $$

for k = 1,2, ..., N ₂. The expected total number of connections (out of N ₁ N ₂ possible) is $N_1 \sum_{k=1}^{N_2} k f(k)$ so that the expected number of connections onto any neuron in population 2 is

$$ \label{eq:Wfromf} W^1= \frac{N_1}{N_2} \sum\limits_{k=1}^{N_2} k f(k). $$

(21)

If f(k) was given by a one-parameter family of distributions, prescribing W ¹ would determine the particular distribution as a function of population sizes N ₁ and N ₂.

Neuron j in population 1 has d _j connections onto neurons in population 2, which we assume are assigned randomly to neurons in population 2. Then, conditioned on this value of d _j , the probability of a connection from neuron j onto any given pair (indexed by i ₁ and i ₂) of neurons in population 2 is

$$ \Pr(\hat{W}_{i_1j}=1 \,\&\, \hat{W}_{i_2j}=1 \,|\, d_j) = \frac{\binom{N_2-2}{d_j-2}}{\binom{N_2}{d_j}} = \frac{d_j(d_j-1)}{N_2(N_2-1)} $$

Multiplying by the probability distribution of d _j (i.e., the outgoing degree distribution) and summing over all possible values of d _j , we determine that the probability a given neuron j in population 1 is connected to a given pair of neurons in population 2 is

$$\begin{array}{lll} \sum\limits_{k=1}^{N_2} \Pr(\hat{W}_{i_1j}&=&1 \,\&\, \hat{W}_{i_2j}=1 \,|\, d_j=k) \Pr(d_j=k)\\\\ & = &\sum\limits_{k=1}^{N_2}\frac{k(k-1)}{N_2(N_2-1)}f(k). \end{array} $$

To calculate W ², the total expected number of shared inputs from all N ₁ neurons, we simply need to multiply by N ₁:

$$ W^2=\sum\limits_{k=1}^{N_2}\frac{N_1k(k-1)}{N_2(N_2-1)}f(k). $$

Dividing by W ¹ Eq. (21) gives the fraction of shared input parameter β in terms of the degree distribution

$$ \beta = \frac{\sum_{k=1}^{N_2}k(k-1)f(k)} {(N_2-1) \sum_{k=1}^{N_2} k f(k)}. $$

(22)

For a random network, the degree distribution is a binomial distribution

$$ f(k) = \binom{N_2}{k} p^k (1-p)^{N_{2}-k} $$

so that

$$ W^1=N_1p \qquad \text{and} \qquad \beta = p, \label{eq:Wbeta_rand} $$

(23)

agreeing with the results from Section 3.2.2. If the degree distribution is given by a power law with maximum degree d _max ≤ N ₂

$$ f(k) = \begin{cases} k^{-\gamma}/ \sum_{n=1}^{d_{\mathrm{max}}} n^{-\gamma}& \text{if $k \le d_{\mathrm{max}}$,}\\\\ 0 & \text{otherwise,} \end{cases} $$

the expressions for W ¹ and β become

$$\begin{array}{lll} \label{eq:Wbeta_pl} W^1 &=& \frac{N_1\sum_{k=1}^{d_{\mathrm{max}}} k^{1-\gamma}} {N_2 \sum_{k=1}^{d_{\mathrm{max}}} k^{-\gamma}} \qquad \text{and}\\\\ \beta &=& \frac{\sum_{k=1}^{d_{\mathrm{max}}}(k-1)k^{1-\gamma}} {(N_2-1) \sum_{k=1}^{d_{\mathrm{max}}} k^{1-\gamma}}. \end{array} $$

(24)

If the degree distribution is given by a Gaussian

$$ f(k) = \frac{ e^{-k^2/2\sigma^2}}{\sum_{n=1}^{N_2} e^{-n^2/2\sigma^2}} $$

the expressions for W ¹ and β become

$$\begin{array}{lll} \label{eq:Wbeta_g} W^1 &=& \frac{N_1\sum_{k=1}^{N_2} ke^{-k^2/2\sigma^2}} {N_2 \sum_{k=1}^{N_2} e^{-k^2/2\sigma^2}} ~~\text{and}\\\\ \beta &=& \frac{\sum_{k=1}^{N_2}k(k-1)e^{-k^2/2\sigma^2}} {(N_2-1) \sum_{k=1}^{N_2} ke^{-k^2/2\sigma^2}}. \end{array} $$

(25)

In our simulations, we set N ₁ = N ₂ = N and solved the first equation of Eq. (23), Eq. (24) or Eq. (25) for the p, γ or σ, respectively that gave the chosen W ¹. Then,we used the second equation of Eq. (23), Eq. (24) or Eq. (25) to calculate β.