Equivalence between synaptic current dynamics and heterogeneous propagation delays in spiking neuron networks

TL;DR

This paper demonstrates a formal equivalence between synaptic current dynamics and propagation delays in spiking neuron networks, providing a unified framework to understand their impact on collective neural behavior.

Contribution

It introduces a mean-field reduction that shows local synaptic filtering and transmission delays are formally equivalent in affecting network dynamics.

Findings

Derived an effective Markovian model for population dynamics.

Proved equivalence between synaptic filtering and delays across various conditions.

Applicable to finite neuron networks and different firing regimes.

Abstract

Message passing between components of a distributed physical system is non-instantaneous and contributes to determine the time scales of the emerging collective dynamics like an effective inertia. In biological neuron networks this inertia is due in part to local synaptic filtering of exchanged spikes, and in part to the distribution of the axonal transmission delays. How differently these two kinds of inertia affect the network dynamics is an open issue not yet addressed due to the difficulties in dealing with the non-Markovian nature of synaptic transmission. Here, we develop a mean-field dimensional reduction yielding to an effective Markovian dynamics of the population density of the neuronal membrane potential, valid under the hypothesis of small fluctuations of the synaptic current. The resulting theory allows us to prove the formal equivalence between local and distributed inertia, holding for any synaptic time scale, integrate-and-fire neuron model, spike emission regimes and for different network states even when the neuron number is finite.

Tables

Table 1. Table 1: Parameters of the VIF neuron networks at drift-dominated (DD) and noise-dominated (ND) regimes.

Parameter	Value (DD)	Value (ND)
Number of excitatory neurons $N$	2000	2000
Mean excitatory synapses per neuron $C$	100	100
Excitatory PSP amplitude due to recurrent spikes $J$ [$v_{\mathrm{thr}}$]	0.0075	0.014
Minimum transmission delay ${\delta }_{min}$ [ms]	10	10
External spike rate $C_{ext}{\nu }_{ext}$ [Hz]	3000	3000
Excitatory PSP amplitude due to external spikes $J_{ext}$ [$v_{\mathrm{thr}}$]	0.0114	0.0726
Current bias $I_{DC}$ [$v_{\mathrm{thr}}/$ms]	$-0.032$	$-0.242$
Refractory period ${\tau }_{ref}$ [ms]	0	0
Firing threshold $v_{\mathrm{thr}}$	1	1
Resting potential $V_E$ [$v_{\mathrm{thr}}$]	0	0
Reset potential $v_{\mathrm{res}}$ [$v_{\mathrm{thr}}$]	0	0
Firing rate $\nu \ast$ [Hz]	10	10

Table 2. Table 2: Parameters of the simulated LIF neuron networks at drift-dominated (DD) and noise-dominated (ND) regimes.

Parameter	Value (DD)	Value (ND)
Number of excitatory neurons $N$	2000	2000
Mean excitatory synapses per neuron $C$	100	100
Excitatory PSP amplitude due to recurrent spikes $J$ [mV]	0.140	0.213
Minimum transmission delay ${\delta }_{min}$ [ms]	3	3
External spike rate $C_{ext}{\nu }_{ext}$ [Hz]	12000	12000
Excitatory PSP amplitude due to external spikes $J_{ext}$ [mV]	0.272	1.518
Current bias $I_{DC}$ [nA]	$-1.23$	$-9.16$
Membrane time constant ${\tau }_m$ [ms]	20	20
Refractory period ${\tau }_{ref}$ [ms]	0	0
Firing threshold $v_{\mathrm{thr}}$ [mV]	20	20
Membrane capacitance $C_m$ [pF]	500	500
Resting potential $V_E$ [mV]	0	0
Reset potential $v_{\mathrm{res}}$ [mV]	0	0
Firing rate $\nu \ast$ [Hz]	40	40

Table 3. Table 3: Parameters of the EIF neuron networks at drift-dominated (DD) and noise-dominated (ND) regimes.

Parameter	Value (DD)	Value (ND)
Number of excitatory neurons $N$	2000	2000
Mean excitatory synapses per neuron $C$	100	100
Minimum transmission delay ${\delta }_{min}$ [ms]	3	3
External spike rate $C_{ext}{\nu }_{ext}$ [Hz]	12000	12000
Excitatory PSP amplitude due to external spikes $J_{ext}$ [mV]	0.0804	1.76
Spike rate from external neurons $C_{ext}$ [Hz]	12000	12000
Current bias $I_{DC}$ [pA]	$-8.46$	$-223$
Membrane time constant ${\tau }_m$ [ms]	10	10
Refractory period ${\tau }_{ref}$ [ms]	1.7	1.7
Spike slope factor ${\Delta }_T$ [mV]	3.48	3.48
Firing threshold $v_{\mathrm{thr}}$ [mV]	$-59.9$	$-59.9$
Membrane capacitance $C_m$ [pF]	500	500
Resting potential $V_m$ [mV]	$-65$	$-65$
Leak potential $E_L$ [mV]	$-65$	$-65$
Reset potential $v_{\mathrm{res}}$ [mV]	$-68$	$-68$
Leak conductance $g_L=C_m/{\tau }_m$ [nS]	50	50
Excitatory PSP amplitude due to recurrent spikes $J$ [mV]	0.134	0.272
$V_{peak}$ [mV]	0	0
Firing rate $\nu \ast$ [Hz]	40	40

Table 4. Table 4: Parameters of the LIF neuron networks used for the Hopf bifurcation test (Fig.  5 ) when τ s = τ d = 4 subscript 𝜏 s subscript 𝜏 d 4 {\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}=4~{} ms (T4) and τ s = τ d = 16 subscript 𝜏 s subscript 𝜏 d 16 {\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}=16~{} ms (T16).

Parameter	Value (T4)	Value (T16)
Number of inhibitory neurons $N$	20000	20000
Mean inhibitory synapses per neuron $C$	2000	2000
Excitatory PSP amplitude due to external spikes $J_{ext}$ [mV]	0.1	0.1
Membrane time constant ${\tau }_m$ [ms]	20	20
Refractory period ${\tau }_{ref}$ [ms]	0	0
Firing threshold $v_{\mathrm{thr}}$ [mV]	20	20
Membrane capacitance $C_m$ [pF]	500	500
Resting potential $V_E$ [mV]	0	0
Reset potential $v_{\mathrm{res}}$ [mV]	10	10
Firing rate $\nu \ast$ [Hz]	100	40

Table 5. Table 5: Set of J 𝐽 J and I D ​ C subscript 𝐼 𝐷 𝐶 I_{DC} values used in the network of inhibitory LIF neurons described in Tab.  4 to produce the results plotted in Fig.  5 .

$J$ [mV]	$I_{DC}$ [pA]
$-0.0025$	$-2020$
$-0.0040$	$-1990$
$-0.0063$	$-1940$
$-0.0100$	$-1870$
$-0.0158$	$-1760$
$-0.0250$	$-1570$
$-0.0628$	$-838$
$-0.0995$	$-134$
$-0.158$	$970$
$-0.250$	$2700$

Equations

$$\dot{V}=-f(V)/{\tau }_m+I+I_{ext},{\tau }_s\dot{I}=-I+{\mu }_I+{\sigma }_I\xi ,$$

$$\begin{array}[]{r}\mu_{I}(t)=J\,C\,\nu(t)\\ \sigma_{I}^{2}(t)=J^{2}\,C\,\nu(t)\end{array}$$

$$\begin{array}[]{rcl}\partial_{t}P&=&-\partial_{x}\mathcal{F}_{x}(x,y)-\partial_{y}\mathcal{F}_{y}(y,\nu)\\ &\equiv&\frac{1}{\tau_{m}}\,\partial_{x}\left[\left(f(x)-y-\mu_{ext}\right)+\frac{1}{2}\sigma_{ext}^{2}\partial_{x}\right]P+\\ &&\frac{1}{\tau_{s}}\,\partial_{y}\left[(y-\mu_{I})+\frac{1}{2}\sigma_{I}^{2}\partial_{y}\right]P\\ &\equiv&\left(\mathcal{L}_{x}+\mathcal{L}_{y}\right)P\,,\end{array}$$

$$\nu (t)={\int }_{-\infty }^{\infty }{\mathcal{F}}_x(v_{\mathrm{thr}},y)dy\equiv {\int }_{-\infty }^{\infty }\nu (y,t)dy,$$

$$\begin{array}[]{rcl}\dot{a}_{n}=\lambda_{n}a_{n}&+&\displaystyle\sum_{m}\langle\psi_{n}|\mathcal{L}_{y}\,a_{m}\phi_{m}\rangle\\ &-&\displaystyle\frac{\sigma^{2}(y)}{2\tau_{m}}\sum_{m}a_{m}\langle\psi_{n}|\partial_{x}^{2}\phi_{m}\rangle\,.\end{array}$$

$${\dot{\rho }}_y={\mathcal{L}}_y{\rho }_y-{\rho }_y\frac{{\sigma }^2(y)}{2{\tau }_m}{\int }_{v_{\min }}^{v_{\mathrm{thr}}}{\partial }_x^2{\varphi }_0(x,y)dx.$$

$${\dot{\rho }}_y={\mathcal{L}}_y{\rho }_y={\partial }_y\frac{(y-{\mu }_I){\rho }_y}{{\tau }_s}+\frac{{\sigma }_I^2}{2{\tau }_s}{\partial }_y^2{\rho }_y.$$

$$\nu (y,t)={\rho }_y(y,t)\Phi (y,\nu )+\sum _{n\ne 0}{a}_n(y,t)f_n(y,\nu ),$$

$${\tau }_s{\dot{\mu }}_y=-{\mu }_y+{\mu }_I,{\tau }_s\dot{{\sigma }_y^2}=-2{\sigma }_y^2+{\sigma }_I^2,$$

$$\begin{array}[]{rcl}\nu&=&\displaystyle\int_{-\infty}^{\infty}\nu(y,t)\,dy\\ &=&\displaystyle\Phi(\mu_{y},\nu)\,\int_{-\infty}^{\infty}\rho_{y}(y,t)\,dy\,+\\ &&\displaystyle\sum_{n\neq 0}f_{n}(\mu_{y},\nu)\,\int_{-\infty}^{\infty}a_{n}(y,t)\,dy+\mathcal{O}(\sigma_{I})\\ &\equiv&\displaystyle\Phi_{0}(\nu)+\sum_{n\neq 0}a_{0n}\,f_{0n}(\nu)+\mathcal{O}(\sigma_{I})\,,\end{array}$$

$$\begin{array}[]{rcl}\dot{a}_{0n}&=&\displaystyle\int_{-\infty}^{\infty}\lambda_{n}(y)\,a_{n}(y,t)\,dy\,+\\ &&\displaystyle\sum_{m}\int_{-\infty}^{\infty}\langle\psi_{n}|\partial_{y}(\frac{y-\mu_{I}}{\tau_{s}}a_{m}\phi_{m})\rangle\,dy+\mathcal{O}(\sigma_{I}^{2})\\ &=&\lambda_{0n}\,a_{0n}-\displaystyle\frac{\mu_{y}-\mu_{I}}{\tau_{s}}\sum_{m}\langle\partial_{y}\psi_{n}|\phi_{m}\rangle|_{y=\mu_{y}}a_{0m}+\mathcal{O}(\sigma_{I})\,,\end{array}$$

$${\dot{a}}_{0n}={\lambda }_{0n}{a}_{0n}+{\dot{\mu }}_y\sum _m⟨{\partial }_y{\psi }_n|{\varphi }_m⟩{|}_{y={\mu }_y}{a}_{0m}+\mathcal{O}({\sigma }_I).$$

$$\{\begin{array}{l}{\dot{\vec{a}}}_0=\left({\mathbf{\Lambda }}_0+{\mathbf{W}}_0{\dot{\mu }}_y\right){\vec{a}}_0+{\vec{w}}_0{\dot{\mu }}_y\\ {\dot{\mu }}_y=\left({\mu }_I(\nu )-{\mu }_y\right)/{\tau }_s\\ \nu ={\Phi }_0+{\vec{f}}_0\cdot {\vec{a}}_0\end{array}.$$

$$\tilde{\nu }=\frac{{\mu }_y}{JC}$$

$$\{\begin{array}{l}\dot{\vec{a}}=\left({\mathbf{\Lambda }}_0+{\mathbf{W}}_{0}JC\dot{\tilde{\nu }}\right)\vec{a}+{\vec{w}}_{0}JC\dot{\tilde{\nu }}\\ \dot{\tilde{\nu }}=\left(\nu -\tilde{\nu }\right)/{\tau }_{\delta }\\ \nu ={\Phi }_0+{\vec{f}}_0\cdot \vec{a}\end{array},$$

$$\rho (\delta )=\frac{1}{{\tau }_{\mathrm{d}}}{e}^{-\delta /{\tau }_{\mathrm{d}}}\Theta (\delta )$$

$$P_{\nu }(\omega )=\frac{1+2\mathrm{Re}\left[{\vec{f}}_0\cdot (i\omega \mathbf{I}-{\mathbf{\Lambda }}_0{)}^{-1}{\vec{\psi }}_0\right]}{{\left|1-\left[{\Phi }_0^{\prime }+i\omega {\vec{f}}_0\cdot (i\omega \mathbf{I}-{\mathbf{\Lambda }}_0{)}^{-1}{\vec{w}}_0\right]\rho (i\omega )\right|}^2}\frac{{\nu }^{\ast }}{N}.$$

$$\{\begin{array}{l}{\dot{\vec{a}}}_0=\left({\mathbf{\Lambda }}_0+{\mathbf{W}}_0{\dot{\mu }}_y\right){\vec{a}}_0+{\vec{w}}_0{\dot{\mu }}_y+{\vec{\psi }}_0{\eta }_N\\ {\dot{\mu }}_y=\left({\mu }_I({\nu }_N)-{\mu }_y\right)/{\tau }_s\\ \nu ={\Phi }_0+{\vec{f}}_0\cdot {\vec{a}}_0\\ {\nu }_N=\nu +{\eta }_N\end{array},$$

$$P_{\nu }(\omega )=\frac{1+2\mathrm{Re}\left[{\vec{f}}_0\cdot (i\omega \mathbf{I}-{\mathbf{\Lambda }}_0{)}^{-1}{\vec{\psi }}_0\right]}{{\left|1-\left[{\Phi }_0^{\prime }+i\omega {\vec{f}}_0\cdot (i\omega \mathbf{I}-{\mathbf{\Lambda }}_0{)}^{-1}{\vec{w}}_0\right]\rho (i\omega )\right|}^2}\frac{{\nu }^{\ast }}{N},$$

$$\frac{\mathrm{d}}{\mathrm{d}\nu }\Phi (\mu ,{\sigma }^2)=CJ\left(\frac{\partial \Phi }{\partial \mu }+J\frac{\partial \Phi }{\partial {\sigma }^2}\right)$$

$$\{\begin{array}{l}\mu =JC{\nu }^{\ast }+J_{ext}{C}_{ext}{\nu }_{ext}+I_{DC}\\ {\sigma }^2=J^{2}C{\nu }^{\ast }+J_{ext}^2{C}_{ext}{\nu }_{ext}\end{array}.$$

Full text

I Introduction

Large distributed systems like brain neuronal networks often have to satisfy both timing and space constraints irrespective of their size. Indeed, they have to be compact in order to be portable, and at the same time their components (i.e., the neurons) have to communicate always to the same pace, which in turn is determined by the environment Buzsáki et al. (2013). The parallel non-instantaneous communication of messages like the spikes emitted by the neurons can be seen as an effective inertia and can be increased acting locally, by slowing down the dynamics of the single elements composing the networks, or globally, by distributing inertia across the system elements and their connections. A large local inertia implies a long memory (i.e., a long decay time bringing to forget initial conditions) and usually a big size, as in the case of a large capacitance in an electronic circuit or a large mass in a given physical system. This aspect usually obstacles the miniaturization of a many-body system made up of elements (with a given memory capacity) with large local inertia. Indeed, the bigger the overall memory capacity, the larger the required space and/or the consumed energy per element Mead (1989); Sarpeshkar (1998); Laughlin and Sejnowski (2003). Alternatively, a global inertia can be implemented by delaying in time the interaction between the network components, thus implementing a distributed memory of the collective activity spread across the whole system, possibly hindering fast synchronizations and hence making the system slower and more stable Atay (2003); Mattia and Del Giudice (2003). Also this strategy has drawbacks, such as the implementation of delay lines with their related costs. For instance, a delayed distribution of the past activity/state of the system requires a memory resource or an additional energy demand in order to keep alive the messages to be exchanged between the network components Laughlin and Sejnowski (2003); Merolla et al. (2014). Maybe it is not by chance that both these strategies coexist in a huge cellular network like the one composing our brain. Indeed, neuronal networks have both local low-pass filters at the level of synapses – which select preferentially slow-varying information and provide local inertia – and a global spread of their instantaneous activity through the axon and dendrites, which propagate emitted spikes with suited delays and provide global inertia. Then a natural question arises: are these two kinds of inertia affecting differently the network dynamics?

This is a challenging question. Indeed, synapses in neuron networks filter incoming spikes with a wide variety of time scales, affecting the stability of various collective dynamics Wang (1999), the selectivity in transmitting information Abbott and Regehr (2004) and the reactivity to suddenly appearing exogenous stimuli Fourcaud and Brunel (2002). Even for extremely simplified models of spiking neurons and network connectivities, theoretical approaches including non-instantaneous transmission rely on approximations valid for relatively small time scales Brunel and Sergi (1998); Brunel et al. (2001); Fourcaud and Brunel (2002); Schuecker et al. (2015), or for quasi-adiabatic dynamical regimes Moreno-Bote and Parga (2004, 2006), or for neurons working with low-noise supra-threshold inputs Lindner (2004); Schwalger et al. (2015). Despite the long history of attempts in statistical physics to work out an effective one-dimensional Markovian description of the dynamics of a non-Markovian system Risken (1984); Hänggi and Jung (1995) (i.e., a suitable representation of the membrane potential dynamics of the neurons when synaptic filtering is incorporated), a theoretical framework valid for any correlation time of the synaptic input and any neuronal activity regime is still lacking if the boundary conditions describing the spike emission process are taken into account.

Here, we address this issue by developing a theoretical framework in which from small to large synaptic time scales the same dynamical model holds for the instantaneous firing rate of spiking neuron networks. Starting from the population density dynamics of single-neuron state variables under mean-field approximation, we extend to the colored-noise case the spectral expansion of the associated Fokker-Planck (FP) equation previously derived under white-noise assumption Mattia and Del Giudice (2002). Resorting to a kind of central moment closure method Levermore (1996), a zero-th order approximation of the population density dynamics is worked out, recovering the same theoretical description as for the white-noise case. In the new model, the equation coefficients are modulated by the low-pass filtered mean synaptic current. This dimensional reduction provides results showing a remarkable agreement with microscopic simulations and finally allows us to provide an answer to the challenging question asked above. This is done by proving a formal equivalence between local and distributed inertia: a suited distribution of spike transmission delays between neurons (global inertia) allows to fully reproduce the dynamical features of the same network having instead non-instantaneous synaptic transmission (local inertia). This equivalence is independent of integrate-and-fire (IF) neuron models (VIF, LIF, EIF), spike emission regimes (sub- and supra-threshold) and dynamical state (stable asynchronous fixed point and limit cycles) chosen for the network, even when finite-size fluctuations are taken into account.

II Local inertia in neuronal networks

As paradigmatic example of local inertia in a many-body distributed system, we use a network composed of $N$ integrate-and-fire (IF) neurons, each receiving an input current $I(t)$ resulting from the low-pass filtered linear combination of the spikes emitted at time $t_{j}$ by a subset of presynaptic cells in the network Destexhe et al. (1998); Brunel et al. (2001): $\tau_{s}\,\dot{I}=-I+J\,\sum_{j}\delta(t-t_{j})$. For simplicity, here we assume a first-order dynamics for synaptic current $I$ with decay time $\tau_{s}$ and constant efficacy $J$ (Fig. 1a-b). The membrane potential $V(t)$ of an IF neuron evolves according to the general equation $\tau_{m}\,\dot{V}=-f(V)+R\,(I+I_{ext})$, with membrane decay constant $\tau_{m}$, neuron resistance $R$ (henceforth set to 1 for simplicity) and a voltage-dependent leakage drift $f(V)$ which depends on the model neuron type. For instance, $f(V)$ is a constant drift for the perfect IF (PIF) neuron and $f(V)=V$ for the leaky IF (LIF) Burkitt (2006). Neurons may receive an additional current $I_{ext}$ modeling external sources, like incoming synaptic input from other networks. Once $V$ crosses a threshold value ${v_{\mathrm{thr}}}$, a spike is emitted and potential $V$ is reset to ${v_{\mathrm{res}}}<{v_{\mathrm{thr}}}$. Due to these boundary conditions, even under stationary regimes single-neuron dynamics is not analytically tractable. This is particularly apparent looking at the probability density $P(V,I)$ of states $(V,I)$ from a long time series worked out numerically by modeling presynaptic activity as a Poissonian spike train (Fig. 1c). The absorbing barrier in ${v_{\mathrm{thr}}}$ and the reentering flux of realizations in ${v_{\mathrm{res}}}$ modeling the emission of spikes, make $P(V,I)$ asymmetric Brunel and Sergi (1998); Fourcaud and Brunel (2002); Apfaltrer et al. (2006) and correlation between $V$ and $I$ non-monotonic, a feature even more apparent for relatively slow synaptic filtering (Fig. 1d). Under diffusion approximation, holding for large rate of incoming spikes each only mildly affecting $V$ Tuckwell (1988); Burkitt (2006), membrane potential dynamics is described by the following system of Langevin equations Risken (1984)

[TABLE]

where for simplicity currents are scaled by a factor $\tau_{m}$ (i.e., measured as voltage per unit time) and $I_{ext}=\mu_{ext}+\sigma_{ext}\,\xi_{ext}$ is a Gaussian white noise homogeneous in time.

When the mean driving force alone (i.e., the deterministic part of the input $I+I_{ext}$) is not enough to make the membrane potential $V$ cross the threshold $v_{thr}$, the neurons are evolving in a noise-dominated (or subthreshold) regime and irregular firing occurs, due to the stochastic part of the input. In the opposite case, the emission of an action potential can occur also in the absence of noisy afferent currents, and the neurons are in a drift-dominated (or suprathreshold) regime of activity, characterized by the regular emission of spike trains Fusi and Mattia (1999); Brunel (2000); Lansky and Ditlevsen (2008).

Equation (1) describes the time course of a generic realization $(V,I)$, which under mean-field approximation represents the dynamics of the statistically identical neurons of the network. According to this theoretical ansatz, infinitesimal moments

[TABLE]

characterize the white noise $\xi$ driving the activity-dependent synaptic current Amit and Brunel (1997) as a function of the network emission rate $\nu(t)$ (i.e., the average number of spikes each neuron emits per time unit) and the average number of presynaptic contacts $C$ ($\leq N$) per neuron. The probability density $P(x,y,t)$ to find neurons at time $t$ with membrane potential $V(t)=x$ and synaptic current $I(t)=y$ follows the FP equation Risken (1984)

[TABLE]

a continuity equation in which the density changes are determined by the divergence of the probability current $(\mathcal{F}_{x},\mathcal{F}_{y})$ (i.e., the flux of realizations) together with few specific boundary conditions (see Fig. 1): i) spike emission as an absorbing barrier at $x={v_{\mathrm{thr}}}$ (only for $y>0$ if $\sigma_{ext}=0$), ii) reentering flux of absorbed realizations at the reset potential ${v_{\mathrm{res}}}$ and iii) lower bound for the membrane potential implemented as a reflecting barrier at ${v_{\mathrm{min}}}$ Knight et al. (1996); Brunel and Sergi (1998); Brunel and Hakim (1999); Mattia and Del Giudice (2002).

In general, Eq. (3) is hard to solve and only approximated solutions (mainly under stationary conditions) have been worked out Risken (1984); Doering et al. (1987); Hänggi and Jung (1995); Fourcaud and Brunel (2002). This is both due to the peculiar boundary conditions and to its intrinsic nonlinearity, as $\nu(t)$ depends on $P(x,y,t)$, being the flux of realizations crossing ${v_{\mathrm{thr}}}$:

[TABLE]

where $\nu(y,t)=\mathcal{F}_{x}({v_{\mathrm{thr}}},y)$ is the ‘partial’ rate of spikes emitted by subpopulation $y$. Here we derive a solution for this two-dimensional problem without making any assumption like stationarity in time and restricted ranges for the synaptic filtering timescale. We extend the spectral expansion approach exploited to study the one-dimensional case Knight et al. (1996); Mattia and Del Giudice (2002), by describing $P(x,y,t)$ as the superposition of an infinite set of one-dimensional interacting sub-populations, labeled by $y$. Sub-population dynamics are studied, as in Mattia and Del Giudice (2002), by projecting $P$ at fixed $y$ on the eigenfunctions $\{|\phi_{m}\rangle\}$ of a suited one-dimensional FP operator, which here we set to be $\mathcal{L}_{xy}=\frac{1}{\tau_{m}}\partial_{x}[f(x)-y-\mu_{ext}]+\frac{1}{2\tau_{m}}\left[\sigma^{2}_{ext}+\sigma^{2}(y)\right]\partial_{x}^{2}$, with an additional diffusion term $\sigma^{2}(y)$ with respect to $\mathcal{L}_{x}$ of Eq. (3). This allows to rewrite this equation in terms of the expansion coefficients $a_{n}(y,t)=\langle\psi_{n}|P\rangle=\int_{{v_{\mathrm{min}}}}^{v_{\mathrm{thr}}}\psi_{n}(x,y)\,P(x,y,t)\,dx$, which now explicitly depend on the auxiliary variable $y$ and where $\langle\psi_{n}|$ is the $n$-th eigenfunction of the adjoint operator of $\mathcal{L}_{xy}$. The stationary mode with eigenvalue $\lambda_{0}=0$ has $\langle\psi_{0}|=1$, thus $a_{0}(y,t)$ is the marginal probability density $\rho_{y}(y,t)=\int_{{v_{\mathrm{min}}}}^{v_{\mathrm{thr}}}P(x,y,t)\,dx$.

The dynamics of $a_{n}$ can be worked out by deriving $a_{n}$ with respect to time Knight et al. (1996); Mattia and Del Giudice (2002), and rewriting $\mathcal{L}_{x}=\mathcal{L}_{xy}-\frac{1}{2\tau_{m}}\sigma^{2}(y)\partial_{x}^{2}$:

[TABLE]

For $n=0$, the dynamics of $\rho_{y}(y,t)=a_{0}(y,t)$ can be worked out, as both $\mathcal{L}_{y}$ and $a_{n}$ do not depend on $x$:

[TABLE]

The integral on the r.h.s. can be solved by parts and turns out to be 0 thanks to both the conservation of the flux of realizations exiting from ${v_{\mathrm{thr}}}$ and re-entering in ${v_{\mathrm{res}}}$ ($\partial_{x}\phi_{0}|_{x={v_{\mathrm{thr}}}}=\partial_{x}\phi_{0}|_{x={v_{\mathrm{res}}}^{+}}-\partial_{x}\phi_{0}|_{x={v_{\mathrm{res}}}^{-}}$) and the reflecting barrier condition in ${v_{\mathrm{min}}}$, provided that the reasonable assumption $\phi_{0}({v_{\mathrm{min}}},y)=0$ holds (this is surely true for ${v_{\mathrm{min}}}\to-\infty$, as $\phi_{0}$ is a probability density). Thus, the $\rho_{y}$ dynamics reduces to

[TABLE]

Relying on the same spectral expansion as in Mattia and Del Giudice (2002), each ‘partial’ rate $\nu(y,t)$ in Eq. (4) can be decomposed as follows:

[TABLE]

where the stationary mode contribution is separated from the others. Here, $\Phi(y,\nu)\equiv\mathcal{F}_{x}({v_{\mathrm{thr}}},y)|_{P=\phi_{0}}$ is the emission rate of the stationary one-dimensional density $\phi_{0}$ at fixed $y$, while the other nonstationary modes contribute with the fluxes $f_{n}(y,\nu)\equiv\mathcal{F}_{x}({v_{\mathrm{thr}}},y)|_{P=\phi_{n}}$. Both $\Phi$ and $f_{n}$ depend on the total emission rate $\nu(t)$ through the input current moments $\mu_{I}$ and $\sigma_{I}$.

III Dimensional reduction of the emission rate dynamics

Equations (4-7) together are a mere reformulation of the original problem (3), not a solution. However, in this framework a perturbative approach can be envisaged. From Eq. (6), synaptic current $y$ results to have time-dependent mean $\mu_{y}(t)=\langle y\rangle$ and variance $\sigma_{y}^{2}(t)=\langle\left[y-\mu_{y}(t)\right]^{2}\rangle$

[TABLE]

like an Ornstein-Uhlenbeck process with nonstationary input moments $\mu_{I}(t)$ and $\sigma_{I}(t)$ Risken (1984). Hence, synaptic current displacement in time is $|y-\mu_{y}|=\mathcal{O}(\sigma_{I})$ and, for small $\sigma_{I}$, we can assume a negligible role for a subpopulation labeled with $y$ distant enough from $\mu_{y}$. This assumption of a narrow $\rho_{y}$ centered around $y=\mu_{y}$ allows to simplify the integrals across the $y$ domain weighted by the expansion coefficients $a_{n}(y,t)$. Indeed, approximations of these integrals can be worked out by expanding in Taylor’s series the $y$-dependent functions in the integrands. By combining Eqs. (4) and (7), the total emission rate can be then rewritten as

[TABLE]

where the [math]-th order terms of the Taylor expansions are defined as $\Phi_{0}(\nu)\equiv\Phi(\mu_{y},\nu)$ and $f_{0n}(\nu)\equiv f_{n}(\mu_{y},\nu)$, together with the introduction of the new integral expansion coefficients $a_{0n}(t)\equiv\int_{-\infty}^{\infty}a_{n}(y,t)\,dy$. The next step to find a self-consistent approximated spectral expansion of Eq. (3) is to work out the dynamics of these $a_{0n}(t)$. This can be done by integrating both sides of Eq. (5) and constraining the additional diffusion term in $\mathcal{L}_{xy}$ to be as small as the fluctuation size of $y$: $\sigma^{2}(y)=\mathcal{O}(\sigma_{I}^{2})$. The dynamics of the integral expansion coefficients reduces to

[TABLE]

where we integrated by parts the second integral and set $\lambda_{0n}\equiv\lambda_{n}(\mu_{y})$. This equation for $a_{0n}$ can be further recast as follows, by taking into account Eq. (8):

[TABLE]

This is a dimensional reduction of Eq. (5) as the coefficients $a_{0n}$ do not depend on $y$, and it holds provided that current $I$ has a narrow distribution across neurons at each time $t$ (i.e., small $\sigma_{I}(t)$), leaving $\tau_{s}$ unconstrained.

To work out a self-consistent equation for emission rate $\nu(t)$, we still need to choose a suited $\sigma(y)$ in $\mathcal{L}_{xy}$. To this purpose, the case of instantaneous synaptic transmission $\tau_{s}=0$ can be used as reference, as Eq. (3) reduces to a one-dimensional FP equation with the only operator $\mathcal{L}_{x0}=-1/\tau_{m}\,\partial_{x}\left[\mu_{I}+\mu_{ext}-f(x)\right]+\left(\sigma_{ext}^{2}+\sigma_{I}^{2}\right)/(2\tau_{m})\partial_{x}^{2}$ Brunel and Hakim (1999); Fusi and Mattia (1999); Mattia and Del Giudice (2002). Hence, $\mathcal{L}_{xy}|_{y=\mu_{y}}$ must tend to $\mathcal{L}_{x0}$ for vanishing $\tau_{s}$, which is the case if $\sigma^{2}(y)=J\,y$. Indeed, in this limit $\sigma^{2}(\mu_{y})=J^{2}\,C\,\nu=\sigma_{I}^{2}$, being $\mu_{y}=\mu_{I}=J\,C\,\nu$ from Eq. (8). This implies that eigenvalues $\lambda_{0n}$ and eigenfunctions in Eq. (10) are the same as in the one-dimensional case with $\delta$-correlated synaptic input, but with collective emission rate $\nu(t)=\mu_{y}(t)/(J\,C)$.

In conclusion, this dimensional reduction extends the emission rate equation previously worked out for $\tau_{s}=0$ Mattia and Del Giudice (2002), and, by combining Eqs. (8–10) and neglecting $\mathcal{O}(\sigma_{I})$ terms, it results to be

[TABLE]

Here, the infinite vectors $\vec{a}_{0}=\{a_{0n}\}$ and $\vec{f}_{0}(\nu)=\{f_{0n}(\nu)\}$ are introduced together with the eigenvalue matrix $\mathbf{\Lambda}_{0}=\operatorname{diag}(\lambda_{0n})$. Separating stationary from non-stationary modes, the synaptic coupling vector $\vec{w}_{0}=\{\langle\partial_{y}\psi_{n}|\phi_{0}\rangle|_{y=\mu_{y}}\}$ and matrix $\mathbf{W}_{0}=\{\langle\partial_{y}\psi_{n}|\phi_{m}\rangle|_{y=\mu_{y}}\}_{m\neq 0}$ are also included, respectively. For the above chosen $\sigma(y)$, we remark that all these terms depend on the total emission rate $\nu(t)$ only through $\mu_{y}(t)$.

IV Local and distributed inertia equivalence

So far, we considered the filtering activity operated by local synapses transmitting incoming spikes as a post-synaptic potential non-instantaneous in time. Here, we recall the known dynamics of a network of spiking neurons where synaptic transmission is instantaneous but a distribution of axonal delays is taken into account Brunel and Hakim (1999); Mattia and Del Giudice (2002, 2003). In this one-dimensional case ($\tau_{s}=0$), the network dynamics can be simply worked out by replacing $\nu$ in $\mu_{I}$ and $\sigma_{I}^{2}$ (see Eqs. (2)) with the instantaneous rate $\tilde{\nu}(t)=\int_{0}^{\infty}\nu(t-\delta)\rho(\delta)\mathrm{d}\delta$ of spikes received by neurons when a distribution $\rho(\delta)$ of axonal transmission delays $\delta$ is taken into account. Interestingly, the coefficients of the related emission rate equation result to have a straightforward relationship with those in Eq. (11). Indeed, the elements of the synaptic coupling vector and matrix now are $\langle\partial_{\tilde{\nu}}\psi_{n}|\phi_{m}\rangle=J\,C\langle\partial_{y}\psi_{n}|\phi_{m}\rangle|_{y=\mu_{y}}$ for any $n$ and $m$, since

[TABLE]

for what shown above. Due to this, the searched emission rate equation reduces to

[TABLE]

where the specific delay distribution

[TABLE]

has been taken into account. This leads $\tilde{\nu}$ to be a version of the collective emission rate $\nu$ smoothed in time by a first-order low-pass filter with decay time $\tau_{\delta}$. $\Theta(\delta)$ is the Heaviside function as only positive transmission delays are admitted.

A comparison between Eqs. (11) and (12) remarkably points out the equivalence of having, in a spiking neuron network, a non-instantaneous synaptic transmission or a suited distribution of axonal delays. In Eq. (12) the expansion coefficients in $\vec{a}$ refer to the one-dimensional density $P(x,t)$, while in Eq. (11) $\vec{a}_{0}$ is only an effective representation of the two-dimensional density $P(x,y,t)$. Nevertheless, provided that ${\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}$, both the dynamics seen from the perspective of $\nu(t)$ are the same, as the functions $\Phi_{0}$, $\vec{f}_{0}$, $\vec{w}_{0}$, $\mathbf{W}_{0}$ and $\mathbf{\Lambda}_{0}$ are the same. In other words, under mean-field approximation and for not too large synaptic current fluctuations $\sigma_{I}$, having a local synaptic filtering with cut-off frequency $1/\tau_{s}$ is equivalent to have random axonal delays with exponential distribution (13) and decay constant ${\tau_{\mathrm{d}}}={\tau_{\mathrm{s}}}$.

V Testing inertia equivalence in simulation

The theory developed in the previous sections has been tested through extensive numerical simulations by using NEST Gewaltig and Diesmann (2007) and the high-performance custom simulator implementing the event-based approach described in Mattia and Del Giudice (2000). The parameters for the numerical simulations have been identified following the procedure described in Appendix VII.2. The network parameters corresponding to the results proposed in this section are summarized in Appendix VII.3.

V.1 Equivalence in the asynchronous state of finite-size networks

A first evaluation of the effectiveness of the low-dimensional description derived above is possible by inspecting the activity in the presence of an endogenous noise in a network of a finite number $N$ of neurons trapped into a stable asynchronous state. As previously done in Mattia and Del Giudice (2002), finite-size fluctuations can be incorporated into Eqs. (11) and (12) as an additive forcing term to the expansion coefficient dynamics, giving rise to an equation for the emission rate $\nu_{N}(t)$ of a finite pool of neurons (see Appendix VII.1 for details). The resulting dynamics can be linearized around the fixed-point $\nu^{*}$, allowing to compute the Fourier transform $\nu_{N}(\omega)$, and from it the power spectral density $P_{\nu}(\omega)=|\nu_{N}(\omega)|^{2}$, which turns out to be

[TABLE]

Here, $\mathbf{I}$ is the identity matrix, all the coefficients depending on $\nu(t)$ are now constants computed at $\nu(t)=\nu^{*}$, $\Phi^{\prime}_{0}=\partial_{\nu}\Phi_{0}$ and $\rho(i\omega)$ is the ratio between the Fourier transforms of $\mu_{y}$ and $\mu_{I}$, which results to be $\rho(i\omega)=1/(1+i\omega{\tau_{\mathrm{s}}})$. Note that $\rho(i\omega)$ is the same as for the Fourier transform of the delay distribution (12) with ${\tau_{\mathrm{d}}}={\tau_{\mathrm{s}}}$. If we consider an additional fixed transmission delay $\delta$, this transform generalizes to $\rho(i\omega)=\exp(-i\omega\delta)/(1+i\omega{\tau_{\mathrm{s}}})$ Mattia and Del Giudice (2002, 2003).

We compared this theoretical result with the $P_{\nu}(\omega)$ estimated from simulations of networks composed of simple IF neurons with synaptic filtering (see Appendix VII.4 for details), finding a remarkable agreement (Fig. 2). Here we resorted to the VIF neuron model Fusi and Mattia (1999), an extended version of the widely used PIF neuron Gerstein and Mandelbrot (1964), for its amenability to analytical treatment Mattia and Del Giudice (2002) (see Appendix VII.3 for details). Besides the good matching between theory and simulations, it is interesting to note how resonant peaks are differently affected by an increasing ${\tau_{\mathrm{s}}}$. Not surprisingly, high-$\omega$ peaks are dampened due to the low-pass filtering of synaptic transmission. These resonances are related to the synaptic reverberation of spiking activity, which gives rise to the so-called transmission poles in the linearized dynamics Treves (1993); Mattia and Del Giudice (2002). In an excitatory network as in Fig. 2, these are expected to be found at frequencies multiples of the inverse of the average time needed by a presynaptic spike to affect postsynaptic potential, here roughly the same as $1/(\delta+{\tau_{\mathrm{s}}})$. On the other hand, the low-$\omega$ peaks do not display any shift in frequency, although their power is dampened as well. This is because such peaks are due to the so-called diffusion poles in the linearized dynamics Spiridon and Gerstner (1999); Mattia and Del Giudice (2002). They occur when neurons emit spikes at drift-dominated (suprathreshold) regime. In this regime, the distribution of the inter-spike intervals is narrow and resonant peaks emerge at $\omega/2\pi$ multiples of $\nu^{*}$ (equal to 10 Hz in Fig. 2). We remark that, due to the shifting at low-$\omega$ of the transmission peaks, a constructive interference with diffusion peaks may occur. This explains the non-monotonic change of power of the second diffusion peak at 20 Hz in Fig. 2. Indeed, its power decreases when ${\tau_{\mathrm{s}}}$ is increased from 0 ms to 4 ms, as expected, whereas due to the mentioned interference it is heightened for ${\tau_{\mathrm{s}}}=32$ ms. This is because the lowest transmission peak in this case is expected to be found at $\omega/2\pi=23.8$ Hz, not too far from $2\,\nu^{*}=20$ Hz.

V.2 Equivalence under noise- and drift-dominated regime

To further test the equivalence between local and distributed inertia, we investigated whether this match worked well not only for networks of neurons under drift-dominated spiking regime, as shown in Fig. 2, but also under fluctuation-driven regimes. Indeed, synaptic filtering may have significantly different effects in these two regimes Brunel et al. (2001); Moreno-Bote and Parga (2004).

Figure 3 shows for comparison the power spectral densities $P_{\nu}(\omega)$ obtained by varying ${\tau_{\mathrm{s}}}$ and ${\tau_{\mathrm{d}}}$ in the same excitatory VIF neuron networks in which either non-instantaneous synaptic transmission or a distribution of transmission delays was incorporated (Fig. 3a and b, respectively). Networks with neurons spiking at drift-dominated regimes (top panels) were the same as those shown in Fig. 2. A remarkable agreement between simulations with local and distributed inertia was confirmed as iso-power curves tightly overlapped (see Fig. 3c) for a wide range of filter timescales. The same good match between these two network types was apparent also under noise-dominated regime (bottom panels).

As expected in this regime, the diffusion poles become real numbers Mattia and Del Giudice (2002) and resonant peaks at frequency $\omega/2\pi$ multiples of firing rate $\nu^{*}$ disappear, although the high-$\omega$ resonances corresponding to transmission poles remain almost unaffected by this change of regime. It interesting to point out here that widening the filtering time window (i.e., increasing ${\tau_{\mathrm{s}}}$ and ${\tau_{\mathrm{d}}}$) leads to almost flat power spectra, compatibly with the flattening of the response amplitude of isolated neurons receiving colored noise as input currents Brunel et al. (2001); Schuecker et al. (2015).

We remark that in Fig. 3a-c comparisons rely on normalized spectra, i.e. $P_{\nu}(\omega)\,N/\nu_{0}$. On the one hand, this means that power spectra shapes does not depend on the particular inertia implemented in the network. On the other hand, this does not guarantee that average firing rates $\nu_{0}$ are the same and do not change with the inertia-related timescales. To address this issue, in Fig. 3d the measured $\nu_{0}$ are compared. As expected, the average firing rate when distribution of delay are incorporated does not change by varying ${\tau_{\mathrm{d}}}$. Mean-field fixed-point $\nu^{*}$ mildly overestimates $\nu_{0}$ and this is due to the assumptions underlying the diffusion approximation. Indeed, the difference $|\nu_{0}-\nu^{*}|$ shrinks as the number of synaptic contacts per neuron and the network size $N$ increase by keeping unchanged drift and diffusion coefficients of the related Langevin equation (this result is not shown in the figures). But what is important to note here is rather that the average firing rate $\nu_{0}$ measured for local and distributed inertia are different, and the difference increases with ${\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}$. This trend is particularly apparent under noise-dominated regime, and it is not completely unexpected. Indeed, in the small-${\tau_{\mathrm{s}}}$ limit Brunel and Sergi (1998); Moreno-Bote and Parga (2004); Schuecker et al. (2015) and in the long synaptic timescale condition Moreno-Bote and Parga (2004) macroscopic differences are known to exist. This phenomenon is not captured by our 0-th order perturbative approach, since the total flux of realizations crossing the emission threshold ${v_{\mathrm{thr}}}$, which contributes to the total emission rate, depends on the shape of the probability current distribution $\mathcal{F}_{x}({v_{\mathrm{thr}}},y)=\partial_{x}P|_{x={v_{\mathrm{thr}}}}=\nu(y,t)$ Brunel and Sergi (1998); Schuecker et al. (2015). Indeed, as shown in Figure 1c-d, the section of $P(x,y,t)$ close to emission threshold ${v_{\mathrm{thr}}}$ shrinks when ${\tau_{\mathrm{s}}}$ increases from 4 ms to 64 ms.In our perturbative approach this dependence is completely neglected, as the distribution is assumed to be a Dirac’s $\delta$.

V.3 Independence from spiking neuron models

The developed theory is of general applicability, i.e. it applies to a wide range of spiking neuron models. As a result, the equivalence between local and distributed inertia proved for a network of simplified VIF neurons together with the theoretical expression for the power spectra of $\nu(t)$ are both expected to hold also for networks of more realistic single-cell models like the leaky IF (LIF, Tuckwell (1988)) and the exponential IF (EIF, Fourcaud-Trocmé et al. (2003)) neurons.

To verify this expectation, we directly simulated networks of LIF and EIF neurons with a mean-field stable equilibrium point at $\nu^{*}=40\,\mathrm{Hz}$. In Fig. 4 the mean $P_{\nu}(\omega)\,N/\nu_{0}$ for these kinds of networks is compared both under drift-dominated (top panels) and noise-dominated (bottom) regimes. A remarkable agreement is apparent for both LIF and EIF neuron networks (Fig. 4a and b, respectively), confirming the generality of the developed approach. To further test the absence of any bias due to the specific neuron model chosen, we computed the relative differences between the power spectral densities when non-instantaneous synaptic filtering ($P_{\mathrm{syn}}$) and transmission delay distribution ($P_{\mathrm{del}}$) are taken into account. The cumulative distributions of these discrepancies across all tested inertia time scales (${\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}$) and Fourier frequencies ($\omega/2\pi$) are shown in Fig. 4c. Interestingly, no significant differences are visible for the three chosen models (VIF, LIF, EIF), further proving that in all regimes (drift- and noise-dominated) our dimensional reduction does not depend on the specific single neuron dynamics.

Starting from this, it is not surprising to note here that spectral features similar to those highlighted in Fig. 3 for VIF neuron networks are also displayed by LIF and EIF neuron networks. More specifically, resonant peaks at multiple frequencies of $\nu_{0}$ under drift-dominated regimes and resonances at higher-$\omega$ due to the transmission poles – i.e., those related to the average spike transmission delay – are also visible in both Fig. 4a and b in LIF and EIF neuron networks, respectively.

V.4 Equivalence beyond the asynchronous state

So far, we tested the validity of the equivalence between non-instantaneous synaptic transmission and delay inhomogeneity in linearizable dynamical regimes like the asynchronous state. To further assess the generality of our theoretical results, we simulated an inhibitory network of LIF neurons (see Appendix VII.3) undergoing a transition from a point-attractor dynamics to a regime of collective oscillations through a supercritical Hopf bifurcation as in Brunel and Hakim (1999).

Figure 5a shows three examples of population activity $\nu(t)$ for increasing magnitude $|J|$ of the inhibitory coupling when a non-instantaneous synaptic transmission with decay time ${\tau_{\mathrm{s}}}=4\,\mathrm{ms}$ is incorporated. For relatively weak couplings (lowest $|J|$), the network is in an asynchronous state and its activity fluctuates around the equilibrium point $\nu^{*}=10$ Hz (black line). By increasing the synaptic inhibition above a threshold value ($|J|>0.04$ mV), the network starts to display coherent global oscillations of the emission rate with increasing amplitude (light and dark gray lines).

We spanned the same range of synaptic strength $J$ making the Hopf bifurcation apparent in networks with both synaptic filtering and delay distribution, having care to use the same mean-field parameters (Fig. 5b). For each $J$ we carried out the related external excitatory current $I_{ext}$ needed to keep unchanged the fixed-point at $\nu^{*}=10$ Hz regardless of its stability (see details in Appendix VII.3). We compared the dynamics of the two network types by measuring the standard deviation $\sigma_{\nu}$ across time of the emission rate $\nu(t)$ under stationary conditions (i.e., discarding the first second of the time series). As a remarkable result, we found an almost perfect match between $\sigma_{\nu}$ measured in networks with both local and distributed inertia. Indeed, the relative mismatch between $\sigma_{\nu}$ is always lower than $1\%$ (compare solid and dashed lines). This also implies that regardless of the inertia type, the transitions across regimes roughly occurs at the same critical value ($J\simeq-0.04$ mV) with a discontinuity of the derivative $\partial\sigma_{\nu}/\partial J$, highlighting in both cases the same supercritical Hopf bifurcation. Similar results were obtained for longer filtering timescales (${\tau_{\mathrm{s}}}={\tau_{\mathrm{d}}}=16$ ms, Fig. 5c), although in this case significant $\sigma_{\nu}$ differences can be measured for large enough $J$.

Altogether, these results further prove the generality of the theoretical equivalence derived in Sec. IV, clearly pointing out that local and distributed inertia give rise to the same collective dynamics not only in a tameable condition like the linearizable asynchronous state, but also for nonstationary network states like the global oscillatory regime.

VI Discussion

Local and distributed inertia in networks of spiking neurons are equivalent, provided that non-instantaneous post-synaptic currents from incoming spikes (local inertia) have the same time course as the distribution of spike transmission delays across axons and dendrites (distributed inertia). This is one of the main findings of this work, and we proved such equivalence by developing a novel population density approach to the network dynamics of neurons with non-Markovian membrane potentials. This general method consists of a dimensional reduction of the two-dimensional population density dynamics arising from the Markovian embedding Risken (1984); Hänggi et al. (1985) of the membrane potential evolution. The effective one-dimensional description we obtained relies on both an extended mean-field approximation Amit and Brunel (1997) and the assumption of a relatively low synaptic noise, without requiring any additional constraint on the correlation time due to the incorporated non-instantaneous synaptic transmission. This is our second main result, from which we obtained the same emission rate equation (ERE) as the one arising from the spectral expansion of the one-dimensional FP equation for the population density in the absence of synaptic filtering Mattia and Del Giudice (2002), but with an important difference; drift and diffusion coefficients, related to mean and variance of the input currents, respectively, now depend on a low-pass filtered version of the instantaneous firing rate of the network. We tested both the generality of the developed theory and the equivalence between local and distributed inertia through extensive numerical simulations, finding a remarkable agreement with theoretical expectations in a variety of dynamical scenarios, including noise- and drift-dominated regimes, equilibrium states and collective oscillations, excitatory and inhibitory synaptic connections, and for a wide set of spiking neuron models.

VI.1 Comparison with previous studies

Reducing into an effective one-dimensional Markovian problem the dynamics of a non-Markovian system – which in turn can be embedded into a larger set of coupled Markovian processes – has a long history in statistical physics Risken (1984); Hänggi and Jung (1995). Of course, the resulting approximation schemes have unavoidable limitations in representing the true dynamics of the original non-Markovian problem. Nevertheless, many successes were accumulated in the late 80s of the past century Sancho et al. (1982); Hänggi et al. (1985); Fox (1986); Grigolini (1986); Jung and Hänggi (1987); Masoliver et al. (1987). More specifically, many of them focused on the nonstationary dynamics underlying the specific problem of the first-passage time (FPT), thus taking into account an absorbing barrier, similarly to the spike emission threshold in neuronal modeling. Remarkably, also in this FPT framework stochastic processes driven by ‘colored noise’ can be effectively approximated by one-dimensional Markovian processes, provided that the related FP equations incorporate a state-dependent nonlinear diffusion coefficient Hänggi et al. (1985); Masoliver et al. (1987). This coefficient is not time-varying only in the limit of small noise and/or small correlation times ${\tau_{\mathrm{s}}}$ Fox (1986); Grigolini (1986). Other perturbative approaches have been worked out to cover the long-${\tau_{\mathrm{s}}}$ range Jung and Hänggi (1987); H’walisz et al. (1989), but the effectiveness of such Markovian approximations is mainly limited by the lack of a proper management of the boundary conditions Doering et al. (1987); Hänggi and Jung (1995); Kłosek and Hagan (1998). Intriguingly, here by including the proper absorbing barrier prescription due to the presence of a forcing white noise with arbitrarily small size $\sigma_{ext}$ (similarly to Moreno et al. (2002)), the theoretical derivation of the approximated ERE (11) eventually recovers a representation with time-varying state-dependent diffusion coefficients (see for instance Section V.B of Hänggi and Jung (1995)). Indeed, all the coefficients of this ERE depend on a filtered version $\mu_{y}(t)$ of the total emission rate $\nu(t)$, which from Eq. (4) is in turn a nonlinear functional of the probability density $P(x,y,t)$.

The approximated population density approach here developed brought also another interesting result. Typically, power spectra and Fourier transfer functions of the population activity in network of spiking neurons and synaptic filtering are characterized outside biologically relevant regimes. This is because theoretical approaches are perturbative and target extreme conditions, in which either fast Brunel et al. (2001); Fourcaud and Brunel (2002); Schuecker et al. (2015) or slow Moreno-Bote and Parga (2004, 2006) synapses are considered, or in which neurons work under strong drift-dominated (low-noise) regimes Lindner (2004); Schwalger et al. (2015). Here, we bridge this gap by presenting a theoretical description valid for both noise- and drift-dominated regimes, and for the whole spectrum of synaptic time scales and Fourier frequencies (Figs. 2–4).

It is also interesting to remark that the proved equivalence between local and distributed inertia allows to map a many-body system with non-Markovian coupled elements onto another many-body system where many Markovian units have delayed interactions. In the latter framework, different types of collective dynamics and, in particular, synchronization have been widely studied Boccaletti et al. (2006). The related theoretical results are then expected to be exportable to the case in which an equivalent synaptic filtering is taken into account. An intriguing example is the use of effective time-delayed feedback loops to control neuronal pathological states like the synchronous firing observed in Parkinson’s disease or epilepsy Schiff et al. (1994); Rosenblum and Pikovsky (2004). The same strategy could be implemented by relying on a suitably slow synaptic feedback, by taking advantage of the natural plasticity of the inter-neuronal couplings. Due to this neurophysiological adaptability, the optimal delayed feedback needed to control pathological dynamics could be in principle ‘learned’, similarly to what occurs in ‘echo-state networks’ and ‘liquid-state machines’ Jaeger and Haas (2004); Maass et al. (2007); Sussillo and Abbott (2009).

VI.2 Dimensional reduction

Other dimensional reductions of the Fokker-Planck dynamics arising when synaptic filtering is incorporated in network of spiking neurons have been proposed in the past to obtain effective kinetic representations amenable to numerical integration Cai et al. (2004); Rangan and Cai (2006); Ly and Tranchina (2007). They mainly rely on the centered moment closure method usually applied to kinetic problems in statistical physics Levermore (1996). Drawbacks of this approach are related to their limited applicability to noise-dominated regimes Ly and Tranchina (2007); Ly (2013), which are of particular relevance in neuroscience, being associated to activity states with balanced excitation and inhibition van Vreeswijk and Sompolinsky (1996); Shu et al. (2003); Yizhar et al. (2011). Even when some of these limitations are removed by introducing modifications to the standard mean-field approximation Ly (2013), this kind of dimensional reduction mainly remains a computational method in which, in addition to the dynamics of the centered moments, a one-dimensional partial differential equation has to be numerically integrated. On the contrary, our approach, in which a spectral expansion of the two-dimensional FP operator is operated together with a suited centered moment closure, is not only amenable to numerical integration. For instance, theoretical insights on the dynamics of the investigated system can be effectively carried out Mattia and Del Giudice (2002, 2004); Schaffer et al. (2013) by focusing on the slowest modes of the expansion.

VI.3 Limitations

The developed theoretical description of the network dynamics with synaptic currents is approximated, and as such it is subject to several limitations. Our starting point is a population density approach which, like others Knight et al. (1996); Brunel and Hakim (1999); Mattia and Del Giudice (2002), relies on both a diffusion and a mean-field approximation. The underlying hypotheses require that each neuron receives a large amount of spikes per unit time, and that postsynaptic currents due to single spikes induce small changes of the membrane potential. These conditions are well satisfied in neocortex Tuckwell (1988); Amit and Brunel (1997).

To reduce the dimensionality of the population density dynamics, we further resorted to two simplifications. Firstly, to safely manage the boundary conditions in the spectral expansion, we introduced a forcing white noise current with an arbitrarily small size $\sigma_{ext}$ (as in Moreno et al. (2002)). The intrinsic stochasticity of ionic channels of neuronal membrane potential can be the source of such additional input current Linaro et al. (2011); Goldwyn and Shea-Brown (2011). But, even in the limit $\sigma_{ext}\to 0$, we expect to find a good agreement between our theory and simulations. Indeed, although $\sigma_{ext}$ widely varied in the various comparisons presented in Figs. 2–5, the agreement between simulations with synaptic filtering and delay distribution did not show any significant difference. However, it is important to remark that, in the absence of this external noise, the spectral expansion we used is no more valid, as it relies on the eigenfunctions of a FP operator with $P({v_{\mathrm{thr}}},y,t)=0$. This boundary condition is not valid if $\sigma_{ext}=0$ Doering et al. (1987); Kłosek and Hagan (1998); Brunel and Sergi (1998); Schuecker et al. (2015), in which case a different basis for the probability density decomposition has to be adopted.

The second strong simplification we implemented was to consider small enough the fluctuation size $\sigma_{I}$ of the synaptic currents. Such a low-noise regime gave us the possibility to assume the marginal distribution of the current $I(t)$ narrowly distributed around its time-varying mean $\mu_{y}(t)$. As a result, the $V-I$ dynamics were reduced to a one-dimensional FP equation centered around $I(t)=\mu_{y}(t)$. Rather surprisingly, such a rough approximation for which the shape of the marginal distribution $\rho_{y}(y,t)$ is completely neglected, seems to work remarkably well. This is a result that we suspect being due to the symmetry of the narrow distribution $\rho_{y}(y,t)$, such that contributions of the terms with $I>\mu_{y}$ are almost perfectly compensated by those with $I<\mu_{y}$. The only main discrepancy we found is in the mean firing rate, which our approximated theory systematically overestimated (see Fig. 3d). This is an expected error, since synaptic filtering is known to mildly reduce the emission rate of spikes in both the small- Brunel and Sergi (1998); Moreno-Bote and Parga (2004); Schuecker et al. (2015) and the large-${\tau_{\mathrm{s}}}$ Moreno-Bote and Parga (2004) limit. Such discrepancy should disappear in our approximated development if higher-orders of the centered moment closure were incorporated. Indeed, this would allow to take into account also other features of $\rho_{y}(y,t)$, such as its standard deviation and skewness.

Acknowledgements.

Work partially supported by the University of Genoa (M.S.) and by the EU Horizon 2020 Research and Innovation Programme under HBP SGA1 (grant no. 720270 to M.M.).

VII Appendix

VII.1 Derivation of the power spectral density $P_{\nu}(\omega)$ of $\nu_{N}(t)$

In a network of a finite number $N$ of neurons, the instantaneous firing rate $\nu_{N}(t)$ display endogenous fluctuations. Indeed, due to the finite number of spikes emitted by such network in relatively small time intervals, the counting Poissonian statistics makes $\nu_{N}(t)$ well described by a stochastic variable whose variance scales as $1/N$ Spiridon and Gerstner (1999); Brunel and Hakim (1999); Mattia and Del Giudice (2002). These activity fluctuations can be seen as an ongoing stimulation of the infinite-size network and can be obtained by introducing an effective stochastic driving force different for each eigenmode of the FP operator Mattia and Del Giudice (2002). Following this approach and generalizing Eq. (11) to the case of a finite-size network, we obtain the following stochastic emission-rate equations

[TABLE]

where $\eta_{N}(t)$ models the finite-size fluctuations of the instantaneous firing rate $\nu(t)$ in the infinite-size limit. For large enough $N$, $\eta_{N}$ is well approximated by a Gaussian memoryless white noise, with zero mean and variance $\nu(t)/N$. In the above equations the coefficients depend on $\nu_{N}(t)$ instead of $\nu(t)$, as the infinitesimal mean $\mu_{I}$ and variance $\sigma^{2}_{I}$ are functions of the instantaneous spike rate of the presynaptic neurons. The additional coefficients $\vec{\psi_{0}}=\{\psi_{n}({v_{\mathrm{res}}},y)|_{y=\mu_{y}}\}$ result from having incorporated finite-size fluctuations to the boundary condition on the flux conservation of the realizations $V(t)$ exiting from ${v_{\mathrm{thr}}}$ and reentering in ${v_{\mathrm{res}}}$ Mattia and Del Giudice (2002).

In the $N\to\infty$ limit, the network dynamics can be set into an asynchronous state such that $\nu(t)=\nu^{*}$ is a fixed-point [$\nu^{*}=\Phi_{0}(\nu^{*})$] with local stability [$\Phi_{0}^{\prime}(\nu^{*})<1$] and the single-neuron spiking activity is asynchronous Amit and Brunel (1997); Brunel (2000). In this state, finite-size fluctuations bring $\nu_{N}(t)$ to wander around the fixed-point $\nu^{*}$, and Eq. (15) can be linearized by neglecting the terms of order higher than $\mathcal{O}(\eta_{N})$ in the Taylor’s series expansion around $\nu(t)=\nu^{*}$. From this linearized dynamics, the Fourier transform $\nu_{N}(\omega)$ can be obtained (see Mattia and Del Giudice (2002, 2004) for details), and the power spectral density $P_{\nu}(\omega)=|\nu_{N}(\omega)|^{2}$ turns out to be

[TABLE]

which is the same expression detailed in Eq. (14).

VII.2 Identification of mean-field parameters

To find the parameters for the numerical simulations (with NEST Gewaltig and Diesmann (2007) and the high-performance custom simulator implementing the event-based approach described in Mattia and Del Giudice (2000)), the following procedure has been adopted.

An emission rate fixed point $\nu^{*}$ was fixed a priori, and the contour line in the $(\mu,\sigma^{2})$ plane defined by $\Phi(\mu,\sigma^{2})=\nu^{*}$ was numerically determined. $\Phi(\mu,\sigma^{2})$ was computed analytically for the VIF Fusi and Mattia (1999) and LIF Ricciardi and Sacerdote (1979); Amit and Brunel (1997) neuron models. For the EIF neuron model, in the absence of an analytically expression for $\Phi(\mu,\sigma^{2})$ valid in all the considered conditions, we used a numerical cubic interpolation (using the Matlab – The MathWorks, Natick, MA – function interp2) passing through samples obtained from NEST simulation data. Depending on the regime of interest (noise- or drift-dominated), a proper point along this iso-frequency line was chosen.

Once determined $\mu$ and $\sigma^{2}$, also $J$ can be determined by imposing the value of $\frac{\mathrm{d}}{\mathrm{d}\nu}\Phi(\mu,\sigma^{2})$ at the fixed point $\nu^{*}$, which is directly related to its degree of stability Mattia and Del Giudice (2002). Indeed, by taking into account Eq. (2), it is sufficient to solve the following second-order equation in $J$:

[TABLE]

where both $\partial\Phi/\partial\mu$ and $\partial\Phi/\partial\sigma^{2}$ can be suitably computed from $\Phi(\mu,\sigma)$. Of the possible solutions, we take the one corresponding to the stable fixed point $\Phi(\mu,\sigma^{2})=\nu^{*}$ with the highest firing rate, since it is related to the proper range of $\nu$ values.

The final step is to determine mean $\mu_{ext}$ and variance $\sigma^{2}_{ext}$ of the external current $I_{ext}$ to be added to the recurrent synaptic contribution in order to obtain the chosen $\mu$ and $\sigma^{2}$. $I_{ext}$ is the sum of a constant current bias $I_{DC}$ and a Poissonian spike train from $C_{ext}$ independent sources firing at rate $\nu_{ext}$ through an instantaneous synaptic transmission with efficacy $J_{ext}$. Under diffusion approximation, $I_{ext}$ is a memoryless Wiener process and by imposing the spike rate $C_{ext}\,\nu_{ext}$ from the external neurons, the two remaining parameters $I_{DC}$ and $J_{ext}$ are uniquely determined, in the distribution of delays case, by solving the following system

[TABLE]

When the non-instantaneous synaptic transmission is incorporated, the recurrent synaptic efficacy has to be simply rescaled by the time constant ${\tau_{\mathrm{s}}}$, i.e. $J\to J/{\tau_{\mathrm{s}}}$. This is the approach used to design all the VIF, LIF and EIF neuron networks analyzed in Figs. 2-4.

For the results shown in Fig. 5, we used the following protocol. We started with the simulation of a LIF inhibitory neuron network in an asynchronous state, thus determining $J$, $I_{DC}$ and $J_{ext}$ as stated above. Then, we increased $|J|$ and kept unchanged the (stable or unstable) fixed-point $\nu=\nu^{*}$ by increasing $I_{DC}$ accordingly.

VII.3 Spiking neuron network simulations

All the simulated IF neuron networks under asynchronous state (Figs. 2-4) are modeled by considering a sparse recurrent connectivity with connection probability $C/N$, being $C$ the average number of synapses per neuron. Each neuron receives an external input from independent Poissonian spike trains with frequency $C_{ext}\,\nu_{ext}$ through instantaneous synapses with efficacy $J_{ext}$. Spikes are always delivered to post-synaptic neurons with a minimum delay $\delta_{min}$. The values for both the synaptic time constant ${\tau_{\mathrm{s}}}$ and the decay constant ${\tau_{\mathrm{d}}}$ of the exponential delay distribution are ${0,1,2,4,8,16,32,64}$ ms.

The VIF neuron networks are composed of the ‘VLSI’ integrate-and-fire neurons introduced in Fusi and Mattia (1999). This model neuron is an extended version of the standard ‘perfect integrate-and-fire’ (PIF) neuron introduced in Gerstein and Mandelbrot (1964): in addition to a constant leakage current $f(V)=\beta$, which here we consider as a part of the current bias $I_{DC}$, a reflecting barrier at $V=0$ is set to avoid a divergent diffusion towards negative membrane potentials. Contrary to the PIF neuron, this makes the VIF neuron capable of having a non-zero (positive) mean firing rate also under subthreshold regimes, i.e., for negative drifts ($\mu<0$). For the VIF neuron, input-output gain function $\Phi(\mu,\sigma)$, eigenfunctions and eigenvalues of the related FP operator have explicit analytical expressions Fusi and Mattia (1999); Mattia and Del Giudice (2002, 2004), which make this model particularly suited for matching theory and simulations. The parameters for the VIF neuron networks used in this work are listed in Tab. 1. For this model the natural unit of measure for the membrane voltage is the emission threshold ${v_{\mathrm{thr}}}$.

For the LIF Tuckwell (1988) and the EIF Fourcaud-Trocmé et al. (2003) neuron networks, the parameters at drift- (DD) and noise-dominated (ND) regimes are listed in Tabs. 2 and 3, respectively.

The LIF neuron networks for the Hopf bifurcation test (Fig. 5) are similar to the one described in Tab. 2, but with a doubled connection probability $C/N=0.1$ and an increased number of neurons, which are now inhibitory. The list of changed parameters is in Tab. 4. Only two values of synaptic time constants were considered: ${\tau_{\mathrm{s}}}=4$ ms (T4) and ${\tau_{\mathrm{s}}}=16$ ms (T16). Table 5 lists the values of $J$ and the corresponding values of $I_{DC}$ set to keep constant the fixed point $\nu^{*}$, either stable or not.

VII.4 Comparison between and characterization of simulated network activities

The time-dependent emission rate $\nu(t)$ of the simulated networks was computed as an histogram of the spike emission times of the neurons, with bin size $dt=0.5\mathrm{ms}$. The power spectral density $P_{\nu}(\omega)$ was estimated from the network emission rate by means of the standard Welch method, using a window size of 4096 samples with an overlap of 2048 samples.