Exact firing time statistics of neurons driven by discrete inhibitory noise

TL;DR

This paper derives exact statistical measures for neuron firing times under inhibitory noise, accounting for both small and large synaptic inputs, and applies to heterogeneous neuron populations.

Contribution

It provides exact analytical results for firing rate, variability, and spectrum of neurons influenced by inhibitory synaptic inputs, including finite amplitude effects.

Findings

Exact firing rate and variability formulas derived

Results applicable to finite and large synaptic inputs

Methods extend to heterogeneous neuron populations

Abstract

Neurons in the intact brain receive a continuous and irregular synaptic bombardment from excitatory and inhibitory pre-synaptic neurons, which determines the firing activity of the stimulated neuron. In order to investigate the influence of inhibitory stimulation on the firing time statistics, we consider Leaky Integrate-and-Fire neurons subject to inhibitory instantaneous post-synaptic potentials. In particular, we report exact results for the firing rate, the coefficient of variation and the spike train spectrum for various synaptic weight distributions. Our results are not limited to stimulations of infinitesimal amplitude, but they apply as well to finite amplitude post-synaptic potentials, thus being able to capture the effect of rare and large spikes. The developed methods are able to reproduce also the average firing properties of heterogeneous neuronal populations.

Equations

$$\frac{dv}{dt}=-\frac{v-{\mu }_0}{\tau }+I(t).$$

$$I(t)=\sum _{\{t_e\}}{a}_e\delta (t-t_e)+\sum _{\{t_i\}}{a}_i\delta (t-t_i),$$

$$\frac{\partial P}{\partial t}+\frac{\partial J}{\partial v}=r(t)\left[\delta (v-v_{re})-\delta (v-v_{th})\right]+\delta (t)\delta (v-v_{re}).$$

$$J=-\left(\frac{v-{\mu }_0}{\tau }\right)P+J_e+J_i;$$

$$J_x=R_x{\int }_{-\infty }^{v}dwP(w,t){\int }_{v-w}^{\infty }da{A}_x(a).$$

$$J_e(v_{th},t)=r(t),J_i(v_{th},t)=0,P(v_{th},t)=0;$$

$${\int }_{-\infty }^{v_{th}}P(v,t)dv=1.$$

$$\frac{\partial J_e}{\partial v}+\frac{J_e}{a_e}$$

$${\tilde{J}}_e(s,t)=\tilde{P}(s,t)Q_e(s)-r(t)S_e(s).$$

$$Q_e(s)=\frac{R_e}{1-s{a}_e}{S}_e(s)=\frac{{\mathrm{e}}^{s{v}_{th}}}{1-s{a}_e}.$$

$$\frac{d{Z}_0}{ds}=\tau \left(\frac{{\mu }_0}{\tau }{Z}_0+{\tilde{J}}_e+{\tilde{J}}_i\right);$$

$$\frac{1}{Z_0}\frac{d{Z}_0}{ds}$$

$$Z_0$$

$$H_e=\{\begin{array}{ll}(1-s{a}_e{)}^{-\tau R_e}& \text{mbox}for{A}_e(a)=\Theta (a)\exp (-a/a_e)/a_e\\ 1& \text{mbox}forDCcurrent\end{array}$$

$$\frac{d{\tilde{P}}_0}{ds}={\tilde{P}}_0\tau \left(\frac{{\mu }_0}{\tau }+Q_e+Q_i\right)-\frac{r_0}{s}\left(e^{s{v}_{re}}-e^{s{v}_{th}}+s{S}_e\right).$$

$$\frac{1}{Z_0}\frac{d\tilde{P}}{ds}=\frac{\tilde{P}}{Z_0^2}\frac{d{Z}_0}{ds}-\tau \frac{r_0}{s{Z}_0}F(s).$$

$$\frac{\tilde{P}(s)}{Z_0(s)}=\tau r_0{\int }_s^{\bar{x}}\frac{dc}{c{Z}_0(c)}F(c)+\frac{\tilde{P}(\bar{x})}{Z_0(\bar{x})}.$$

$$\frac{1}{\tau r_0}={\int }_0^{\bar{x}}\frac{dc}{c{Z}_0(c)}F(c),$$

$$\bar{x}=\{\begin{array}{ll}1/a_e& \text{mbox}for{H}_e=(1-s{a}_e{)}^{-\tau R_e}\\ \infty & \text{mbox}for{H}_e=1.\end{array}$$

$$\frac{d\hat{P}(s,\omega )}{ds}-\tau \left(\frac{\mu }{\tau }+Q_e+Q_i-\frac{i\omega }{s}\right)\hat{P}(s,\omega )=-\frac{\hat{\rho }(\omega )}{s}F(s)+\frac{e^{s{v}_{re}}}{s};$$

$$\hat{\rho }(\omega )=\frac{{\int }_0^{\bar{x}}ds{s}^{i\omega \tau }\frac{d}{ds}A(s)}{{\int }_0^{\bar{x}}ds{s}^{i\omega \tau }\frac{d}{ds}B(s)};$$

$$A(s)=\frac{e^{s{v}_{re}}}{Z_0};B(s)=\frac{F(s)}{Z_0}.$$

$$\hat{q}(\omega )=\frac{\hat{\rho }(\omega )}{1+\hat{\rho }(\omega )};\hat{C}(\omega )=r_0\frac{1-|\hat{q}(\omega ){|}^2}{(1-\hat{q}(\omega ){)}^2}+2\pi r_0^2\delta (\omega ),$$

$$\frac{{\partial }^n\hat{q}}{\partial {\omega }^n}{|}_{\omega =0}=(-\mathrm{i}{)}^n⟨t^n⟩.$$

$$\hat{\rho }(\omega )\approx \frac{n_0+n_1\omega +n_2{\omega }^2}{d_0+d_1\omega +d_2{\omega }^2};$$

$$n_1=i{\int }_0^{\bar{x}}ds\log s\frac{d}{ds}A(s)d_2={\int }_0^{\bar{x}}ds\frac{\log s}{s}B(s).$$

$$⟨t⟩=-i{d}_1=1/r_0⟨t^2⟩=2\frac{d_1^2-d_2{n}_0+d_1{n}_1}{n_0^2};$$

$$\text{CV}=\frac{\sqrt{⟨t{⟩}^2-⟨t^2⟩}}{⟨t⟩}.$$

$$A_i(a)=-\frac{\Theta (-a)}{a_i}{\mathrm{e}}^{-a/a_i}.$$

$$\frac{\partial J_i}{\partial v}+\frac{J_i}{a_i}=R_{i}P(v,t).$$

$${\tilde{J}}_i=\frac{R_i}{1-s{a}_i}\tilde{P}.$$

Full text

Introduction

Neurons in the neocortex in vivo are subject to a continuous synaptic bombardment reflecting the intense network activity [1]. In the so-called high-input regime, in which neurons receive hundreds of synaptic inputs during each interspike interval [2], the firing statistics of model neurons is usually obtained in the context of the diffusion approximation (DA) [3, 4].

Within such an approximation the post-synaptic potentials (PSPs) are assumed to have small amplitudes and high arrival rates, therefore the synaptic inputs can be treated as a continuous stochastic process characterized simply by its average and variance, while the shape of the distribution of the amplitudes of the PSPs is irrelevant [5]. However several experimental studies have revealed that rare PSPs of large amplitude can have a fundamental impact in the network activity [6, 7]. Furthermore, the experimentally measured synaptic weight distributions display, both for excitatory and inhibitory PSPs, a long tail towards large amplitudes and a peak at low amplitudes [8, 9, 6, 10, 7, 11, 12]. The effect of rare and large excitatory post-synaptic potentials (EPSPs) has been recently examined for generalized leaky integrate-and-fire (LIF) models with generic EPSP distributions [13] and in balanced sparse networks for conductance based LIF neurons with log-normal EPSP distributions [14]. The presence of few strong synapses induce faster and more reliable responses of the network even for small inputs [13, 14]. Interestingly, in [14] it has been shown that a single neuron driven by random synaptic inputs log-normally distributed reveals a clear aperiodic stochastic resonance [15, 16], which is not evident for Gaussian distributed EPSPs.

However, even for the simple case of LIF neurons exact analytic results are still lacking for large PSPs with generic synaptic weight distributions, apart for the case of the exponentially distributed PSPs reported in [17]. In particular, Richardson and Swarbrick have been able to obtain the statistics of interspike interval (ISI) for LIF neurons receiving balanced excitatory and inhibitory Poissonian spike trains with exponentally distributed synaptic weights [17]. Furthermore, results for generic EPSP distributions have been obtained in  [13] by developing a semi-analytic approach to solve the continuity equation for the membrane potential distribution.

In this paper, we report exact analytic results for the firing time statistics of neurons receiving inhibitory Poissonian spike trains for various synaptic weights distributions. Namely, we estimate the firing time statistics for LIF neurons subject to inhibitory post-synaptic potentials (IPSPs) (with instananeous rise and decay time) characterized by constant amplitude, as well as for uniform and truncated Gaussian IPSP distributions. Furthermore, we apply the developed formalism to sparse inhibitory networks with heterogeneous neuronal properties.

Models and Methods

Model and population-based formalism

We will consider the firing statistics of a LIF neuron  [18, 19] subject to a constant external DC current $\mu_{0}$ and to a synaptic drive $I(t)$, in this case the dynamical evolution of the membrane potential $v$ is given by the following equation

[TABLE]

where $\tau=20$ ms is the membrane time constant. The neuron fires whenever the membrane potential reaches the threshold value $v_{th}=10$ mV, afterwards the potential is reset to the value $v_{re}=5$ mV. The synaptic current $I(t)$ accounts for the linear superposition of the instantaneous excitatory and inhibitory PSPs and it can be written as

[TABLE]

where $a_{x}$ denote the amplitudes of EPSPs $x=e$ ($a_{e}>0$) and IPSPs $x=i$ ($a_{i}<0$), while the variables $t_{x}$ represent their respective arrival times, which are assumed to be Poissonian distributed with rates $R_{x}$. Eqs. (1) and (2) are equivalent to the Stein’s model [20] with pulse amplitudes randomly drawn from distributions $A_{x}(a)$. his model, despite its extreme simplicity, has been shown to be able to provide optimal predictions for the average firing rates and spike times of cortical neurons [21, 22].

The aim of this paper is to provide exact analytic expressions for the first two moments of the stationary firing statistics, namely the average firing rate $r_{0}$ and the associated coefficient of variation CV, as well as for the spike-train spectrum (STS) $\hat{C}(\omega)$ [23]. To obtain such results we follow the approach developed in [17], in particular within a population-based formalism we introduce the probability density $P(v,t)$ of the membrane potentials together with the associated flux $J(v,t)$. The continuity equation relating these two quantities can be written as:

[TABLE]

On the r.h.s. of the above equation are reported the sink (source) term for the neuronal population associated to the membrane threshold (reset), with $r(t)$ being the instantaneous firing rate of the population. The last term on the r.h.s. takes into account the initial distribution of the membrane potentials, which are assumed to be all equal to the reset value at $t=0$ [24].

The flux $J(v,t)$ can be decomposed in three terms as follows

[TABLE]

where the first term on the r.h.s is the average drift, while $J_{e}$ ($J_{i}$) represents the excitatory (inhibitory) fluxes originating from the Poissonian synaptic drives. The fluxes can be written as a convolution of the distribution of the membrane potentials with the synaptic amplitude distribution, namely

[TABLE]

The previous set of equations is complemented by the following boundary conditions:

[TABLE]

and by the requirement that the membrane potential distribution is properly normalized at any time, i.e

[TABLE]

Analytical method to obtain the exact firing time statistics

The estimation of the firing statistics for a LIF neuron subject to shot noise has proven to be a problem analytically hard to solve [25, 26]. The reason is related to the overshoots over the threshold $v_{th}$ induced by the finite amplitude of the PSPs, which renders difficult the estimation of the membrane potential distribution. However, it is well known that one of the few cases in which the first passage time problem can be solved, is represented by exponentially distributed PSP weights, thanks to the memory-less property associated with exponential distributions [27]. Richardson and Swarbrick [17] made use of this unique property to derive the exact solution of the firing rate for the case in which both inhibitory and excitatory kick amplitudes are exponentially distributed. The fact that the only boundary relevant for the first passage time is $v_{th}$, and considering that no trajectory can cross it from above, implies that inhibitory kicks do not contribute to the overshoot and therefore no restriction over the distribution of their amplitudes should be in principle imposed in order to obtain an analytic solution of the problem.

In the following, using the Laplace transform method, we will derive the analytic expressions for the firing rates, the coefficient of variation and for the spike-train spectrum for various distributions $A_{i}(a)$ of the inhibitory amplitudes. For what concerns the excitatory synaptic input, we will limit our investigation to two analytic solvable cases: namely, to exponentially distributed synaptic weights, where $A_{e}=\Theta(a)\exp(-a/a_{e})/a_{e}$, and to constant excitatory synaptic drive, encompassed in an external DC current $\mu_{0}>v_{th}$.

Let us first consider the excitatory term for exponentially distributed $a_{e}$, in this case the integral equation for the excitatory flux Eq. (5) can be rewritten in a differential form as

[TABLE]

where the last term on the r.h.s. of the above equation accounts for the absorbing boundary condition at threshold. The bidirectional Laplace transform (from now on only Laplace transform) $\tilde{f}(s)=\int_{-\infty}^{\infty}{\rm e}^{sv}f(v)dv$ of Eq. (7) can be written as a combination of linear functions in $\tilde{P}(s,t)$ and $r(t)$, namely

[TABLE]

For the particular case of the exponentially distributed EPSP amplitudes

[TABLE]

Whenever the excitatory input is simply given by the DC current $\mu_{0}$, we will assume $J_{e}=Q_{e}=S_{e}=0$ and apply the same formulation that we will expose in the following.

For the distribution of inhibitory amplitudes, the only restriction that we will impose is that, one should be able to write the Laplace transform of the inhibitory flux as a linear function of the probability density function, namely $\tilde{J}_{i}(s,t)=\tilde{P}(s,t)Q_{i}(s)$.

Steady State Firing Rate

Under the above assumptions, we can estimate the Laplace transform of the sub-threshold voltage distribution $Z_{0}\equiv Z_{0}(s)$, which corresponds to the generating function for the sub-threshold voltage moments. Therefore, $Z_{0}$ can be estimated as the Laplace transform of $P(v,t)$ when $v_{th}\to\infty$, implying also that $J=r(t)=0$. In particular, by taking the Laplace transform of Eq. (4) together with the assumption that $J=0$, we obtain

[TABLE]

where we set $Z_{0}=\tilde{P}$.

Since $\tilde{J}_{e}$ and $\tilde{J}_{i}$ are linear in $\tilde{P}$, we can rewrite Eq. (10) and solve it as:

[TABLE]

where the excitatory contribution is encompassed in the term

[TABLE]

Once we have calculated the generating function, we can solve the stationary case corresponding to $\partial P/\partial t=0$, performing the Laplace transform of the continuity equation (3), which reads as

[TABLE]

In our notation, the variables with a zero subscript denote stationary quantities. As expected, for $r_{0}=0$ the function $\tilde{P}_{0}$ satisfies Eq. (11), therefore we can rewrite the previous equation as

[TABLE]

where we have indicated with $F(s)$ the function multiplying the term $r_{0}/s$ in the r.h.s of Eq. (14). The straightforward solution of Eq. (15) is

[TABLE]

Notice that, although we have derived an expression for $Z_{0}$, we have no knowledge of the functional form of $\tilde{P}$. Whenever it is possible to identify an integration limit $\bar{x}$, where the term $1/Z_{0}(\bar{x})$ vanishes, the following exact analytic expression for the stationary firing rate can be obtained

[TABLE]

where we have made use of the normalization condition of the probability densities, i.e. $\tilde{P}_{0}(0)=Z_{0}(0)=1$.

Our analysis is limited to the two previously reported types of excitatory drive, because in these two cases an integration limit $\bar{x}$, where $Z_{0}^{-1}(\bar{x})=0$, can be easily found to be

[TABLE]

It should be remarked that in presence of both sources of excitatory drive, the integration limit can still be identified whenever $\mu_{0}<v_{th}$ and it corresponds to the first one in (18), while we have been unable to solve the case when both, the excitatory drift and the excitatory spike train, can lead the neuron to fire (see conclusions section for a discussion of this point).

First and Second Moment of the First Passage Time Distribution

Let us now focus on the time dependent evolution of the continuity equation, in this case, the equation can be solved by performing the Fourier transform in time and the Laplace transform in the membrane potential of Eq. (3) and (4). Namely, we obtain

[TABLE]

where $\hat{f}(s,\omega)=\int_{-\infty}^{\infty}{\rm e}^{-2\pi i\omega t}dt\int_{-\infty}^{\infty}{\rm e}^{sv}f(v,t)dv$, and $\hat{\rho}(\omega)$ is the Fourier transform of the spike-triggered rate. Dividing both sides by $Z_{0}$ and integrating the functions over the interval $s=[0,\bar{x}]$, the l.h.s of Eq. (19) vanishes and an analytic expression for $\hat{\rho}(\omega)$ can be obtained, namely

[TABLE]

where

[TABLE]

Equation (20) diverges exactly at $\omega\equiv 0$, and in that case it should be complemented with the expression $\hat{\rho}(\omega=0)=r_{0}\pi\delta(\omega)$. The spike-triggered rate (also called the conditional firing rate) provides all the information on the spike train statistics. For instance, the spike train spectrum $\hat{C}(\omega)$, which is the Fourier transform of the auto-correlation function of the spike train, is related to the spike triggered rate via the formula $\hat{C}(\omega)=r_{0}[1+2\text{Re}(\hat{\rho}(\omega))]$ [23].

Moreover, $\hat{C}(\omega)$ and $\hat{\rho}(\omega)$ are also related with the Fourier transform of the first-passage time density $\hat{q}(\omega)$ as follows:

[TABLE]

where $\lim_{\omega\to 0}\hat{C}(\omega)=r_{0}$ and $\lim_{\omega\to\infty}\hat{C}(\omega)=\text{CV}^{2}r_{0}$  [28].

Therefore the $n-th$ moment of the first passage time distribution is given by

[TABLE]

For the estimation of the first two moments it is sufficient to expand to the second order in $\omega$ the terms entering in Eq. (20), in particular $s^{i\omega\tau}\approx 1+i\tau\omega\log(s)-\frac{1}{2}(\omega\tau\log s)^{2}+\mathcal{O}(\omega^{3})$. Finally, $\hat{\rho}(\omega)$ can be approximated as a polynomial in $\omega$, that reads as

[TABLE]

where $n_{0}=-1$, $d_{1}=i/r_{0}$, $d_{0}=0$ and

[TABLE]

It is worth to notice that, despite the expansion is limited to the second order, the solutions for the first and second moments are exact since terms of higher order disappear when evaluated at $\omega=0$.

From the expression (23) it is easy to verify that the first and second moments of the first passage time distribution are given by

[TABLE]

and from these it is straightforward to estimate the coefficient of variation CV of the ISI, namely

[TABLE]

As can be seen by this general solution, the two relevant quantities that define uniquely the stationary firing statistics are the functions $Z_{0}(s)$ and $F(s)$ defined in the Laplace space, where $F(s)$ depends only on the considered excitatory drive, while $Z_{0}(s)$ on the whole sub-threshold input.

Explicit expressions of $Z_{0}$ for selected distributions

In this manuscript, we focus on the exact solutions for four relevant types of distributions of the inhibitory synaptic weights $A_{i}(a)$: namely, exponential (ED), $\delta$ (DD), uniform (UD) and truncated Gaussian distribution (TGD). The exact expressions for $Z_{0}(s)$ for each of these four cases are reported in the following, together with the average and variance of the corresponding distributions $A_{i}(a)$.

Exponential Distribution (ED): In this case the synaptic weight distribution is of the form:

[TABLE]

For this distribution, the equations for the inhibitory flux (Eq. (5)) can be rewritten in a differential form analogous to Eq. (7), namely

[TABLE]

where the term accounting for the absorbing boundary is not present since we are considering the inhibitory shot noise.

The Laplace transform of this equation reads as

[TABLE]

For this choice of distribution, it is possible to obtain an explicit expression for the generating function the sub-threshold voltage moments, namely

[TABLE]

The mean and variance due to the distribution of the inhibitory synaptic weights for the exponential case are $\langle a_{i}\rangle=a_{i}$ and $var[a_{i}]=a_{i}^{2}$; while the absolute value of the skewness is two.

$\delta$** Distribution (DD):** In the case in which the inhibitory population delivers spikes with a constant amplitude $a_{i}$, the distribution becomes simply a $\delta$-function:

[TABLE]

In this case the equation for the inhibitory flux becomes the following

[TABLE]

and the associate Laplace transform reads as

[TABLE]

For this simple distribution the generating function $Z_{0}$ reads as:

[TABLE]

where $E_{I}(y)=\int_{y}^{\infty}dye^{y}/y$ is the exponential integral and $C_{0}$ accounts for the normalization condition requiring that $Z_{0}(0)=1$ and its explicit expression is

[TABLE]

where $\Gamma\approx 0.577731$ is the Euler-Mascheroni constant.

For this distribution the mean and variance of the inhibitory process are simply $\langle a_{i}\rangle=a_{i}$ and $var[a_{i}]=0$, and also the skewness is zero.

Uniform Distribution (UD): For uniformly distributed synaptic weights with support $[l_{1}\;l_{2}]$, we can express the UD as follows

[TABLE]

The variation of the flux respect to the potential is given by:

[TABLE]

and the corresponding Laplace transform reads as

[TABLE]

This leads to the generating function together with the normalization constant:

[TABLE]

The mean and variance of the inhibitory shot noise are $\langle a_{i}\rangle=(l_{2}+l_{1})/2$ and $var[a_{i}]=(l_{2}-l_{1})^{2}/12$, respectively; while the skewness is zero.

Truncated Gaussian Distribution (TGD): A biologically relevant distribution is the Gaussian distribution [38], which is peaked at $a_{p}$ and with a standard deviation equal to $\sigma_{G}$. For notation simplicity we write the Gaussian distribution as

[TABLE]

Since we are only interested in the inhibitory kicks, we truncate the original distribution and we impose the support $(-\infty,0]$. The distribution of the synaptic weights can be written as

[TABLE]

$\Phi(0)$ is the cumulative distribution of the normal distribution evaluated at the upper limit of the support according to the equation

[TABLE]

The flux of inhibitory probability takes the form

[TABLE]

The corresponding Laplace transform is:

[TABLE]

The complicated expression appearing in Eq. (46) does not allow us to obtain the explicit expression for $Z_{0}$. However it is easy to integrate numerically Eq. (12) once we know $\tilde{J}_{i}$. When dealing with a width of the Gaussian distribution $\sigma_{G}>1$, the exponential term in Eq. (46) can grow very rapidly while the complementary error function tends to 0, generating numerical problems due to machine precision, specially in the evaluation of the terms $s>1$. In such cases one can use the first order expansion of the complementary error function: $\text{erfc}(y)\approx\exp(-y^{2})/\sqrt{\pi}y$ when $y>>1$ and Eq. (46) can be simplified in such cases to:

[TABLE]

For the TGD, the mean value and variance of the membrane potentials associated to the inhibitory shot noise take the following form

[TABLE]

The value of the negative skewness grows with $\sigma_{G}$, namely it passes from $\simeq-0.22$ for $\sigma_{G}=|a_{i}/2|$ to $\simeq-0.92$ for $\sigma_{G}=|5a_{i}|$ and it tends to the value $-\sqrt{2}(4-\pi)/(\pi-2)^{3/2}\simeq-0.9992$ for $\sigma_{G}\to\infty$.

Effective Input and Synaptic Noise Intensity

In order to perform meaningful comparisons between the neuronal response for different synaptic weights distributions, we will consider the responses obtained for the same effective average input $\mu_{T}$ and noise intensity $\sigma$.

For a neuron receiving an inhibitory Poissonian spike train at a rate $R_{i}$ and with an average synaptic weight $\langle a_{i}\rangle$ the effective average input is

[TABLE]

where $a_{i}<0$ and $\mu_{e}$ is the average excitatory input. When only the drift term is present $\mu_{e}=\mu_{0}$, while for a neuron receiving also an excitatory Possonian spike train of rate $R_{e}$ and with average synaptic weight $\langle a_{e}\rangle$, it becomes $\mu_{e}=\mu_{0}+\tau R_{e}\langle a_{e}\rangle$. On the one hand, for an effective sub-threshold input $\mu_{T}<v_{th}$ the dynamics of the LIF neuron is characterized by two timescales: the relaxation time from the reset value $v_{re}$ to the resting value $\mu_{T}$ and the activation time associated to the escape process from the resting state to the threshold induced by the fluctuations in the input [28]. On the other hand, for a supra-threshold LIF neuron for which $\mu_{T}>v_{th}$, in absence of a refractory state, the only characteristic time is the tonic firing rate

[TABLE]

Moreover, in the set-up that we are studying, the neuron is in general subjected to two sources of randomness. A first source due to the variability of the arrival times of the Poisson process, and a second source due to the distribution of the amplitudes of the synaptic weights. Therefore, the total noise intensity $\sigma^{2}$ associated to these two uncorrelated processes is given by the sum of the variances of each process, namely

[TABLE]

From Eqs. (50) and (52) one can see that the values of $\mu_{T}$ and $\sigma$ can be independently tuned with an appropriate selection of the parameter $\mu_{e}$ and $R_{i}$ for any fixed average inhibitory amplitude. Whenever the excitatory drive is given simply by a DC term $\mu_{0}$, for any choice of the couple of values {$\mu_{T},\;\sigma\}$ one has an unique solution, which however does not always correspond to a firing activity. For instance, at small values of $\sigma^{2}$ one can find values of $\mu_{0}<v_{th}$, meaning that the neuron will never fire. This implies the existence of a firing onset threshold for this class of systems which depends on the shape of the distribution. Conversely, when exponentially distributed excitatory kicks are present, two parameters (amplitude and the rate of arrival of the excitatory kicks) define $\mu_{e}$ and therefore a whole family of solutions exist for any fixed $a_{i}$. Therefore in this case there is not a finite firing onset threshold, allowing for arbitrarily small rate even for small noise intensities.

Throughout this paper we will compare the results obtained within the shot noise framework against the widely used diffusion approximation [3]. Within such approximation the particular shape of the synaptic weight distributions are irrelevant and indeed the only relevant quantities defining the stationary firing rate and the CV are $\mu_{T}$ and $\sigma^{2}$ through the formulas

[TABLE]

where $y_{\theta}=(v_{th}-\mu_{T})/\sigma$ and $y_{r}=(v_{re}-\mu_{T})/\sigma$.

By following [17], Eqs. (53) and (54) can be obtained from the shot noise formulation, namely from Eq. (26) and Eq. (22) by considering the expansion of $Z_{0}$ limited to the first two voltage moments:

[TABLE]

Results

In this Section we will apply the developed formalism to estimate the response of a single LIF neuron as well as the firing characteristics of a sparse inhibitory neural network for different synaptic weight distributions. In particular, to verify the limits of applicability of the approach we will compare the theoretical estimations with numerical data and with the DA.

Usually, the firing time statistics of a neuron subject to a noisy uncorrelated input is theoretically estimated within the so-called DA [3, 29, 4]. This approximation is only valid, however, when the PSP amplitudes are small compared with the reset-threshold voltage distance and the arrival frequencies are sufficiently high. Outside of such limits the DA fails to reproduce the numerical data and in particular it is unable to capture the differences due to different synaptic weight distributions [17, 13]. Furthermore, the DA has been employed to reproduce network activity of recurrent networks, in such a case one should assume that the spike trains impinging on the neuron are temporally uncorrelated, a condition usually fulfilled in sparse networks [30, 29]. However this should be considered as a first approximation, indeed correlations are present even in sparse balanced networks and they can be captured by driving a single neuron with a colored Gaussian noise self-consistently generated [31].

Influence of the IPSP distributions on the firing statistics

As a first aspect, let us consider the dynamics of a single LIF neuron subject to an excitatory DC current $\mu_{e}=\mu_{0}>v_{th}$ plus the inhibitory contribution given by a Poissonian train of IPSPS with constant amplitude $a_{i}$. In particular, we examine the neuronal response by increasing the noise intensity $\sigma^{2}$ at a constant effective input current, namely $\mu_{T}=9$ mV, corresponding to a sub-threshold case. In particular, the noise intensity is varied by adjusting the rate of the inhibitory train $R_{i}$. The results, reported in Fig. 1 (a), show that the neuron fires only for sufficiently large noise and the firing rate $r_{0}$ increases with $\sigma^{2}$ as expected. Furthermore, as shown in Fig. 1 (b) the coefficient of variation exhibits a clear minimum at an intermediate noise amplitude $\sigma_{m}^{2}$, an effect known as coherence resonance (CR) and widely studied in the context of excitable systems [32]. The emergence of CR is related to the presence of at least two competing time scales which depend differently on the noise amplitude, for the sub-threshold LIF neuron these two time scales are the relaxation and escape times [28].

For sufficiently small IPSP amplitudes the agreement between the DA, the numerical results and the analytic expression reported in Eq. (17) is almost perfect as evident from Fig. 1 (a-b), where the blue symbols/curve refer to $|a_{i}|=0.1$ mV. However, by increasing the IPSP amplitude the agreement between numerical results and the estimation given by the shot-noise result (17) remains very good, while the DA is unable to capture the effect of large IPSP. In particular, for $|a_{i}|=1$ mV the DA fails in reproducing the onset of the firing activity as well as the position and height of the minimum of the CV, and in general the firing statistics for low noise amplitudes (see red curves/symbols in Fig. 1 (a-b)). Unitary IPSPs of amplitude $\simeq 1$ mV have been measured experimentally in the hippocampus of Guinea-pig in vitro [8].

As shown in Fig. 1 (a-b), the increase of the IPSP amplitude has a noticeable effect on the neuronal response, even if the effective input $\mu_{T}$ and the noise intensity are the same as in the case $|a_{i}|=0.1$ mV. Increasing IPSP amplitudes leads to a decrease in the firing activity and it induces an increase of the maximal coherence observable at $\sigma_{m}^{2}$, which now occurs at larger noise amplitude with respect to the case $|a_{i}|=0.1$ mV. This can be explained by the fact that the relaxation times to the equilibrium value $\mu_{T}$ are longer for larger IPSP and this induces a deeper minimum in the CV as reported in [33]. Furthermore the irregularity in the emitted spikes, as measured by the CV, increases for $\sigma^{2}<\sigma_{m}^{2}$ and becomes more regular at larger noise intensities.

The effect of the IPSP amplitude can be better appreciated by performing a different test, namely by maintaining constant both the noise intensity and $\mu_{T}$ while increasing $|a_{i}|$. These two quantities can be independently tuned with a suitable selection of $R_{i}$ and $\mu_{0}$ in Eqs. (50) and (52). The results of this analysis are reported in Fig. 1 (c-d), the firing rate exhibits a dramatic decrease for increasing $|a_{i}|$, as expected due to the increase of average inhibitory current $\mu_{i}$. This effect is well reproduced by Eq. (17), but it is absolutely not captured by the DA as shown in Fig. 1 (c). The increase of the IPSP amplitude leads to a small variation of the CV revealing a minimum at an intermediate value $|a_{i}|\simeq 0.9$ mV (see Fig. 1 (d)). Nevertheless, the DA provides a constant value for the CV in the whole examined range. The origin of the minimum can be understood observing Fig. 1 (b): by increasing $|a_{i}|$ the overall minimum of the CV curve shifts towards larger noise amplitudes, and at the same time the CV values decrease (increase) for $\sigma^{2}>\sigma_{m}^{2}$ ($\sigma^{2}<\sigma_{m}^{2}$). In Fig. 1 (d), we consider a noise intensity $\sigma^{2}=2$ mV2 for all the simulations. For small (large) $|a_{i}|$ the maximal coherence is observable at $\sigma_{m}^{2}<2$ mV2 ($\sigma_{m}^{2}>2$ mV2), thus the CV value at $\sigma^{2}=2$ mV2 decreases (increases) with $|a_{i}|$. The minimum in Fig. 1 (d) occurs exactly when the CV displays its absolute minimum at $\sigma_{m}^{2}=2$ mV2.

For the moment we have considered only the case of constant IPSP amplitudes, now we will examine the influence of different distributions $A_{i}$ on the response of the single LIF neuron. In general, we observe that the shape of the IPSP distribution can noticeably influence the firing rate and the CV. In order to verify this observation, we consider the firing statistics of a neuron subject to the same average input and noise intensity obtained by considering inhibitory spike trains with IPSP distributions of different shapes, but with the same average amplitude $\langle|a_{i}|\rangle$. In particular, more asymmetric distributions, characterized by higher skewness and presenting longer tails towards larger IPSP amplitudes, induce lower firing rates, as shown in Fig. 2 (a) and (c) for supra-threshold and sub-threshold neurons, respectively. In the sub-threshold case this implies that passing from $\delta$-distributions, to uniform, to truncated Gaussian and to exponential ones, the firing onset will occur for larger and larger $\sigma^{2}$. For sub-threshold neurons the finite IPSP amplitude enhance the coherence resonance effect with respect to infinitesimal IPSP (corresponding to the DA), while the long inhibitory tails induce a shift of the minimum in the CV towards larger noise amplitudes (see Fig. 2 (b)). For supra-threshold neurons the increased asymmetry in the distributions simply induces more regular firing, as shown in Fig. 2 (d). It is quite peculiar that the TGD and the UD give almost identical results, despite the fact that the TGD is more asymmetric and characterized by a larger standard deviation, namely 0.61 mV for the TGD and 0.29 mV for the UD. Differences are instead seen with respect to the ED which reveals extremely long tails and a standard deviation of 1 mV and to the DD where there is no variability in the IPSP amplitude.

A much more detailed characterization of the firing time statistics, beyond the first two moments that we have considered so far, can be achieved by evaluating the spike train spectrum ${\hat{C}}(\omega)$ which is directly related with the first-passage time density (22). The comparison of the theoretical estimations with the numerical findings is reported in Fig. 3, showing a very good agreement for all the reported cases. We report each spectra normalized by the corresponding average firing rate $r_{0}$, in order to emphasize the changes produced by the shape of the IPSP distribution, rather than the changes due to the different values of the firing rates. From Fig. 3 (a), we observe that for $\delta$-distributed synaptic weights, the increase in the kick amplitude $a_{i}$ induces a higher peak in the spectrum at a lower frequency. Therefore, for increasing $a_{i}$ not only the rate decreases, as previously reported, but also the peak of the ISI distribution shifts from 41 msec for $|a_{i}|=0.1$ mV to 54 msec for $|a_{i}|=1$ mV, thus suggesting that the entire dynamics slows down for increasing IPSP amplitudes.

Furthermore, as shown in Fig. 3 (b) the shape of the distributions of the IPSPs has also an influence on the spectra. In particular, the increase in the asymmetry of the distributions induces peaks at lower and lower frequencies, while their height increases. The neuron activity is slowed down for longer and longer tails in the IPSP distributions.

These conclusions are supported by the data reported in Fig. 4 (a). In the figure are reported the average firing rate $r_{0}$ for TGDs with different standard deviation $\sigma_{G}$ and fixed position of the maximum at $a_{p}=-1$ mV. For increasing $\sigma_{G}$ the neuronal activity decreases for corresponding noise amplitudes. Since the value of the skewness increases with $\sigma_{G}$, more asymmetric IPSP distributions induce a lower neuronal firing rate. Furthermore, a larger asymmetry is also responsible for a shift of the position of the minimum of CV, associated to the CR phenomenon, towards larger $\sigma^{2}$ and for a more regular firing activity observable at $\sigma^{2}>\sigma_{m}^{2}$, as shown in Fig. 4 (b).

So far we considered as excitatory input a constant DC term, however the analytic approach here presented can be applied also for exponentially distributed excitatory amplitudes, namely $A_{e}=\exp(-a/a_{e})/a_{e}$. We have reported the results for this case for different inhibitory distributions in Fig. 5, the analytic estimations reproduce very well the numerical results both for the firing rate and the CV in an ample range of noise intensities. However, for small values of the noise intensity, the analytic results slightly deviates from the numerical values due to the fact that in the limit of small rates the numerical evaluation of the integrals entering in the expresssion of the generating function $Z_{0}$ have problems of convergece. The effect due to the different IPSP distributions is analogous to the one observed with a constant DC excitatory term. The DA, conversely, is unable to capture the precise values of these two quantities. Once more, for the reported case in Fig. 5 where ${\mu_{T}}=9$ mV, the CR effect is present. It is important to remark that, in the specific case of exponentially distributed amplitudes of the EPSP, the approach discussed in this article fails when an additional supra-threshold DC current is applied, namely for ${\mu_{0}}>v_{th}$. A brief discussion in this regard will be provided in the concluding section.

Heterogeneous Sparse Networks

The previous theoretical analysis of single neuron response to an external Poissonian input can find application also in the analysis of the dynamics of recurrent LIF heterogeneous networks with random sparse connectivity. For sparse networks, the spike-trains impinging a certain neuron can be assumed to be uncorrelated and Poissonian [29, 35]. Furthermore, similarly to what done in  [29], we can assume that the spike-trains in the networks can be self-consistently described as Poissonian processes with firing rates $r_{0}(\mu(j))$ related to the neuronal excitability $\mu(j)$ of each single neuron.

For the sake of simplicity, we will consider a network of inhibitory neurons, where each neuron is characterized by a different level of excitability $\{\mu(j)\}$ encompassing any excitatory external drive as well as the specific characteristic of the considered neuron. Therefore, the dynamics of the $j$-th neuron in the network can be written as

[TABLE]

where $K<<N$ is the number of synaptic neighbors, $C_{jk}$ is the connectivity matrix with entries 1 (0) if the $k$-th neuron is connected (not connected) to neuron $j$. The amplitudes of the IPSP $a_{i}(k)=g_{i}(k)/K<0$ associated to the firing of the $k$-th neuron is assumed to be randomly distributed following some of the distributions previously introduced.

For the population dynamics we will give an estimate of the average firing rate via a self-consistent approach by assuming that in average each neuron receive a single Poissonian spike train with a rate $R_{i}=\bar{r}_{0}K$ and where each IPSP has a random amplitude $a_{i}$ taken from a distribution $A_{i}(a)$. An estimation of the average firing rate can be obtained by solving the following implicit equation

[TABLE]

which is an extension to an heterogeneous network of Eq. (17) and where, as we have previously discussed, $F(c)$ depends on the excitatory input and $Z_{0}$ on the chosen IPSP distribution. It is important to stress that usually in inhibitory networks not all neurons are firing, but just a certain fraction $n^{*}$ will be active [34], therefore the integral reported in (57) is limited to these neurons, which are the one with higher excitability $\mu(j)\in\{\mu_{A}\}$ . Furthermore $\Delta=\int_{\{\mu_{A}\}}d\mu P(\mu)$ is the support of the active neurons. Once the firing rates have been obtained self consistently we can derive the coefficient of variation of each single neuron according to Eqs. (26) and (27) and then to perform the population average as follows

[TABLE]

In  [34] a theoretical approach to obtain self-consistently $n^{*}$, $\bar{r}_{0}$ and the average coefficient of variation $\overline{\text{CV}}$ has been developed for constant IPSP, i.e. for $\delta$-distribution. Here we extend such approach to more generic IPSP distributions, however we will limit to obtain the analytic estimations of $\bar{r}_{0}$ and $\overline{\text{CV}}$. The values of $n^{*}$, entering in the expressions for the average population rate and coefficient of variation, will be considered as parameter values obtained directly from the simulations. We have verified that the comparison with the numerical findings is very good for all the four considered IPSP distributions and for the whole considered range of average synaptic inputs $\langle a_{i}\rangle$. For illustration purposes we present in Fig. 6 (a) and (b) the average firing rate $\bar{r}_{0}$ and $\overline{\text{CV}}$, for the two distributions that have consistently presented larger differences between them, namely DD and ED, and for the DA. While it is evident that the DA overestimates $\bar{r}_{0}$ and $\overline{\text{CV}}$ already for $\langle|a_{i}|\rangle>0.5$ mV, we observe that in the case of heterogeneous networks the differences among the various IPSP distributions are quite limited at the level of the average firing rate $\bar{r}_{0}$ and $\overline{\text{CV}}$. In order to investigate more in details the influence of different IPSP distributions on the neuronal response, we have numerically estimated the probability distribution functions $P(r_{0})$ of the single neuron firing rate $r_{0}$ for DD and ED distributions. As shown in Fig. 6 (c) and (d), the $P(r_{0})$ obtained for the DD display a higher peak at low firing rate $r_{0}$ with respect to the ED, while for higher firing rates they essentially coincide. This difference should be ascribed to the fact that sub-threshold neurons subject to IPSP trains with ED have a firing onset at definitely larger noise amplitudes than the ones subject to IPSPs with DD, while for supra-threshold neurons the differences are definitely less evident, as previously reported in Fig. 2 (a) and (c). This explains also why the exam of the average firing rate does not reveal large difference between the two distributions, since the neurons with low firing rate contribute negligibly to the average activity. The influence of synaptic weight distributions on neuronal population dynamics has been usually examined in the context of homogeneous neuronal populations [17, 13], where the single neuron response represents a good mean field approximation of the network dynamics. However, from the present analysis it emerges that the heterogeneity in the single neuron excitability can render extremely difficult to distinguish among stimulations with different IPSP distributions.

As a final remark, we would like to stress that the reported approach works very well also for heterogeneous sparse networks, provided that the collective dynamics is asynchronous. Otherwise, the presence of partial synchronization or of collective oscillations can induce correlations in the input spike trains, which cannot be accounted for with this approach. In particular, to avoid phase locking among the neurons the distribution $P(\mu)$ of the excitabilities should be sufficiently wide.

Conclusions

In this paper we have reported a theoretical methodology to obtain exact firing statistics for leaky integrate-and-fire neurons subject to discrete inhibitory noise, accounting for Poissonian trains of uncorrelated post-synaptic potentials. Our results represent an extension to generic synaptic weights distribution of the approach developed in [17] for exponentially distributed post-synaptic potentials. In particular, we report explicit results for the firing rate, the coefficient of variation and the spike train spectrum.

The comparison with numerical simulations reveals a very good agreement for all the considered distributions over all the reported ranges of shot noise amplitude. Moreover, the method is also able to reproduce the average activity of an heterogeneous inhibitory neural network with sparse connectivity, by making use of a self-consistent mean field formulation. Conversely, the diffusion approximation [3] (the most used theoretical approach), gives a reasonable estimate of the firing time statistics for sufficiently small IPSP amplitudes, but the agreement rapidly degrades, and reveals large discrepancies already for amplitudes $>0.5$ mV, corresponding to physiologically relevant values [8, 36, 37].

As a general result we observe that the firing statistics of single neurons is strongly influenced by the shape of the IPSP distributions. Distributions with longer tails lead to smaller firing rate and in general to more regular spike trains (for sufficiently large inhibition). For heterogeneous networks the value of the average IPSP has a noticeable influence on the firing activity of the neuronal population, while the shape of the synaptic weight distributions seem to have a really limited impact on the average properties of the network. However, differences induced by the IPSP distributions are still observable at the level of single neuron.

We have shown that the method works for any choice of IPSP amplitude distributions, however one has to be careful when dealing with the excitatory drive. We have shown that the agreement with numerical simulations is very good when the firing of the neuron is either promoted exclusively by a suprathreshold DC current or provided by the stochastic arrivals of exponentially distributed EPSP amplitudes. In this latter case, an excitatory DC current can as well be present in addition to the EPSP trains, however the DC contribution must be strictly sub-threshold, namely, $\mu_{0}<v_{th}$. Whenever the two firing mechanisms are active at the same time (i.e. $\mu_{0}>v_{th}$), the formulation reported in this paper is no more valid due to the inconsistency in the choice of the proper integration limit $\bar{x}$ in Eq. (17). This limitation is illustrated in Fig. 7, where it is reported the response of a LIF neuron, subject to a supra-threshold excitatory DC current $\mu_{0}=12$ mV as well as to Poissonian spike trains of ED excitatory and DD inhibitory synaptic inputs. For small values of the noise intensity, the theoretical approach dramatically fails. In such a region the activity is almost exclusively current driven; i.e, the firing of the neuron promoted by the supra-threshold DC current is much faster than the arrival rate $R_{e}$ of the EPSPs. This can be confirmed by the fact that in this regime the firing rate coincides with $r_{0}^{t}$ in Eq. (51) for a tonic firing LIF neuron subject to a constant DC current $\mu_{T}=\mu_{0}$. As the noise intensity grows, corresponding to an increased rate $R_{e}$ of arrival of the EPSPs, the numerical data approach the theoretical prediction obtained for exponentially distributed EPSPs. The crossover occurs for $R_{e}\approx 2r_{0}^{t}$, indicating that the firing activity of the neuron is now mainly driven by the stochastic component.

An interesting semi-analytic approach has been recently reported for excitatory shot noise in [13]. However, further progresses are required to achieve exact firing time statistics going beyond the diffusion approximation for generic distributions of instantaneous excitatory PSPs, as well as for more realistic neuronal models like the Exponential Integrate and Fire (EIF) neuron [1] able to reproduce quite accurately the dynamics of cortical neurons [2]. In particular, as shown in the Supplementary Information, the LIF neuron response is a limit case of the EIF dynamics both within the diffusion approximation as well as for shot noise. Furthermore, the diffusion approximation fails also for the EIF to capture the neuronal firing statistics for sufficiently large IPSP amplitudes in particular at low noise intensities.

Acknowledgments

We thank for useful discussions B. Lindner, G. Mato, A. Politi, G. Martelloni, MJE Richardson, M. Timme. This work has been partially supported by the European Commission under the program “Marie Curie Network for Initial Training”, through the project N. 289146, “Neural Engineering Transformative Technologies (NETT)” and by the A∗MIDEX grant (No. ANR-11-IDEX-0001-02) funded by the French Government “Programme Investissements d’Avenir” (D.A.-G. and A.T.). AT gratefully acknowledges a “VELUX Visiting Professor Grant” from Villum Fonden that supported his stay at Aarhus University during the initial stages of this work in 2011. AI is supported by the Danish Council for Independent Research and the Villum Fonden. This work has been completed at Max Planck Institute for the Physics of Complex Systems in Dresden (Germany) as part of the activity of the Advanced Study Group 2016 entitled ”From Microscopic to Collective Dynamics in Neural Circuits”.

Author contributions statement

S.O and D.A-G performed the numerical simulations. D.A-G and A.T prepared the manuscript. All the authors developed the theoretical methods and reviewed the manuscript.

Additional information

Competing financial interests Authors declare to have no competing financial interests that might have influenced this manuscript.

Supplementary Information

Exponential Integrate and Fire

The Exponential Integrate-and-Fire (EIF) model is a simple non-linear integrate and fire neuronal model introduced by Fourcaud-Trocmé et al. [1] able to reproduce quite accurately the dynamics of cortical neurons [2]. The model can be written as follows

[TABLE]

where $\mu_{0}$ represents an external DC current and $I(t)$ the synaptic drive. In this model, the spike generation occurs in a finite time controlled by the parameter $\Delta_{T}$. In particular, once the membrane potential has reached the threshold value $v_{th}$ this will rapidly grow towards infinity in a finite time interval, the parameter $\Delta_{T}$ establishes how fast the infinite limit is reached. In the limit ${\Delta_{T}\rightarrow 0}$, the spike generation is instantaneous and the LIF model is recovered. As in the usual LIF model, once the neuron has fired its membrane potential is resetted to the value $v_{re}=5$ mV, we also set $\tau=20$ ms and $v_{th}=10$ mV, as in the LIF model studied in the article.

As a first analysis, we will examine the response of the EIF neuron subject to small Gaussian noise of zero average and intensity $\sigma$, this can be obtained by solving the associated continuity equation for the probability $P(v,t)$ of finding the membrane voltage between $v$ and $v+dv$ at time $t$, which reads as [1]:

[TABLE]

where $J=J(v,t)$ is the associated flux. . In this case, since the effective threshold is located at infinity, the steady firing rate can be evaluated as

[TABLE]

where $J(v)$ is the stationary solution of the continuity equation for the flux.

In particular, we made use of the threshold-integration method [3] to calculate the firing rate $r_{0}$ of the EIF neuron subject to inhibitory inputs and compare it with the diffusion approximation (DA) for the LIF and the corresponding shot noise solution for $\delta$-distributed IPSP amplitudes with $|a_{i}|=0.1$ mV. The results are shown in Fig. 8 (a), where it is clearly shown that for $\Delta_{T}\to 0$ the EIF results converge to the LIF solution, both for the DA and the the shot noise results found for small IPSP. Furthermore direct simulations of the EIF and the corresponding shot noise solution of the LIF for large IPSP amplitudes, namely $|a_{i}|=1$ mV, show that also in this case the LIF limit is recovered for $\Delta_{T}\to 0$. However it is clear that also in the case of the EIF, the diffusion limit is unable to capture the onset of the activity and the firing rate at small intensities.