Population density equations for stochastic processes with memory kernels

We present a method for solving population density equations (PDEs)–-a mean-field technique describing homogeneous populations of uncoupled neurons—where the populations can be subject to non-Markov noise for arbitrary distributions of jump sizes. The method combines recent developments in two different disciplines that traditionally have had limited interaction: computational neuroscience and the theory of random networks. The method uses a geometric binning scheme, based on the method of characteristics, to capture the deterministic neurodynamics of the population, separating the deterministic and stochastic process cleanly. We can independently vary the choice of the deterministic model and the model for the stochastic process, leading to a highly modular numerical solution strategy. We demonstrate this by replacing the master equation implicit in many formulations of the PDE formalism by a generalization called the generalized Montroll-Weiss equation—a recent result from random network theory—describing a random walker subject to transitions realized by a non-Markovian process. We demonstrate the method for leaky- and quadratic-integrate and fire neurons subject to spike trains with Poisson and gamma-distributed interspike intervals. We are able to model jump responses for both models accurately to both excitatory and inhibitory input under the assumption that all inputs are generated by one renewal process.

Figure 1

Dynamics of LIF (left) and QIF neurons (right). Time $t$ is shown as a function of $V_0$ on an interval $[V_{\min },V_{\max }]$.

Figure 2

A geometric grid for LIF neurons. (The bin $[0,V^{*}]$ is a fiducial bin used to avoid having exponentially many bins near $V=0$.)

Figure 3

A shift of probability mass through the geometric grid is sufficient to capture the evolution of the density profile due to the deterministic neuronal dynamics. For LIF neurons (top) mass moves in the direction of the equilibrium point. Mass that enters the fiducial bin surrounding the equilibrium point remains there. For QIF neurons with $I<0$ (bottom), there are two equilibrium points: one stable and one unstable. Movement towards the stable equilibrium is similar to the LIF case. Movement away from the unstable equilibrium towards the threshold is upwards. Probability mass reaching threshold must be removed from the system and reinserted at a potential $V_{\mathrm{reset}}$, possibly after observing a refractive period.

Figure 4

The shift of probability mass can largely be replaced by an index update. The black bars at the top indicate a geometric grid. The array V stores the values of the grid boundaries; the relationship between the contents of V and the grid boundaries are indicated symbolically by red lines. This relationship is immutable; it remains constant throughout simulation. P is an array representing the probability mass at time $t$. The top figure represents the situation at $t=0$: Each element of $P$ contains a numerical value representing the amount of probability mass. For a given element, the blue arrow indicates the potential interval containing this; together the V and P arrays represent a discretized density profile. The evolution of the density profile is realized by updating the relationship between the P and the V arrays, as indicated in the bottom panel representing the density profile at $t=\mathrm{\Delta }t$, by the change of the blue arrows. The contents of the V array remain unchanged, as do the contents of P, with the exception of two bins. The situation depicted in this figure shows decay towards a steady state that is represented by the third and fourth potential bins from the left. Mass represented by the two outward-pointing arrows on the extreme left and right represents mass that has cycled from the equilibrium bins to potential values at the extreme end of the potential interval. If the stationary point is stable, this is undesirable, and this mass should be removed and added to the elements of the mass array that currently point to the equilibrium bins; this is indicated by the dashed arrows. This figure shows a neural model that has a single stable equilibrium point, such as the LIF model. For a model with more than one stationary point, such as the QIF neuron for $I<0$, the V and P arrays must be separated into isolated strips, each with their own relationship between the P and V arrays. By updating the relationship between elements of P and V, the shift of mass is captured almost entirely without moving data around.

Figure 5

Coefficients for the Poisson master equation are purely determined by synaptic efficacy $h$ and the bin boundaries. Two grids are shown above, with the bottom displaced by $h$. The transition matrix is determined by what fraction the displaced grid covers each bin of the original grid. For example, the bin labeled “bin $k$” has length $L_k$. After displacement by $h$, it overlaps bin $k$ by some interval $l_{k,k}$ and bin $k+1$ by an interval $l_{k,k+1}$. Therefore, the transition matrix would have entries $m_{k,k}=\frac{l_{k,k}}{L_k},m_{k,k+1}=\frac{l_{k,k+1}}{L_k},$ and so on. At the same time, part of bin $k+1$ is pushed above threshold, which corresponds to spiking. The probability mass pushed above threshold is placed in the bin containing $V_{\mathrm{reset}}$ at the next time step (or after a delay if a refractory period is desired). This procedure for bin $k$ allows us to calculate the $k\mathrm{th}$ row of the transition matrix, so repeating the procedure for all bins gives us the full matrix.

Figure 6

Firing rates of the LIF neuron with inputs from a $\mathrm{\Gamma }(\alpha ,\nu )$ distribution. Lines are calculated using our method, while markers are from Monte Carlo simulations of $10000$ neurons. $h=0.03,V_{\mathrm{th}}=1$ in all cases. Solid line and crosses: $\alpha =1,\nu =800$, i.e., a Poisson process with rate 800. Dashed line and circles: $\alpha =2,\nu =1600$. Dotted line and triangles: $\alpha =3,\nu =2400$. ($\nu$ is varied such that the expectation of the input process remains the same across all cases.)

Figure 7

(a) Firing rates of the LIF neuron with inputs from a $\mathrm{\Gamma }(\alpha ,\nu )$ distribution. Lines are calculated using our method, while markers are from Monte Carlo simulations of $10000$ neurons. $h=0.1$ in all cases. Solid line and crosses: $\alpha =1,\nu =150$, i.e., a Poisson process with rate 150. Dashed line and circles: $\alpha =2,\nu =300$. Dotted line and triangles: $\alpha =3,\nu =450$. ($\nu$ is varied such that the expectation of the input process remains the same across all cases.) (b) The steady-state density profiles for the different shape factors.

Figure 8

Steady-state density of the generalized OU process. Lines are calculated using our method, while markers are from Monte Carlo simulations of $20000$ neurons. $h=0.1$ in all cases. Solid line and crosses: $\alpha =1,\nu =10$, i.e., a Poisson process with rate 10. Dashed line and circles: $\alpha =2,\nu =20$. Dotted line and diamonds: $\alpha =3,\nu =30$.

Figure 9

Gain curves for $h=0.15$. Solid line: shape = 1 (Poisson process). Dashed line: shape = 2. Dotted line: shape = 3. Input firing rate is expected input $\nu /\alpha$.

Figure 10

The transition matrix $\mathcal{M}$ as an operator for moving probability mass from bin $i$ to bin $j$. Synaptic noise removes neurons from their current position, resulting in a loss term along the diagonal. Neurons undergoing an excitatory jump move up in potential and thereby end up in the mass array with a higher bin number. Neurons undergoing inhibition end up at a lower bin number, at a larger distance from the diagonal, reflecting $h_i=-4h_e$. The complex shape of the two bands is a result of using a geometric grid: Near the reversal potential the same jump in potential covers more bins. The reversal bin is larger than the neighboring geometric bins so on translation covers a large number of them. This represents the straight part of the bands.

Figure 11

Firing rates of the LIF neuron with inputs from a $\mathrm{\Gamma }(\alpha ,\nu )$ distribution. Lines are calculated using our method, while markers are from Monte Carlo simulations of $10000$ neurons. In all cases, an input spike has an 0.8 probability of being excitatory ($h=0.05$) and an 0.2 probability of being inhibitory ($h=-0.2$). Solid line and crosses: $\alpha =1,\nu =2000$, i.e., a Poisson process with rate 2000. Dashed line and circles: $\alpha =2,\nu =4000$. Dotted line and triangles: $\alpha =3,\nu =6000$. ($\nu$ is varied such that the expectation, $\nu /\alpha$, of the input process remains the same across all cases.)

Figure 12

(a) Firing rates of the QIF neuron with inputs from a $\mathrm{\Gamma }(\alpha ,\nu )$ distribution. Lines are calculated using our method, while markers are from Monte Carlo simulations of $10000$ neurons. The population has a time constant $\tau =0.01$ and a constant current $I=-1$, and the stochastic input has a synaptic efficacy $h=0.2$ in all cases. Solid line and crosses: $\alpha =1,\nu =500$, i.e., a Poisson process with rate 500. Dashed line and circles: $\alpha =2,\nu =1000$. Dotted line and triangles: $\alpha =3,\nu =1500$. ($\nu$ is varied such that the expectation of the input process remains the same across all cases.) (b) The steady-state density profiles for the different shape factors. The membrane potential has been renormalized so the $V_{\mathrm{th}}=1.0$ and is dimensionless.

Figure 13

The firing rate of a LIF population with two input processes, as a function of excitatory synaptic efficacy. The solid lines are the simulation of Poisson noise with our method, while the dash-dotted line is the solution of the Montroll-Weiss equation for the gamma processes. The Monte Carlo simulations are the triangles and crosses with error bars. $\tau =0.05, \nu =500$. The inhibitory synaptic efficacy $h_I$ is fixed at 0.15.