Introduction

Neurons in the neocortex in vivo are subject to a continuous synaptic bombardment reflecting the intense network activity1. In the so-called high-input regime, in which neurons receive hundreds of synaptic inputs during each interspike interval2, 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 irrelevant5. However several experimental studies have revealed that rare PSPs of large amplitude can have a fundamental impact in the network activity6, 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 amplitudes6,7,8,9,10,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 distributions13 and in balanced sparse networks for conductance based LIF neurons with log-normal EPSP distributions14. The presence of few strong synapses induce faster and more reliable responses of the network even for small inputs13, 14. Interestingly, in ref. 14 it has been shown that a single neuron driven by random synaptic inputs log-normally distributed reveals a clear aperiodic stochastic resonance15, 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 ref. 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 weights17. Furthermore, results for generic EPSP distributions have been obtained in ref. 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 neuron18, 19 subject to a constant external DC current μ 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

$$\frac{dv}{dt}=-\,\frac{v-{\mu }_{0}}{\tau }+I(t\mathrm{).}$$
(1)

where τ = 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

$$I(t)=\sum _{\{{t}_{e}\}}{a}_{e}\delta (t-{t}_{e})+\sum _{\{{t}_{i}\}}{a}_{i}\delta (t-{t}_{i}),$$
(2)

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 model20 with pulse amplitudes randomly drawn from distributions A x (a). This 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 neurons21, 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) \(\mathop{C}\limits^{\frown {}}(\omega )\) 23. To obtain such results we follow the approach developed in ref. 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:

$$\frac{\partial P}{\partial t}+\frac{\partial J}{\partial v}=r(t)[\delta (v-{v}_{re})-\delta (v-{v}_{th})]+\delta (t)\delta (v-{v}_{re}\mathrm{).}$$
(3)

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 = 024.

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

$$J=-\,(\frac{v-{\mu }_{0}}{\tau })P+{J}_{e}+{J}_{i};$$
(4)

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

$${J}_{x}={R}_{x}{\int }_{-\infty }^{v}dwP(w,t){\int }_{v-w}^{\infty }da{A}_{x}(a\mathrm{).}$$
(5)

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

$${J}_{e}({v}_{th},t)=r(t),\,{J}_{i}({v}_{th},t)=0\,,\,P({v}_{th},t)=0;$$
(6)

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

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

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 solve25, 26. The reason is related to the overshoots over the threshold v th induced by the finite amplitude of the PSPs, which render 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 distributions27. Richardson and Swarbrick17 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}={\rm{\Theta }}(a)\exp (-a/{a}_{e})/{a}_{e}\), and to constant excitatory synaptic drive, encompassed in an external DC current μ 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

$$\frac{\partial {J}_{e}}{\partial v}+\frac{{J}_{e}}{{a}_{e}}={R}_{e}P(v,t)-r(t)\delta (v-{v}_{th}),$$
(7)

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

$${\tilde{J}}_{e}(s,t)=\tilde{P}(s,t){Q}_{e}(s)-r(t){S}_{e}(s\mathrm{).}$$
(8)

For the particular case of the exponentially distributed EPSP amplitudes

$${Q}_{e}(s)=\frac{{R}_{e}}{1-s{a}_{e}}\,{S}_{e}(s)=\frac{{{\rm{e}}}^{s{v}_{th}}}{1-s{a}_{e}}\mathrm{.}$$
(9)

Whenever the excitatory input is simply given by the DC current μ 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 Z0 ≡ Z0(s), which corresponds to the generating function for the sub-threshold voltage moments. Therefore, Z0 can be estimated as the Laplace transform of P(v, t) when v th  → ∞, 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

$$\frac{d{{\rm{Z}}}_{0}}{ds}=\tau (\frac{{\mu }_{0}}{\tau }{{\rm{Z}}}_{0}+{\tilde{J}}_{e}+{\tilde{J}}_{i});$$
(10)

where we set \({{\rm{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:

$$\frac{1}{{{\rm{Z}}}_{0}}\frac{d{{\rm{Z}}}_{0}}{ds}=\tau (\frac{{\mu }_{0}}{\tau }+{Q}_{e}+{Q}_{i}),$$
(11)
$${{\rm{Z}}}_{0}={H}_{e}\exp ({\mu }_{0}s+\tau {\int }_{0}^{s}ds{Q}_{i}),$$
(12)

where the excitatory contribution is encompassed in the term

$${H}_{e}=(\begin{matrix}{(1-s{a}_{e})}^{-\tau {R}_{e}} & {\rm{for}}\,{A}_{e}(a)={\rm{\Theta }}(a)\exp (\,-\,a/{a}_{e})/{a}_{e}\\ 1 & {\rm{for}}\,{\rm{DC}}\,{\rm{current}}\,\end{matrix}$$
(13)

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

$$\frac{d{\tilde{P}}_{0}}{ds}={\tilde{P}}_{0}\tau (\frac{{\mu }_{0}}{\tau }+{Q}_{e}+{Q}_{i})-\frac{{r}_{0}}{s}({e}^{s{v}_{re}}-{e}^{s{v}_{th}}+s{S}_{e}).$$
(14)

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

$$\frac{1}{{{\rm{Z}}}_{0}}\frac{d\tilde{P}}{ds}=\frac{\tilde{P}}{{{\rm{Z}}}_{0}^{2}}\frac{d{{\rm{Z}}}_{0}}{ds}-\tau \frac{{r}_{0}}{s{{\rm{Z}}}_{0}}F(s).$$
(15)

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

$$\frac{\tilde{P}(s)}{{{\rm{Z}}}_{0}(s)}=\tau {r}_{0}{\int }_{s}^{\bar{x}}\frac{dc}{c{{\rm{Z}}}_{0}(c)}F(c)+\frac{\tilde{P}(\bar{x})}{{{\rm{Z}}}_{0}(\bar{x})}.$$
(16)

Notice that, although we have derived an expression for Z0, 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/{{\rm{Z}}}_{0}(\bar{x})\) vanishes, the following exact analytic expression for the stationary firing rate can be obtained

$$\frac{1}{\tau {r}_{0}}={\int }_{0}^{\bar{x}}\frac{dc}{c{{\rm{Z}}}_{0}(c)}F(c),$$
(17)

where we have made use of the normalization condition of the probability densities, i.e. \({\tilde{P}}_{0}(0)={{\rm{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 \({{\rm{Z}}}_{0}^{-1}(\bar{x})=0\), can be easily found to be

$$\bar{x}=(\begin{matrix}1/{a}_{e} & {\rm{for}}\,{H}_{e}={(1-s{a}_{e})}^{-\tau {R}_{e}}\\ \infty & {\rm{for}}\,{H}_{e}=1.\end{matrix}$$
(18)

It should be remarked that in presence of both sources of excitatory drive, the integration limit can still be identified whenever μ 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 Eqs (3) and (4). Namely, we obtain

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

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 Z0 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

$$\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)};$$
(20)

where

$$A(s)=\frac{{e}^{s{v}_{re}}}{{Z}_{0}};\,B(s)=\frac{F(s)}{{Z}_{0}}\mathrm{.}$$
(21)

Equation (20) diverges exactly at ω ≡ 0, and in that case it should be complemented with the expression \(\hat{\rho }(\omega =\mathrm{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 \(\mathop{C}\limits^{\frown {}}(\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{\rm{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:

$$\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 ),$$
(22)

where \({\mathrm{lim}}_{\omega \to 0}\hat{C}(\omega )={r}_{0}\) and \({\mathrm{lim}}_{\omega \to \infty }\hat{C}(\omega )={{\rm{CV}}}^{2}{r}_{0}\) 28.

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

$${\frac{{\partial }^{n}\hat{q}}{\partial {\omega }^{n}}|}_{\omega =0}={(-{\rm{i}})}^{n}\langle {t}^{n}\rangle \mathrm{.}$$
(23)

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

$$\hat{\rho }(\omega )\approx \frac{{n}_{0}+{n}_{1}\omega +{n}_{2}{\omega }^{2}}{{d}_{0}+{d}_{1}\omega +{d}_{2}{\omega }^{2}};$$
(24)

where n 0 = −1, d 1 = i/r 0, d 0 = 0 and

$${n}_{1}=i{\int }_{0}^{\bar{x}}ds\,\mathrm{log}\,s\frac{d}{ds}A(s)\,{d}_{2}={\int }_{0}^{\bar{x}}ds\frac{\mathrm{log}\,s}{s}B(s\mathrm{).}$$
(25)

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 ω = 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

$$\langle t\rangle =-\,i{d}_{1}=\mathrm{1/}{r}_{0}\,\langle {t}^{2}\rangle =2\frac{{d}_{1}^{2}-{d}_{2}{n}_{0}+{d}_{1}{n}_{1}}{{n}_{0}^{2}};$$
(26)

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

$${\rm{C}}{\rm{V}}=\frac{\sqrt{\langle{t}^{2}\rangle-{\langle t\rangle}^{2}}}{\langle t\rangle}.$$
(27)

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

Explicit expressions of Z0 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), δ (DD), uniform (UD) and truncated Gaussian distribution (TGD). The exact expressions for Z0(s) for each of these four cases are reported in the following, together with the average, variance and skewness of the corresponding distributions A i (a).

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

$${A}_{i}(a)=-\,\frac{{\rm{\Theta }}(\,-\,a)}{{a}_{i}}{{\rm{e}}}^{-a/{a}_{i}}\mathrm{.}$$
(28)

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

$$\frac{\partial {J}_{i}}{\partial v}+\frac{{J}_{i}}{{a}_{i}}={R}_{i}P(v,t).$$
(29)

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

$${\tilde{J}}_{i}=\frac{{R}_{i}}{1-s{a}_{i}}\tilde{P}.$$
(30)

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

$${{\rm{Z}}}_{0}(s)={H}_{e}\frac{{{\rm{e}}}^{{\mu }_{0}s}}{{(1-s{a}_{i})}^{{R}_{i}\tau }}\mathrm{.}$$
(31)

The mean and variance due to the distribution of the inhibitory synaptic weights for the exponential case are 〈a i 〉 = a i and \(var[{a}_{i}]={a}_{i}^{2}\); while the absolute value of the skewness is two.

δ Distribution (DD): In the case in which the neuron receives spikes with a constant amplitude a i , the distribution becomes simply a δ-function:

$${A}_{i}(a)=\delta (a-{a}_{i})\mathrm{.}$$
(32)

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

$$\frac{\partial {J}_{i}}{\partial v}={R}_{i}[P(v,t)-P(v-{a}_{i},t)]\mathrm{.}$$
(33)

and the associate Laplace transform reads as

$${\tilde{J}}_{i}=-\,{R}_{i}\tilde{P}\frac{(1-{{\rm{e}}}^{s{a}_{i}})}{s}\mathrm{.}$$
(34)

For this simple distribution the generating function Z0 reads as:

$${{\rm{Z}}}_{0}(s)={H}_{e}\frac{\exp ({\mu }_{0}s+\tau {R}_{i}{E}_{I}({a}_{i}s))}{{C}_{0}{s}^{\tau {R}_{i}}};$$
(35)

where \({E}_{I}(y)={\int }_{y}^{\infty }dy{e}^{y}/y\) is the exponential integral and C 0 accounts for the normalization condition requiring that Z0(0) = 1 and its explicit expression is

$${C}_{0}=\exp (\tau {R}_{i}({\rm{\Gamma }}+\,\mathrm{log}|{a}_{i}|)),$$
(36)

where Γ ≈ 0.577731 is the Euler-Mascheroni constant.

For this distribution the mean and variance of the inhibitory process are simply 〈a i 〉 = 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

$${A}_{i}(a)=\frac{{\rm{\Theta }}(a-{l}_{1})-{\rm{\Theta }}(a-{l}_{2})}{{l}_{2}-{l}_{1}}\mathrm{.}$$
(37)

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

$$\frac{\partial {J}_{i}}{\partial v}=\frac{{R}_{i}}{{l}_{2}-{l}_{1}}{\int }_{-\infty }^{v}dwP(w,t)[(v-w-{l}_{2}){\rm{\Theta }}(v-w-{l}_{2})-(v-w-{l}_{1}){\rm{\Theta }}(v-w-{l}_{1})]$$
(38)

and the corresponding Laplace transform reads as

$${\tilde{J}}_{i}=-\,\tilde{P}\frac{{R}_{i}}{{l}_{2}-{l}_{1}}[\frac{{{\rm{e}}}^{{l}_{1}s}}{{s}^{2}}-\frac{{{\rm{e}}}^{{l}_{2}s}}{{s}^{2}}+\frac{{l}_{2}-{l}_{1}}{s}]\mathrm{.}$$
(39)

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

$${{\rm{Z}}}_{0}(s)={H}_{e}\frac{\exp ({\mu }_{0}s+\frac{\tau {R}_{i}}{{l}_{2}-{l}_{1}}({l}_{2}{E}_{I}({l}_{2}s)-{l}_{1}{E}_{I}({l}_{1}s)+\frac{{{\rm{e}}}^{{l}_{1}s}}{s}-\frac{{{\rm{e}}}^{{l}_{2}s}}{s}))}{{C}_{0}{s}^{\tau {R}_{i}}}$$
(40)
$${C}_{0}=\exp (\frac{\tau {R}_{i}}{{l}_{2}-{l}_{1}}({l}_{1}(1-{\rm{\Gamma }}-\,\mathrm{log}|{l}_{1}|)-{l}_{2}(1-{\rm{\Gamma }}-\,\mathrm{log}|{l}_{2}|)))$$
(41)

The mean and variance of the inhibitory shot noise are 〈a i 〉 = (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 distribution38, which is peaked at a p and with a standard deviation equal to σ G . For notation simplicity we write the Gaussian distribution as

$$\varphi (y)=\frac{1}{\sqrt{2\pi }{\sigma }_{G}}\exp (-\frac{{(y-{a}_{p})}^{2}}{2{\sigma }_{G}^{2}})\mathrm{.}$$
(42)

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

$${A}_{i}(a)=\frac{{\rm{\Theta }}(\,-\,a)\varphi (a)}{{\rm{\Phi }}\mathrm{(0)}};$$
(43)

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

$${\rm{\Phi }}(y)=\frac{1}{2}(1+{\rm{erf}}(\frac{y-{a}_{p}}{{\sigma }_{G}\sqrt{2}}))\mathrm{.}$$
(44)

The flux of inhibitory probability takes the form

$$\frac{\partial {J}_{i}}{\partial v}=\frac{{R}_{i}}{{\rm{\Phi }}\mathrm{(0)}}({\rm{\Phi }}\mathrm{(0)}P(v,t)+{\int }_{-\infty }^{v}dwP(w,t)\varphi (v-w){\rm{\Theta }}(v-w))\mathrm{.}$$
(45)

The corresponding Laplace transform is:

$${\tilde{J}}_{i}=-\,\tilde{P}\frac{{R}_{i}}{{\rm{\Phi }}\mathrm{(0)}s}({\rm{\Phi }}\mathrm{(0)}-\frac{1}{2}{\rm{erfc}}(\frac{s{\sigma }_{G}^{2}+{a}_{p}}{\sqrt{2}{\sigma }_{G}}){{\rm{e}}}^{s(s{\sigma }_{G}^{2}+2{a}_{p})\mathrm{/2}})$$
(46)

The complicated expression appearing in Eq. (46) does not allow us to obtain the explicit expression for Z0. However it is easy to integrate numerically Eq. (12) once we know \({\tilde{J}}_{i}\). When dealing with a width of the Gaussian distribution σ 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: \({\rm{erfc}}(y)\approx \exp (\,-\,{y}^{2})/\sqrt{\pi }y\) when y 1 and Eq. (46) can be simplified in such cases to:

$${\tilde{J}}_{i}\approx -\,\tilde{P}\frac{{R}_{i}}{{\rm{\Phi }}\mathrm{(0)}s}({\rm{\Phi }}\mathrm{(0)}+{\sigma }_{G}\frac{{{\rm{e}}}^{-\frac{{a}_{p}^{2}}{2{\sigma }_{G}^{2}}}}{\sqrt{2\pi }{\sigma }_{G}^{2}s+{a}_{p}})\,{\rm{if}}\,s,{\sigma }_{G} > 1$$
(47)

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

$$\langle {a}_{i}\rangle ={a}_{p}-{\sigma }_{G}^{2}\frac{\varphi \mathrm{(0)}}{{\rm{\Phi }}\mathrm{(0)}}$$
(48)
$$var[a]={\sigma }_{G}^{2}(1+{a}_{p}\frac{\varphi \mathrm{(0)}}{{\rm{\Phi }}\mathrm{(0)}}+{({\sigma }_{G}\frac{\varphi \mathrm{(0)}}{{\rm{\Phi }}\mathrm{(0)}})}^{2});$$
(49)

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

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 μ T and noise intensity σ2.

For a neuron receiving an inhibitory Poissonian spike train at a rate R i and with an average synaptic weight 〈a i 〉 the effective average input is

$${\mu }_{T}={\mu }_{e}+{\mu }_{i}={\mu }_{e}+\tau {R}_{i}\langle {a}_{i}\rangle \mathrm{.}$$
(50)

where a i  < 0 and μ e is the average excitatory input. When only the drift term is present μ e  = μ 0, while for a neuron receiving also an excitatory Possonian spike train of rate R e and with average synaptic weight 〈a e 〉, it becomes μ e  = μ 0 + τR e a e 〉. On the one hand, for an effective sub-threshold input μ 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 μ T and the activation time associated to the escape process from the resting state to the threshold induced by the fluctuations in the input28. On the other hand, for a supra-threshold LIF neuron for which μ T  > v th , in absence of a refractory state, the only characteristic time is the tonic firing rate

$$\frac{1}{{r}_{0}^{t}}=\tau {\rm{l}}{\rm{o}}{\rm{g}}(\frac{{\mu }_{T}-{v}_{re}}{{\mu }_{T}-{v}_{th}}).$$
(51)

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 σ 2 associated to these two uncorrelated processes is given by the sum of the variances of each process, namely

$${\sigma }^{2}=\sum _{x=i,e}{R}_{x}\tau ({\langle {a}_{x}\rangle }^{2}+var[{a}_{x}])\mathrm{.}$$
(52)

From Eqs (50) and (52) one can see that the values of μ T and σ can be independently tuned with an appropriate selection of the parameter μ e and R i for any fixed average inhibitory amplitude. Whenever the excitatory drive is given simply by a DC term μ 0, for any choice of the couple of values {μ T , σ} one has an unique solution, which however does not always correspond to a firing activity. For instance, at small values of σ 2 one can find values of μ 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 μ e and therefore a whole family of solutions exist for any fixed a i . Also, 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 approximation3. 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 μ T and σ 2 through the formulas

$$\frac{1}{{r}_{0}}=\tau \sqrt{\pi }{\int }_{{y}_{r}}^{{y}_{\theta }}dx{{\rm{e}}}^{{x}^{2}}(1+{\rm{erf}}(x))$$
(53)
$${{\rm{CV}}}^{2}=\frac{2\pi {\tau }^{2}}{{r}_{0}^{2}}{\int }_{{y}_{r}}^{{y}_{\theta }}dx{{\rm{e}}}^{{x}^{2}}{\int }_{-\infty }^{x}dy{{\rm{e}}}^{{y}^{2}}{\mathrm{(1}+{\rm{erf}}(y))}^{2}$$
(54)

where y θ  = (v th  − μ T )/σ and y r  = (v re  − μ T )/σ.

By following17, Eqs (53) and (54) can be obtained from the shot noise formulation, namely from Eqs (26) and (22) by considering the expansion of Z0 limited to the first two voltage moments:

$${{\rm{Z}}}_{0}(s)\simeq \exp (s{\mu }_{T}+\frac{{s}^{2}}{4}{\sigma }^{2})\mathrm{.}$$
(55)

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 DA3, 4, 29. 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 distributions13, 17. 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 networks29, 30. 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 generated31.

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 μ e  = μ 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 σ 2 at a constant effective input current, namely μ 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 σ 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 systems32. 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 times28.

Figure 1
figure 1

Firing time statistics for different IPSP amplitudes with δ-distributions (DD). Average firing rate r 0 (a) and coefficient of variation CV (b) as a function of noise intensity σ 2. The red curves/symbols correspond to |a i | = 1 mV and the blue ones to |a i | = 0.1 mV. Firing rate r 0 (c) and CV (d) as a function of the synaptic amplitude |a i | for constant noise intensity, namely σ 2 = 2 mV2. Black dashed curves refer to the diffusion approximation (DA) (Eqs (53) and (54)); solid curves to the theoretical results for shot noise with constant amplitude a i (Eqs (17) and (27)) and symbols to the corresponding numerical simulations. Simulations were performed by exactly integrating Eq. (1) with an event driven scheme (see ref. 34 for details) and the statistics were estimated over ≈107 spikes. For all the panels the effective input is μ T  = 9.0 mV and the excitatory drive is simply a DC term; other parameter values are v re  = 5 mV, v th  = 10 mV, τ = 20 ms.

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 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 μ 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 μ T are longer for larger IPSP and this induces a deeper minimum in the CV as reported in ref. 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 μ T while increasing |a i |. These two quantities can be independently tuned with a suitable selection of R i and μ 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 μ 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 σ 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 σ 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 〈|a i |〉. 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,c) for supra-threshold and sub-threshold neurons, respectively. In the sub-threshold case this implies that passing from δ-distributions, to uniform, to truncated Gaussian and to exponential ones, the firing onset will occur for larger and larger σ 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.

Figure 2
figure 2

Firing time statistics for different IPSP distributions. Firing rate r 0 and CV as a function of noise intensity for sub-threshold (μ T  = 9 mV) (a,b) and supra-threshold neurons (μ T  = 11 mV) (c,d). In all the panels the symbols correspond to numerical simulations, the dashed black lines to the DA, calculated from Eqs (53) and (54), and the solid lines to the theoretical results reported in Eq. (17) and Eq. (27). The average IPSP amplitude is set to 〈|a i |〉 = 1 mV for all the distributions. For the TGD: the peak position |a p | and the width σ G of the distribution are equal (namely, ≈0.7766 mV). For the UD: l 1 = 2a i and l 2 = 0. Other parameters and simulation procedures as in Fig. 1.

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 δ-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.

Figure 3
figure 3

Spike Train Spectra. (a) Normalized spike train spectra \(\hat{C}(\omega )/{r}_{0}\) for the same cases shown in Fig. 1(a,b) for an effective input μ T  = 9 mV and noise intensity σ 2 = 2 mV2; (b) Normalized spike train spectra for the distributions considered in Fig. 2(a,b) for 〈|a i |〉 = 1, μ T  = 9 mV and σ 2 = 4 mV2. In both panels, the symbols correspond to the simulation results, while the continuous line to the theoretical estimations using Eq. (22). Simulated spectra were obtained by calculating the squared modulus of the Fourier transform of the spike train using a time trace of 64 s with 1 ms binning, and averaging over 10,000 realizations.

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 σ G and fixed position of the maximum at a p  = −1 mV. For increasing σ G the neuronal activity decreases for corresponding noise amplitudes. Since the value of the skewness increases with σ 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 σ 2 and for a more regular firing activity observable at \({\sigma }^{2} > {\sigma }_{m}^{2}\), as shown in Fig. 4(b).

Figure 4
figure 4

Effect of the distribution asymmetry on the firing time statistics. (a) Average firing rate r 0 as a function of the noise intensity for different values of the standard deviation of the Gaussian distribution σ G as indicated in the legend. (b) Coefficient of variation CV for the same cases depicted in (a). The reported data refer to TGDs with the peak located at |a p | = 1 mV and to an effective input μ T  = 9 mV. Symbols correspond to numerical simulations, and the solid lines to the theoretical results reported in Eq. (17) and Eq. (27). Other parameters and simulation procedures as in Fig. 1

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 Z0 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 μ 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 μ 0 > v th . A brief discussion in this regard will be provided in the concluding section.

Figure 5
figure 5

Firing time statistics for exponentially distributed EPSPs. Firing rate r 0 (a) and CV (b) as a function of noise intensity for excitatory synaptic weights distributed exponentially with average amplitude 〈a e 〉 = 〈|a i |〉 = 1 mV. In this figure excitatory and inhibitory spike trains are balanced, i.e R e  = R i , while the effective input is sub-threshold, namely μ T  = μ 0 = 9 mV. The symbols and lines have the same definition as in Fig. 2, as well as all the other parameters for the IPSP distributions and the numerical simulations.

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 Poissonian29, 35. Furthermore, similarly to what done in ref. 29, we can assume that the spike-trains in the networks can be self-consistently described as Poissonian processes with firing rates r 0(μ(j)) related to the neuronal excitability μ(j) of each single neuron j.

For the sake of simplicity, we will consider a network of inhibitory neurons, where each neuron is characterized by a different level of excitability {μ(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

$$\dot{v}(j)=-\,\frac{[v(j)-\mu (j)]}{\tau }+\frac{1}{K}\sum _{k}{C}_{jk}{g}_{i}(k)\delta (t-{t}_{k});$$
(56)

where KN 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

$$\tau {\bar{r}}_{0}=\frac{1}{{\rm{\Delta }}}{\int }_{\{{\mu }_{A}\}}d\mu P(\mu ){[{\int }_{0}^{\bar{x}}\frac{dc}{c}\frac{F(c)}{c{Z}_{0}(c,\mu ,{\bar{r}}_{0})}]}^{-1};$$
(57)

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 Z0 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 active34, therefore the integral reported in (57) is limited to these neurons, which are the one with higher excitability μ(j)  {μ A }. Furthermore \({\rm{\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

$$\overline{{\rm{CV}}}=\frac{1}{{\rm{\Delta }}}{\int }_{\{{\mu }_{A}\}}d\mu P(\mu ){\rm{CV}}(\mu )\mathrm{.}$$
(58)

In ref. 34 a theoretical approach to obtain self-consistently n *, \({\bar{r}}_{0}\) and the average coefficient of variation \(\overline{{\rm{CV}}}\) has been developed for constant IPSP, i.e. for δ-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{{\rm{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 〈a i 〉. For illustration purposes we present in Fig. 6(a,b) the average firing rate \({\bar{r}}_{0}\) and \(\overline{{\rm{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{{\rm{CV}}}\) already for 〈|a i |〉 > 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{{\rm{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,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,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 populations13, 17, 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.

Figure 6
figure 6

Firing time statistics for heterogenoeus networks as a function of the average synaptic weight. (a) Average network frequency \({\bar{r}}_{0}\) and (b) average coefficient of variation \(\overline{{\rm{CV}}}\) as a function of the average synaptic weight 〈|a i |〉 for two representative IPSP distributions: DD (blue) and ED (green). Symbols correspond to simulation values and solid lines to the theoretical results from Eqs (57) and (58). Dashed line denote the results of the DA. The numerically estimated probability distribution functions P(r 0) of the single neuron firing rate r 0 are also reported for the two considered IPSP distributions for two values of the average synaptic weight: namely 〈|a i |〉 = 0.04 (c) and 〈|a i |〉 = 0.4 (d). Simulations were performed with N = 400 neurons, K = 20 and a distribution of the input currents uniformly chosen in the interval μ(j)  [10, 11] mV. The time averages are calculated, after discarding an initial transient corresponding to 106 spikes, over the following 106 spikes. Silent neurons (those that do not emit any spike in the considered time lapse) are not included in the statistics. In all cases 〈|a i |〉 denotes the average value of the corresponding distribution.

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(μ) 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 ref. 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 approximation3 (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 values8, 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, μ 0 < v th . Whenever the two firing mechanisms are active at the same time (i.e. μ 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 μ 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 μ T  = μ 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 2{r}_{0}^{t}\), indicating that the firing activity of the neuron is now mainly driven by the stochastic component.

Figure 7
figure 7

Firing time statistics for exponentially distributed EPSPs and supra-threshold DC current. Firing rate r 0 as a function of noise intensity for excitatory synaptic weights distributed exponentially and a supra-threshold DC current. All parameters as in Fig. 5 except for μ T  = μ 0 = 12 mV. The dashed red horizontal line refers to the rate \({r}_{0}^{t}\) of the LIF neuron subject to constant current μ T (Eq. (51)). Other symbols and lines have the same meaning as in Fig. 2, as well as all the other parameters for the IPSP distribution and the numerical simulations.

An interesting semi-analytic approach has been recently reported for excitatory shot noise in ref. 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) neuron39 able to reproduce quite accurately the dynamics of cortical neurons40. 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.