# A ‘Rosetta stone’ for the population dynamics of spiking neuron networks

Gianni V. Vinci  Natl. Center for Radioprotection and Computational Physics,
Istituto Superiore di Sanità, 00169 Roma, Italy  PhD Program in Physics, “Tor
Vergata” University of Rome, 00133 Roma, Italy  Maurizio Mattia
<EMAIL_ADDRESS>Natl. Center for Radioprotection and Computational
Physics, Istituto Superiore di Sanità, 00169 Roma, Italy

(November 2, 2021 – Ver. 7MM)

###### Abstract

Populations of spiking neuron models have densities of their microscopic
variables (e.g., single-cell membrane potentials) whose evolution fully
capture the collective dynamics of biological networks, even outside
equilibrium. Despite its general applicability, the Fokker-Planck equation
governing such evolution is mainly studied within the borders of the linear
response theory, although alternative spectral expansion approaches offer some
advantages in the study of the out-of-equilibrium dynamics. This is mainly due
to the difficulty in computing the state-dependent coefficients of the
expanded system of differential equations. Here, we address this issue by
deriving analytic expressions for such coefficients by pairing perturbative
solutions of the Fokker-Planck approach with their counterparts from the
spectral expansion. A tight relationship emerges between several of these
coefficients and the Laplace transform of the inter-spike interval density
(i.e., the distribution of first-passage times). ‘Coefficients’ like the
current-to-rate gain function, the eigenvalues of the Fokker-Planck operator
and its eigenfunctions at the boundaries are derived without resorting to
integral expressions. For the leaky integrate-and-fire neurons, the coupling
terms between stationary and nonstationary modes are also worked out paving
the way to accurately characterize the critical points and the relaxation
timescales in networks of interacting populations.

##  I INTRODUCTION

As in many other complex systems like macroscopic bodies composed of atoms and
molecules [22], the collective dynamics of a network of neurons can be
effectively described looking at the time evolution of the density function
$p$ of its microscopic state variables. This density changes in time according
to the following continuity equation

$\partial_{t}p=-\vec{\nabla}\cdot\vec{S}_{p}\,,$ (1)

where $\vec{S}_{p}$ is the net flux of particles, i.e., the probability
current of neurons flowing out from an infinitesimal volume of the phase
space. Neurons interact by exchanging action potentials, the so-called
‘spikes’, which are brief and stereotyped messages traveling along their
axons. Given a ionic current $I(t)$ a neuron receives through its membrane,
the timing of the spikes it emits can be faithfully reproduced by generalized
integrate-and-fire (IF) neuron models [36, 42]. The core state variable of IF
neurons is the somatic potential $V(t)$ integrating the input $I(t)$ as a
passive RC circuit:

$\frac{dV}{dt}=F(V)+\frac{I(t)}{{C_{\mathrm{m}}}}\,.$ (2)

Here ${C_{\mathrm{m}}}$ is the membrane capacitance set to 1 by expressing
$I(t)$ in units of voltage per time, and $F(V)$ is the drifting current
determining the model-specific relaxation dynamics towards the equilibrium
potential $V=0$ [19]. $F(V)$ incorporates the decay time ${\tau_{\mathrm{m}}}$
of the free membrane potential. When $V(t)$ crosses the threshold value
${v_{\mathrm{thr}}}$, a spike is emitted and the membrane potential is reset
to $V={v_{\mathrm{res}}}$ for a refractory period ${\tau_{0}}$ past which the
subthreshold dynamics (2) is restored.

Additional state variables such as synaptic conductances or activity-dependent
ionic currents can be in principle incorporated to make more realistic the
single-neuron spiking activity [44, 19]. However, one-dimensional IF neurons
preserve without loss of generality, all the key features of the network
dynamics of interest here. In this framework the current $I(t)$ is given by
the combination
$I(t)=\sum_{i=1}^{K}J_{i}(V)\sum_{k}\delta(t-t_{i,k}-\delta_{i})$, where $K$
is the number of presynaptic neurons and $J_{i}(V)$ is the synaptic efficacy,
i.e., the instantaneous change of $V(t)$ each spike emitted at time $t_{i,k}$
by the presynaptic neuron $i$ induces after the axonal delay $\delta_{i}$. For
mammalian cortical networks it is reasonable to assume a large number of
presynaptic contacts ($K\to\infty$) and relatively small excitatory and
inhibitory synaptic efficacies ($J_{i}\to 0$) [13, 24]. In this limit, $I(t)$
can be linearized ($\partial_{V}J_{i}\simeq 0$) [3] and the diffusion limit
holds [44] such that the current is well approximated by a Gaussian white
noise with infinitesimal mean $\mu(t)$ and variance $\sigma^{2}(t)$

$\begin{array}[]{rcl}\mu(\nu)&=&KJ\nu(t)\\\
\sigma^{2}(\nu)&=&KJ^{2}(1+\Delta^{2})\nu(t)\end{array}\,.$ (3)

For simplicity’s sake, here we considered a homogeneous population of
presynaptic neurons with firing rate
$\nu(t)=1/K\sum_{i=1}^{K}\sum_{k}\delta(t-t_{i,k}-\delta_{i})$ and synaptic
efficacies randomly distributed with mean $J$ and standard deviation $\Delta
J$.

By incorporating the above stochastic current in Eq. (2), a Langevin equation
for the membrane potential of a single neuron results:

$dV=\left[F(V)+\mu\right]dt+\sigma dW\,.$ (4)

Here, $W(t)$ is the standard memoryless Wiener process with 0 mean and unit
variance: $\langle W(t)W(t^{\prime})\rangle=\delta(t-t^{\prime})$. Note that,
as the infinitesimal moments in Eq. (3) depend on time via the firing rate
$\nu(t)$, $V(t)$ is in general a inhomogeneous stochastic process [7, 26]
leading to a faithful description of the out-of-equilibrium population
dynamics. In homogeneous populations of $N$ neurons with similar single-neuron
parameters like ${v_{\mathrm{thr}}}$, ${v_{\mathrm{res}}}$ and
${\tau_{\mathrm{m}}}$, and constant connection likelihood $\epsilon=K/N$, the
membrane potentials of different cells can be seen as independent stochastic
realizations of the same Eq. (4). In the cortical network limit such
independence is granted by having $J\to 0$, while the central limit theorem
with $K\to\infty$ leads to have input currents with the same $\mu$ and
$\sigma^{2}$.

This ‘extended’ mean-field approximation [2] including also the fluctuation
size $\sigma$ of the ‘field’ (i.e., the synaptic current), is the foundation
of the so-called population density approach [21, 7, 31, 20, 26] where the
continuity equation (1) reduces to the following Fokker-Planck equation:

$\partial_{t}p=-\partial_{v}S_{p}=-\partial_{v}[(F+\mu)p]+\frac{1}{2}\sigma^{2}\partial^{2}_{v}p\equiv\mathcal{L}\,p\,.$
(5)

Here $\mathcal{L}=-\partial_{v}(F+\mu)+\frac{1}{2}\sigma^{2}\partial^{2}_{v}$
is the Fokker-Planck operator. From Eq. (3) it is apparent that all neurons
are driven by different currents with same statistics determined by the firing
rate $\nu(t)$ of the whole population. Thus, although two generic neurons are
not directly influencing each other, they both have dynamics conditioned by
the same collective activity of the network.

In this framework, the firing rate is the flux of realizations (i.e., neurons)
crossing the emission threshold ${v_{\mathrm{thr}}}$

$\nu(t)=S_{p}({v_{\mathrm{thr}}})=-\frac{1}{2}\left.\sigma^{2}\partial^{2}_{v}p\right|_{v={v_{\mathrm{thr}}}}\,,$
(6)

resulting from Eq. (5) by setting $p({v_{\mathrm{thr}}},t)=0$ as
${v_{\mathrm{thr}}}$ is an absorbing barrier for the diffusive process $V(t)$
[44].

Finding the general solutions of Eq. (5) is challenging, and this is due to
two main reasons. Firstly, the Fokker-Planck operator depends on $\mu(\nu)$
and $\sigma^{2}(\nu)$, and, as the firing rate depends on $p(v,t)$ from Eq.
(6), it makes $\mathcal{L}=\mathcal{L}(p)$. Eq. (5) is thus a nonlinear
partial differential equation. Secondly, when a neuron emits a spike, it does
not disappear even if it is ‘absorbed’ in ${v_{\mathrm{thr}}}$. Actually, its
dynamics is restored after resetting the membrane potential to
${v_{\mathrm{res}}}$. This guarantees that the number of neurons is conserved,
i.e., $\int_{v_{\mathrm{min}}}^{v_{\mathrm{thr}}}{p(v,t)dv}=1$ at any time
$t$. Such reset is incorporated in Eq. (5) by reentering the exiting flux of
realizations in ${v_{\mathrm{thr}}}$ as a point-like source
$\delta(v-{v_{\mathrm{res}}})\nu(t-{\tau_{0}})$. This source gives rise to a
peculiar boundary condition:
$S_{p}({v_{\mathrm{res}}}^{+})-S_{p}({v_{\mathrm{res}}}^{-})=\nu(t-{\tau_{0}})$.
As a result, the flux $S_{p}(v)$ is non-vanishing even under stationary
conditions [1, 21, 7, 16], leading to a breaking of the so-called ‘potential
condition’ [39] usually holding in equilibrium statistical physics and hardly
relying on the detailed balance.

In addition to the numerical integration of the partial differential equations
(1) and (5) [31, 33, 12, 8], their analytical characterization is mainly
investigated via the linear response theory dealing with the small deviations
from the equilibrium states [1, 43, 7, 16, 20, 6, 37, 38]. Due to its
derivation intricacy only few valuable attempts have been pursued
investigating perturbatively nonlinear dynamical regimes like limit cycles [7,
5].

An alternative approach is to expand the density $p$ as a linear combination
of the eigenfunctions $\phi_{n}(v)$ of the Fokker-Planck operator
$\mathcal{L}$

$\mathcal{L}\phi_{n}=\lambda_{n}\phi_{n}\,,$ (7)

corresponding to the eigenvalues $\lambda_{n}$, i.e., the spectrum of
$\mathcal{L}$. Firstly adopted in [1, 21] to study the population dynamics of
neuronal networks, this is a standard approach in statistical physics [39,
15]. It is reminiscent of the time-dependent Hartree-Fock theory in quantum
mechanics [29] where the basis $\\{\phi_{n}\\}_{n\in\mathbb{Z}}$ moves in time
following the system evolution. Indeed, $\mathcal{L}$ – and hence $\phi_{n}$ –
depends parametrically on the time-dependent firing rate [21, 26].
Interestingly, the nonstationary modes $\phi_{n}$ associated to the non-
vanishing eigenvalues $\lambda_{n}$ of $\mathcal{L}$, display a hierarchy of
timescales [21, 26, 4] such that only the slowest modes of the spectrum
contribute to the firing rate dynamics making it low-dimensional [1, 26, 40].

This elegant ‘spectral expansion’ approach is particularly well suited to
describe the out-of-equilibrium dynamics of neuronal networks as it does not
relies on the perturbations of equilibrium states. However, it comes at a
price: the coefficients of such expansion taking into account recurrent and
external synaptic couplings, involve integrals that only in simplified models
can be analytically solved [20, 26]. Not only, if the firing rate fluctuates
due to a broadband input from upstream neurons – or endogenously expressed in
finite-size networks (i.e., composed of a finite number of neurons) – are
taken into account, the approach requires to manage series which are not
easily summable. This is an issue typically arising when second-order
statistics like the Fourier power-spectrum $P_{\nu}(\omega)=|\nu(\omega)|^{2}$
is studied under stationary conditions [26, 27, 25].

Starting from such remarks, here we address the following question: is there
any way to exploit the advantages of both the perturbative and the spectral
expansion approach to overcome their intrinsic limitations? In what follows we
aim at using the relaxation dynamics of an uncoupled set of neurons and the
response to small perturbations of a coupled network as a ‘Rosetta stone’.
This will allow us to translate some key expressions with a closed form in
both linear response and renewal theory into compact and manageable sum of
series and state-dependent coefficients for the spectral expansion approach.
The results of this effort pave the way to a further exploitation of the
spectral expansion to investigate the out-of-equilibrium dynamics of spiking
neuron networks.

##  II Matching the relaxation of $\nu(t)$

Under stationary conditions and in presence of synaptic current fluctuations
($\sigma\neq 0$), the collective firing rate $\nu(t)$ of a set of uncoupled IF
neurons always approaches an equilibrium point $\nu_{0}$. Due to the reset
mechanism following the emission of a spike, IF neurons are Markovian renewal
processes. As such, inter-spike intervals (ISIs) separating two consecutive
spikes do not depend on the ISIs occurred before. Finding the probability
$\rho(t)dt$ to have an ISI in the interval $[t,t+dt]$ is a first-passage time
problem [44], since it requires to know when $V(t)$ crosses the threshold
${v_{\mathrm{thr}}}$ for the first time starting from
$V(0)={v_{\mathrm{res}}}$. With this initial condition, the ISI density can be
interpreted as the probability of first-spike occurrence
$\rho_{1}(t)\equiv\rho(t)$, leading to recursively define the density
$\rho_{k}(t)$ of the $k$-th emitted spike as

$\rho_{k}(t)=\int_{0}^{t}\rho(t)\rho_{k-1}(t-\tau)d\tau\,.$ (8)

This ‘renewal equation’ [44, 19] allows to derive the so-called ‘hazard’
function, i.e., the mean density of spikes emitted by an isolated IF neuron at
time $t$

$\nu(t)=\sum_{k=1}^{\infty}\rho_{k}(t)\,,$ (9)

which in turn is equivalent to the relaxation to equilibrium of the firing
rate $\nu(t)$.

The firing rate in Eq. (9) has a straightforward Laplace trasform
$\nu(s)\equiv\int_{0}^{\infty}\nu(t)e^{-st}dt$. Indeed, being $\rho_{k}(t)$ a
convolution, Eq. (8) reduces to

$\rho_{k}(s)=\rho(s)\rho_{k-1}(s)=\rho(s)^{k}\,.$

Introducing it in Eq. (9), $\nu(s)$ results to be

$\begin{split}\nu(s)&=\sum_{k=1}^{\infty}\rho_{k}(s)=\sum_{k=1}^{\infty}\rho(s)^{k}=\frac{1}{1-\rho(s)}-1\\\
&=\frac{\rho(s)}{1-\rho(s)}\,.\end{split}$ (10)

This equation establishes a direct nonlocal relationship between the firing
rate of a homogeneous pool of independent neurons and the ISI density of a
single cell [10, 19].

###  II.1 Relaxation dynamics from the spectral expansion

The firing rate $\nu(t)$ of the mentioned set of independent neurons can be
alternatively derived by expanding the density $p(v,t)$ in Eq. (5) as

$p(v,t)=\sum_{n}a_{n}(t)\phi_{n}(v)\,,$

where
$a_{n}(t)=\langle\psi_{n}|p\rangle=\int_{v_{\mathrm{min}}}^{v_{\mathrm{thr}}}\psi_{n}(v)p(v,t)dv$
are the projections of $p$ on the eigenfunctions (i.e., the modes) of the
operator $\mathcal{L}$ defined in Eq. (7). The infinite set
$\\{\psi_{n}(v)\\}$ with $n\in\mathbb{Z}$, is an orthonormal basis of a non-
Hermitian space such that $\langle\psi_{n}|\phi_{m}\rangle=\delta_{nm}$ [21,
26]. Here, $\psi_{n}(v)$ are the eigenfunctions of the adjoint operator
$\mathcal{L}^{\dagger}=(F+\mu)\partial_{v}+\frac{1}{2}\sigma^{2}\partial^{2}_{v}$
defined as
$\langle\mathcal{L}^{\dagger}\psi_{n}|\phi_{n}\rangle\equiv\langle\psi_{n}|\mathcal{L}\phi_{n}\rangle$,
and having the same spectrum as in Eq. (7):
$\mathcal{L}^{\dagger}\psi_{n}=\lambda_{n}\psi_{n}$. The ‘expanded’ $p$ in the
Fokker-Planck equation leads to the following equivalent dynamics for the
projections $a_{n}$:

$\begin{split}\dot{\vec{a}}&=\mathbf{\Lambda}\vec{a}\\\
\nu&=\Phi+\vec{f}\cdot\vec{a}\end{split}\,,$ (11)

where the contribution due to the stationary mode ($n=0$ with $\lambda_{0}=0$)
has been isolated, and a matrix formalism has been adopted such that
$\\{\vec{a}\\}_{n}=a_{n}$ and
$\\{\mathbf{\Lambda}\\}_{nm}=\lambda_{n}\delta_{nm}$ with $m,n\neq 0$ [21,
26]. The current-to-rate gain function

$\Phi(\mu,\sigma)=S_{\phi_{0}}({v_{\mathrm{thr}}})=-\left.\frac{1}{2}\sigma^{2}\partial_{v}^{2}\phi_{0}\right|_{v={v_{\mathrm{thr}}}}$
(12)

is the flux of realization crossing ${v_{\mathrm{thr}}}$ under stationary
condition, i.e., the asymptotic firing rate $\nu_{0}$. The flux due to the
nonstationary modes ($\phi_{n}$ with $n\neq 0$) are instead the elements of
the infinite vector $\vec{f}$, which can be conveniently set to
$f_{n}=1/{\tau_{\mathrm{m}}}$ [14].

The linear system (11) has a straightforward solution having constant
coefficients ($\mu$ and $\sigma$ do not depend on time):

$\begin{split}\vec{a}(t)&=e^{\mathbf{\Lambda}t}\vec{a}(0)\\\
\nu(t)&=\Phi+\vec{f}\cdot\vec{a}\end{split}\,.$ (13)

As all neurons have a starting membrane potential $V(0)={v_{\mathrm{res}}}$,
$p(v,0)=\delta(v-{v_{\mathrm{res}}})$ and the initial value of the projections
result to be

$a_{n}(0)=\langle\psi_{n}|p(v,0)\rangle=\psi_{n}({v_{\mathrm{res}}})\,.$

###  II.2 ISI moments and sums in the spectral expansion

A first sum of coefficients from the spectral expansion results by setting
$t=0$ in the relaxation dynamics (13):

$\Phi=-\vec{f}\cdot{\vec{\psi}_{\mathrm{res}}}\equiv-\frac{1}{{\tau_{\mathrm{m}}}}\sum_{n\neq
0}{\psi_{n}({v_{\mathrm{res}}})}\,,$ (14)

where we took into account that $\nu(0)=0$. To our knowledge, this is a novel
expression to carry out the current-to-rate gain function $\Phi(\mu,\sigma)$
for any IF neuron model.

We anticipate that this is only a particular case of a more general result.
Indeed, the Laplace transform $\nu(s)$ from the renewal theory in Eq. (10) can
be directly compared to the one resulting from the firing rate equation (11).
In this case we have

$\begin{split}s\,\vec{a}(s)-\vec{a}(0)&=\mathbf{\Lambda}a(s)\\\
\nu(s)&=\frac{\Phi}{s}+\vec{f}\cdot\vec{a(s)}\end{split}\,,$

leading to
$\vec{a}(s)=(s\mathbf{I}-\mathbf{\Lambda})^{-1}{\vec{\psi}_{\mathrm{res}}}$,
and eventually to

$\nu(s)=\frac{\Phi}{s}+\vec{f}\cdot(s\mathbf{I}-\mathbf{\Lambda})^{-1}{\vec{\psi}_{\mathrm{res}}}\,.$
(15)

Recalling Eq. (10) and comparing the two transform $\nu(s)$ we obtain:

$\frac{\rho(s)}{1-\rho(s)}-\frac{\Phi}{s}=\vec{f}\cdot(s\mathbf{I}-\mathbf{\Lambda})^{-1}{\vec{\psi}_{\mathrm{res}}}\equiv
h(s)\,,$ (16)

highlighting a tight relationship between the ISI density $\rho(t)$, the
eigenvalues $\lambda_{n}$ and ${\vec{\psi}_{\mathrm{res}}}$.

From Eq. (16) other sums can be worked out recalling that the ISI moments
${\langle(-t)^{k}\rangle}=\lim_{s\to 0}\frac{d^{k}\rho(s)}{ds^{k}}$. We then
perform the same limit on both hand sides of the equation finding that

$\lim_{s\to 0}h(s)=\frac{\Phi{\langle t^{2}\rangle}-2{\langle
t\rangle}}{2{\langle t\rangle}}\,.$

Here, considering that the gain function $\Phi=1/{\langle t\rangle}$ and that
$h(0)=-\vec{f}\cdot\mathbf{\Lambda}^{-1}{\vec{\psi}_{\mathrm{res}}}$, we
obtain

$h(0)=\frac{c_{v}^{2}-1}{2}=-\frac{1}{{\tau_{\mathrm{m}}}}\sum_{n\neq
0}\frac{\psi_{n}({v_{\mathrm{res}}})}{\lambda_{n}}\,,$ (17)

where $c_{v}=\sqrt{{\langle t^{2}\rangle}-{\langle t\rangle}^{2}}/{\langle
t\rangle}$ is the coefficient of variation of the ISIs.

We remark that other sums can be worked out deriving by $s$ both hand sides of
Eq. (16) and taking the limit $s\to 0$. Indeed, in this limit the derivatives
of $h(s)$ reduce to

$\left.\frac{d^{k}h}{ds^{k}}\right|_{s\to
0}=-k!\vec{f}\cdot\mathbf{\Lambda}^{-(k+1)}{\vec{\psi}_{\mathrm{res}}}=-\frac{k!}{{\tau_{\mathrm{m}}}}\sum_{n\neq
0}\frac{\psi_{n}({v_{\mathrm{res}}})}{\lambda_{n}^{k+1}}$

for any $k>0$.

###  II.3 Eigenvalues $\lambda_{k}$ and eigenfunction values
$\psi_{k}({v_{\mathrm{res}}})$

We can further exploit Eq. (16) by making use of the Cauchy’s residue theorem.
Indeed,

$\mathrm{Res}_{s=\lambda_{k}}h(s)=\mathrm{Res}_{s=\lambda_{k}}\frac{\rho(s)}{1-\rho(s)}=\psi_{k}({v_{\mathrm{res}}})$
(18)

for $\lambda_{k}\neq 0$. This allows us to work out from the Laplace transform
of the ISI probability density the value of the $k$-th eigenfunction of
$\mathcal{L}^{\dagger}$ at $v={v_{\mathrm{res}}}$. In this way
$\psi_{k}({v_{\mathrm{res}}})$ can be found without explicitly knowing neither
the analytic expression for $\psi(v,s)$ nor the exact value of the
$\lambda_{k}$. As $s=\lambda_{k}$ are poles of $h(s)$, a generic line integral
of $\rho/(1-\rho)$ around each eigenvalues will give
$\psi_{k}({v_{\mathrm{res}}})$. Similarly, the eigenvalues $\lambda_{k}$ can
be obtained without resorting to a often complicated minimization procedure as

$\frac{\mathrm{Res}_{s=\lambda_{k}}s\,h(s)}{\mathrm{Res}_{s=\lambda_{k}}h(s)}=\frac{\lambda_{k}\psi_{k}({v_{\mathrm{res}}})}{\psi_{k}({v_{\mathrm{res}}})}=\lambda_{k}\,.$
(19)

###  II.4 Spectral equation for $\lambda_{k}$ from $\rho(s)$

The link between the spectrum of $\mathcal{L}^{\dagger}$ and $\rho(s)$ pointed
out by Eq. (16) should not be surprising as $\rho(t)$ is the density of a
first-passage time. As such, $\rho(t)$ conditioned to have
$V(0)={v_{\mathrm{res}}}$ and $V(t)={v_{\mathrm{thr}}}$, solves the backward
Kolmogorov equation [10, 44]

$\partial_{t}\rho=A({v_{\mathrm{res}}})\partial_{v_{\mathrm{res}}}\rho+\frac{1}{2}B({v_{\mathrm{res}}})\partial_{v_{\mathrm{res}}}^{2}\rho=\mathcal{L}^{\dagger}_{v_{\mathrm{res}}}\rho\,,$

with boundary condition
$\left.\rho(t)\right|_{{v_{\mathrm{res}}}={v_{\mathrm{thr}}}}=\delta(t)$. In
our case of study
$A({v_{\mathrm{res}}})=\left.F+\mu\right|_{v={v_{\mathrm{res}}}}$ and
$B({v_{\mathrm{res}}})=\sigma^{2}$. For the sake of clarity, we write
explicitly the dependence of the operator
$\mathcal{L}^{\dagger}_{v_{\mathrm{res}}}$ on the initial value of the
membrane potential.

Performing the Laplace transform of both hand sides of the above equation we
obtain

$\left.-\rho(t)\right|_{t=0}+s\rho(s)=\mathcal{L}^{\dagger}_{v_{\mathrm{res}}}\rho(s)$

with $\rho(s)$ the Laplace transform of the ISI density implicitly depending
on ${v_{\mathrm{thr}}}$. In the case of interest,
${v_{\mathrm{res}}}<{v_{\mathrm{thr}}}$ and no spikes are emitted at $t=0$
such that $\left.\rho(t)\right|_{t=0}=0$ and the above equation reduces to

$(\mathcal{L}^{\dagger}_{v_{\mathrm{res}}}-s)\rho({v_{\mathrm{res}}},s)=0\,,$

where we made explicit the dependence of $\rho$ on the initial potential
${v_{\mathrm{res}}}$. As pointed out in [41], this is exactly the same
equation for the eigenfunctions $\psi_{k}(v)=\psi(v,\lambda_{k})$ provided
that $s=\lambda_{k}$ leading to the equivalence:

$\rho({v_{\mathrm{res}}},s)=a(s)\psi({v_{\mathrm{res}}},s)\,.$

The arbitrary function $a(s)$ can be derived taking into account the boundary
condition $\rho({v_{\mathrm{thr}}},t)=\delta(t)$ leading to
$\rho({v_{\mathrm{thr}}},s)=1$ and hence to
$a(s)=1/\psi({v_{\mathrm{thr}}},s)$, allowing to recover the known
relationship [41]

$\rho(s)=\frac{\psi({v_{\mathrm{res}}},s)}{\psi({v_{\mathrm{thr}}},s)}\,.$
(20)

Following [35], this equation allows to derive a spectral equation to
determine the eigenvalues $\lambda_{k}$. Indeed, from the boundary condition
associated to the conservation of flux of realizations exiting from
${v_{\mathrm{thr}}}$ and reentering in ${v_{\mathrm{res}}}$ after the emission
of spike, we have
$\psi({v_{\mathrm{res}}},\lambda_{k})=\psi({v_{\mathrm{thr}}},\lambda_{k})$
[21, 20, 26]. Incorporating this constraint into Eq. (20) we eventually have

$\rho(\lambda_{k})=1\,,$ (21)

i.e., yet another way to compute the eigenvalues of the Fokker-Planck operator
in Eq. (7).

###  II.5 Normalization factor of eigenfunctions $\psi_{k}(v)$

As summarized at the beginning of this Section, the spectral expansion of the
population density relies on the assumption that the eigenfunctions
$\psi_{k}(v)$ and $\phi_{k}(v)$ are orthonormal, i.e.,

$\langle\psi_{n}|\phi_{m}\rangle=\int_{v_{\mathrm{min}}}^{v_{\mathrm{thr}}}{\psi_{n}(v)\phi_{m}(v)dv}=\delta_{nm}\,.$

Such condition is guaranteed by normalizing the generic solution of the
equation $(\mathcal{L}^{\dagger}-\lambda_{k})\tilde{\psi}_{k}=0$ by a factor
$Z_{k}$ in order to have $\psi_{k}\equiv\tilde{\psi}_{k}/Z_{k}$ and

$Z_{k}=\langle\tilde{\psi_{k}}|\phi_{k}\rangle=\int_{v_{\mathrm{min}}}^{v_{\mathrm{thr}}}{\tilde{\psi}_{k}(v)\phi_{k}(v)dv}\,.$

This integral expression requires to know explicitly the eigenfunctions of
both $\mathcal{L}$ and $\mathcal{L}^{\dagger}$. Not only, even when they are
available, a closed form of this integral is difficult to be obtained. Hence,
its numerical evaluation is required even when the integration domain is
infinite, as usually ${v_{\mathrm{min}}}\to\infty$, eventually making the
computation of $Z_{k}$ demanding.

Even in this case, a way out to this issue is offered by the residue theorem
starting from the knowledge of $\rho(s)$. Indeed, inserting Eq. (21) into Eq.
(18) we have

$\mathrm{Res}_{s=\lambda_{k}}\frac{\tilde{\psi}({v_{\mathrm{res}}},s)}{\tilde{\psi}({v_{\mathrm{thr}}},s)-\tilde{\psi}({v_{\mathrm{res}}},s)}=\psi_{k}({v_{\mathrm{res}}})=\frac{\tilde{\psi}_{k}({v_{\mathrm{res}}})}{Z_{k}}\,,$

where we took into account that the normalization factor does not depend on
$v$ and thus the equivalence
$\psi({v_{\mathrm{res}}},s)/\psi({v_{\mathrm{thr}}},s)=\tilde{\psi}({v_{\mathrm{res}}},s)/\tilde{\psi}({v_{\mathrm{thr}}},s)$
holds. Putting out from the residue computation the numerator which becomes
$\tilde{\psi}({v_{\mathrm{res}}},s)\to\tilde{\psi}_{k}({v_{\mathrm{res}}})$,
the following expression results

$\mathrm{Res}_{s=\lambda_{k}}\frac{1}{\tilde{\psi}({v_{\mathrm{thr}}},s)-\tilde{\psi}({v_{\mathrm{res}}},s)}=\frac{1}{Z_{k}}$
(22)

giving the normalization factor as a limit expression without both computing
any integral and knowing $\phi_{k}(v)$.

###  II.6 ISI density transform for notable neuron models

The simplest model is the perfect IF (PIF) neuron introduced in [18] in which
the membrane potential $V(t)$ is a diffusion Wiener process with drift, i.e.,
$F(V)=\mu$. In this model ISI moments are finite only if $\mu>0$, and the ISI
density $\rho(s)$ is an inverse Gaussian whose Laplace transform is [11, 10,
44]

$\rho_{\mathrm{PIF}}(s)=\mathrm{exp}\left[\frac{{v_{\mathrm{thr}}}-{v_{\mathrm{res}}}}{\sigma^{2}}\left(\mu-\sqrt{\mu^{2}+2\sigma^{2}s}\right)\right]\,.$
(23)

From this, the first and second moments of the ISI result to be ${\langle
t\rangle}=1/\Phi=({v_{\mathrm{thr}}}-{v_{\mathrm{res}}})/\mu$ and ${\langle
t^{2}\rangle}=\sigma^{2}/(\mu^{2}\Phi)$ [44].

A generalization of the PIF neuron model is the one in which a reflecting
barrier at ${v_{\mathrm{min}}}=0$ is incorporated, allowing to have finite
moments of the ISI also for $\mu\leq 0$. The model was introduced as a former
‘neuromorphic’ implementation in VLSI circuits of an IF neuron [30]. This VLSI
IF (VIF) neuron has a closed form of $\rho(s)$ [16]:

$\rho_{\mathrm{VIF}}(s)=\frac{\zeta(s)e^{\xi}}{\zeta(s)\cosh[\zeta(s)]+\xi\sinh[\zeta(s)]}$
(24)

with $\xi={v_{\mathrm{thr}}}\mu/\sigma^{2}$ and
$\zeta(s)=\sqrt{\xi^{2}+2s{v_{\mathrm{thr}}}/\sigma^{2}}$, fixing the reset
potential at the reflecting barrier ${v_{\mathrm{res}}}={v_{\mathrm{min}}}=0$.

The standard, and widely used, leaky integrate-and-fire (LIF) neuron with
$F(V)=-V/{\tau_{\mathrm{m}}}$ has also a closed form for the Laplace transform
of the ISI density [41, 9]:

$\rho_{\mathrm{LIF}}(s)=\frac{\sqrt{e^{{x_{\mathrm{r}}}^{2}}}\mathcal{D}_{-s}(-\sqrt{2}{x_{\mathrm{r}}})}{\sqrt{e^{{x_{\mathrm{t}}}^{2}}}\mathcal{D}_{-s}(-\sqrt{2}{x_{\mathrm{t}}})}$
(25)

with
${x_{\mathrm{t}}}=({v_{\mathrm{thr}}}-\mu{\tau_{\mathrm{m}}})/\sigma\sqrt{{\tau_{\mathrm{m}}}}$
and
${x_{\mathrm{r}}}=({v_{\mathrm{res}}}-\mu{\tau_{\mathrm{m}}})/\sigma\sqrt{{\tau_{\mathrm{m}}}}$.
$\mathcal{D}_{\nu}(z)$ is the parabolic cylinder function solution of the
Weber differential equation [32].

From Eq. (21), the eigenvalues for the PIF neuron result to be [35]

$\lambda_{n}=-2\pi^{2}\frac{\sigma^{2}}{({v_{\mathrm{thr}}}-{v_{\mathrm{res}}})^{2}}n^{2}+i2\pi\Phi
n\,.$ (26)

As under strongly drift-dominated regime ($\mu\gg F({v_{\mathrm{thr}}})-\mu$)
both LIF and VIF neurons are well approximated by a PIF model, it is not
surprising to see that the above expression has been previously worked out in
both [1, 21, 20] and [26], respectively. Equation (26) can then be exploited
as an initial guess of $\lambda_{n}$ in order to identify a suited line
integral to work out the actual eigevalues from Eq. (19). Note that, in the
regimes dominated by the current fluctuations $\sigma$, $\lambda_{n}$ are no
more well approximated by Eq. (26) as they are real values [26, 45].

##  III Linear response as a ’Rosetta stone’

###  III.1 Linear response from spectral expansion

Under the ‘extended’ mean-field approximation [2], the drift and diffusion
coefficients of the Fokker-Planck equation for a network of coupled spiking
neurons are function of the firing rate $\nu$. This because both the
infinitesimal mean and variance of the synaptic current depend on the network
activity: $\mu=\mu(\nu)$ and $\sigma^{2}=\sigma^{2}(\nu)$. In this general
framework the dynamics of the projection coefficients $a_{k}(t)$ and of
$\nu(t)$ result to be [21, 26]

$\begin{split}\dot{\vec{a}}&=\left(\mathbf{\Lambda}+\mathbf{C}\dot{\nu}\right)\vec{a}+\vec{c}\,\dot{\nu}\\\
\nu&=\Phi+\vec{f}\cdot\vec{a}\end{split}\,,$ (27)

where in addition to the coefficients found in Eq. (11), the coupling terms
between modes are defined as
$\\{\mathbf{C}\\}_{nm}=\langle\partial_{\nu}\psi_{n}|\phi_{m}\rangle$ and
$\\{\vec{c}\\}_{n}=\langle\partial_{\nu}\psi_{n}|\phi_{0}\rangle$ for any
$m,n\neq 0$.

The couplings can in general affect both current moments $\mu$ and
$\sigma^{2}$ independently, and perturbations of the input received by the
neurons in the network can be a combination of the following changes

$\begin{split}\mu(t)&=\mu_{0}+\epsilon(t)\\\
\sigma^{2}(t)&=\sigma^{2}_{0}+\epsilon(t)\end{split}\,,$

In what follows, only the input changes $\epsilon(t)$ up to the first order
are taken into account (i.e., $o(\epsilon^{2})\simeq 0$).

Under this approximation, the Fourier transform of Eq. (27) leads to the
Fourier transfer function $H_{\epsilon}(\omega)$ of the firing rate [26, 27]

$H_{\epsilon}(\omega)\equiv\frac{\nu(\omega)}{\epsilon(\omega)}=\partial_{\epsilon}\Phi+i\omega\vec{f}\dot{(}i\omega\mathbf{I}-\mathbf{\Lambda})^{-1}\vec{c}_{\epsilon}\,,$
(28)

where $\nu(\omega)=\int_{-\infty}^{\infty}{\nu(t)e^{-i\omega t}dt}$ is the
Fourier transform of the firing rate fluctuating around
$\nu_{0}=\Phi(\mu_{0},\sigma_{0})$, $\epsilon(\omega)$ is the same transform
for $\epsilon(t)$ and
$\\{\vec{c}_{\epsilon}\\}_{n}=\langle\partial_{\epsilon}\psi_{n}|\phi_{0}\rangle$.
With this notation we generalize the derivatives associated to the input
changes as $d/d\epsilon$, which eventually should be replaced by either
$d/d\mu$ or $d/d\sigma^{2}$.

###  III.2 Eigenfunction $\psi_{n}$ from Fokker-Planck equation

The same transfer function $H(\omega)$ can be worked out directly from the
Fokker-Planck Eq. (5) [17], a perturbative approach adopted in the study of
spiking neuron networks (see for instance [6, 23]). In networks of LIF neurons
for which analytical results can be carried out, we start rescaling the
membrane potential $V$ as

$x=\frac{V-\mu_{0}}{\sigma_{0}}\,,$ (29)

leading to simplify the Fokker-Planck equation (5) into:

$\partial_{t}p(x,t)=\partial_{x}(xp)+\frac{1}{2}\partial_{x}^{2}p\equiv\mathcal{L}_{x}p\,.$

As above, $\mu_{0}$ and $\sigma_{0}$ do not depend on time and set the fixed
point around which perturbations are studied. For the sake of simplicity we
set ${\tau_{\mathrm{m}}}=1$. A different unit time can be adopted simply
multiplying by ${\tau_{\mathrm{m}}}$ all the variables like $\mu_{0}$,
$\sigma_{0}^{2}$, $\nu(t)$ and $\lambda_{n}$.

In general, the spectral equation $(\mathcal{L}_{x}-\lambda_{n})\phi_{n}=0$ is
a second-order differential equation which can be recast into the general form

$k_{2}(x)\frac{d^{2}\phi_{n}}{dx^{2}}+k_{1}(x)\frac{d\phi_{n}}{dx}+k_{0}(\lambda_{n})\phi_{n}=0$
(30)

where $\phi_{n}$ result to be a combination of two fundamental (i.e.,
independent) solutions $f_{1}(x,\lambda_{n})$ and $f_{2}(x,\lambda_{n})$ of
the same equation. The coefficients $a$, $b$ and $d$ of the combination are
determined by solving a linear system determined by the boundary conditions of
the problem (see for instance [14]):

$\phi_{n}(x)=\left\\{\begin{matrix}a\,f_{1}(x,\lambda)+b\,f_{2}(x,\lambda)&{x_{\mathrm{r}}}<x<{x_{\mathrm{t}}}\\\
d\,f_{2}(x,\lambda)&x<{x_{\mathrm{r}}}\end{matrix}\right.\,,$

where we arbitrarily impose $f_{2}(x,\lambda_{n})\to 0$ in the limit
$x\to{x_{\mathrm{min}}}$. The threshold $x_{\mathrm{t}}$ and the reset
${x_{\mathrm{r}}}$ values comes from Eq. (29) by setting
$V={v_{\mathrm{thr}}}$ and $V={v_{\mathrm{res}}}$, respectively.

The eigenfunctions $\psi_{n}(x)=\psi(x,\lambda_{n})$ of the adjoint Fokker-
Planck operator $\mathcal{L}^{\dagger}_{x}$ can be derived from the same
fundamental solutions of Eq. (30). They rely also on the Wronskian
$\mathcal{W}(x)=f_{1}^{\prime}(x)f_{2}(x)-f_{1}(x)f^{\prime}_{2}(x)$ which can
be in general expressed as [32]

$\mathcal{W}(x)=C\,e^{-\displaystyle\int\frac{k_{1}(x)}{k_{2}(x)}dx}\,,$ (31)

being a solution of the differential equation
$\mathcal{W}^{\prime}(x)=-k_{1}(x)/k_{2}(x)$ solvable up to a constant factor
$C$. Note that although $\mathcal{W}$ depends on $f_{1}$ and $f_{2}$ it does
not depend on $\lambda_{n}$ and for the LIF neuron we have
$\mathcal{W}_{\mathrm{LIF}}(x)=e^{-x^{2}}$.

As mentioned above, to derive the adjoint eigenfunctions two boundary
conditions must be taken into account: i)
$\psi_{n}({x_{\mathrm{r}}})=\psi_{n}({x_{\mathrm{r}}})$ and ii) $\psi_{n}(x)$
must be smooth on all the domain of $x$ [21, 20, 26]. With these conditions
the adjoint eigenfunctions are in general

$\psi_{n}(x)=\psi(x,\lambda_{n})=\frac{f_{2}(x,\lambda_{n})}{Z_{n}\,\mathcal{W}(x)}\,,$
(32)

being it a solution of the spectral equation
$(\mathcal{L}^{\dagger}_{x}-\lambda_{n})\psi_{n}=0$ if Eq. (31) is also taken
into account [14]. The normalization factor $Z_{n}$ is the same as the Eq.
(22).

###  III.3 Linear response from Fokker-Planck equation

To derive a closed formula for the transfer function $H_{\epsilon}(\omega)$
directly from the Fokker-Planck equation we follow the perturbative approach
exploited in [7, 37] to the study of IF neuron networks. In this framework Eq.
(5) with the standardized variable $x$ can be decomposed by developing in
Taylor series the probability density

$p(x,t)=\phi_{0}(x)+\epsilon(t)p_{1}(x)+o(\epsilon^{2})\,.$

The eigenfunction $\phi_{0}$ is the stationary probability density when the
input current has moment $\mu_{0}$ and $\sigma_{0}$. As a consequence, the
firing rate $\nu(t)$ is

$\nu(t)=\Phi+\frac{1}{2\sigma_{0}}p_{1}^{\prime}(x_{\mathrm{t}})\epsilon(t)+o(\epsilon^{2})\,,$

from which the transfer function we are looking for results to be

$H_{\epsilon}^{(P)}(\omega)=\frac{\nu(\omega)}{\epsilon(\omega)}=\frac{p_{1}^{\prime}(x_{\mathrm{t}})}{2\sigma_{0}}\,.$
(33)

Note that the factor $1/\sigma_{0}$ comes from having adopted $x$ as state
variable, and due to this $\partial_{v}=1/\sigma_{0}\partial_{x}$.

Considering that also the operator $\mathcal{L}_{x}$ depends on $\epsilon$ we
can write it as
$\mathcal{L}_{x}=\mathcal{L}_{x0}+\epsilon(t)\mathcal{L}_{x1}+o(\epsilon^{2})$
leading to

$\begin{split}\dot{\epsilon}p_{1}&=\mathcal{L}_{x0}\phi_{0}+(\mathcal{L}_{x0}p_{1}+\mathcal{L}_{x1}\phi_{0})\epsilon+o(\epsilon^{2})\\\
&\simeq(\mathcal{L}_{x0}p_{1}+\mathcal{L}_{x1}\phi_{0})\epsilon\end{split}\,.$

The Laplace transform of both hand side of this equation affects only
$\epsilon(t)$, leading to the disappearance of its transform $\epsilon(s)$:

$s\,p_{1}(x)=\mathcal{L}_{x0}p_{1}(x)+F_{0}(x)$ (34)

with $F_{0}(x)=\mathcal{L}_{x1}\phi_{0}$. As in [7], from this equation
$p_{1}(x)$ can be derived as the sum of two solutions: the homogeneous
equation $(\mathcal{L}_{x0}-s)p_{1}=0$, which is the same as our usual
spectral equation with $s=\lambda_{n}$, and a particular solution $Q(x,s)$ to
be determined:

$p_{1}(x)=\left\\{\begin{matrix}a\,f_{1}(x)+b\,f_{2}(x)+Q(x)&{x_{\mathrm{r}}}<x<{x_{\mathrm{t}}}\\\
d\,f_{2}(x)+Q(x)&x<{x_{\mathrm{r}}}\end{matrix}\right.\,.$

For the sake of clarity, the dependency on $s$ here is implicit. As in the
above derivation, the usual boundary conditions lead to a system of linear
equations in $a$, $b$ and $d$ to be solved. After some algebra and referring
to Eq. (32) (see Appendix A for details), we found that:

$p_{1}^{\prime}({x_{\mathrm{t}}})=\frac{\displaystyle\frac{\mathcal{W}_{Q,f_{2}}({x_{\mathrm{t}}})}{\mathcal{W}({x_{\mathrm{t}}})Z}-\frac{\mathcal{W}_{\Delta
Q,f_{2}}({x_{\mathrm{r}}})}{\mathcal{W}({x_{\mathrm{r}}})Z}}{\psi({x_{\mathrm{t}}})-\psi({x_{\mathrm{r}}})}$
(35)

where the Wronskian $\mathcal{W}(x)$ is from Eq. (31), $\Delta
Q(x)=\lim_{\delta\to 0}Q(x+|\delta|)-Q(x-|\delta|)$ and
$\mathcal{W}_{f,g}(x)=f^{\prime}(x)g(x)-f(x)g^{\prime}(x)$. We remark that
this formula is valid in general for any IF neuron model. Besides it depends
only on the particular solution $Q(x)$ and on the eigenfunction $\psi$, as
from Eq. (32) $f_{2}(x)=\psi(x)\,Z\,\mathcal{W}(x)$ with
$Z(s)=Z_{n}|_{\lambda_{n}=s}$. As above all these function implicitly depend
also on $s$ and then on the eigenvalues if $s=\lambda_{n}$.

###  III.4 Matching the transfer functions to derive $c_{n0}$

Focusing now on the specific case of the LIF neuron and considering an input
change affecting only the drift $\mu$, we have
$F_{0}(x)=-\phi_{0}^{\prime}(x)/\sigma_{0}$. The related particular solution
of Eq. (34) can be verified to be [7]

$Q_{\mu}(x,s)=-\frac{\phi_{0}^{\prime}(x)}{\sigma_{0}(s+1)}\,.$

Recalling that $\nu_{0}=\Phi=-\phi_{0}^{\prime}({x_{\mathrm{t}}})/2$, we
remark that $Q_{\mu}({x_{\mathrm{t}}},s)=2\nu_{0}/[\sigma_{0}(s+1)]$
highlighting the direct link between $Q(x)$ and the flux of realizations
$S_{\phi_{0}}(x)$ under stationary condition. Making use of all this
information in Eq. (35) (see Appendix B for details), we can derive
$p_{1}^{\prime}$ and hence from Eq. (28), the transfer function of the
modulated drift $\mu(t)$ of the input current

$H_{\mu}^{(P)}(s)=\frac{\nu_{0}}{\sigma_{0}(s+1)}\frac{\partial_{x}\psi({x_{\mathrm{t}}},s)-\partial_{x}\psi({x_{\mathrm{r}}},s)}{\psi({x_{\mathrm{t}}},s)-\psi_{(}{x_{\mathrm{r}}},s)}\,.$
(36)

With a similar derivation detailed in Appendix B, the transfer function of the
perturbations to the variance $\sigma^{2}(t)$ result to be

$H_{\sigma}^{(P)}(s)=\frac{\nu_{0}}{2\sigma_{0}^{2}(s+2)}\left(2-\frac{{x_{\mathrm{t}}}\partial_{x}\psi({x_{\mathrm{t}}},s)-{x_{\mathrm{r}}}\partial_{x}\psi({x_{\mathrm{r}}},s)}{\psi({x_{\mathrm{r}}},s)-\psi_{(}{x_{\mathrm{r}}},s)}\right)$
(37)

As under the mean-field approximation we know that both the mean and the
variance of the input current depend on the firing rate $\nu$, from these
expressions the rate-to-rate transfer function $H_{\nu}(s)$ can be eventually
carried out:

$H_{\nu}(s)=H_{\mu}(s)\frac{d\mu}{d\nu}+H_{\sigma}(s)\frac{d\sigma^{2}}{d\nu}$

We conclude deriving an expression for the coupling coefficients between
nonstationary and stationary modes
$c_{n0}^{(\epsilon)}=\\{\vec{c}_{\epsilon}\\}_{n}$ found in Eqs. (27) and
(28). To this purpose we note that the transfer function $H_{\epsilon}(s)$
(where $s=i\omega$) derived with the spectral expansion approach, can give

$\mathrm{Res}_{s=\lambda_{k}}H_{\epsilon}(s)=\lambda_{n}\,c_{n0}^{(\epsilon)}\,.$

Due to the equivalence between this $H_{\epsilon}(s)$ and
$H_{\epsilon}^{(P)}(s)$ from the above perturbative approach, we can carry out
an explicit expression for the coupling coefficients as

$c_{n0}^{(\epsilon)}=\frac{1}{\lambda_{n}}\mathrm{Res}_{s=\lambda_{k}}H_{\epsilon}^{(P)}(s)\,.$

To compute the residues from Eqs. (36) and (37), we point out that the
spectral equation Eq. (21) can be expressed directly in terms of adjoint
eigenfunction:
$\psi({v_{\mathrm{res}}},\lambda_{k})=\psi({v_{\mathrm{thr}}},\lambda_{k})$.
From this we have

$\mathrm{Res}_{s=\lambda_{k}}\frac{1}{\psi({v_{\mathrm{thr}}},\lambda_{k})-\psi({v_{\mathrm{res}}},\lambda_{k})}=1\,,$

which straightforwardly leads to

$c^{(\mu)}_{n0}=\frac{\nu_{0}}{\sigma_{0}\lambda_{n}(\lambda_{n}+1)}\left(\psi_{n}^{\prime}({x_{\mathrm{t}}})-\psi_{n}^{\prime}({x_{\mathrm{r}}})\right)$
(38)

and to

$c^{(\sigma)}_{n0}=-\frac{\nu_{0}}{2\sigma_{0}^{2}\lambda_{n}(\lambda_{n}+2)}\left({x_{\mathrm{t}}}\psi_{n}^{\prime}({x_{\mathrm{t}}})-{x_{\mathrm{r}}}\psi_{n}^{\prime}({x_{\mathrm{r}}})\right)\,.$
(39)

Remarkably, as for the LIF neuron we have

$\psi_{n}(x)=\frac{e^{x^{2}/2}}{Z_{n}}\mathcal{D}_{-\lambda_{n}}(-\sqrt{2}x)\,,$
(40)

once we known the eigenvalues $\lambda_{n}$, the expressions for both
$c^{(\mu)}_{n0}$ and $c^{(\sigma)}_{n0}$ are closed formulas, not requiring
any integral on the $x$-domain. Finally, recalling the dependence on $\nu$ of
the moments of synaptic current in Eq. (3), the coupling coefficients in the
firing rate equation Eq. (27) are

$\begin{split}c_{n0}&=\langle\partial_{\nu}\psi_{n}|\phi_{0}\rangle\\\
&=\frac{d\mu}{d\nu}c^{(\mu)}_{n0}+\frac{d\sigma^{2}}{d\nu}c^{(\sigma)}_{n0}\\\
&={\tau_{\mathrm{m}}}KJ\left(c^{(\mu)}_{n0}+J(1+\Delta^{2})c^{(\sigma)}_{n0}\right)\end{split}\,,$
(41)

where we reintroduced explicitly ${\tau_{\mathrm{m}}}$.

##  IV Conclusion

Here, by carrying out and pairing the dynamics of the firing rate $\nu(t)$
derived from two different approaches, we found novel analytic expressions for
the coefficients underlying the spectral expansion of the population density
$p(v,t)$. In the case of an uncoupled set of spiking neurons, we used the
renewal theory in combination with the spectral expansion approach by pairing
the relaxation dynamics of $\nu(t)$ – part of our ‘Rosetta stone’ – eventually
obtaining sums of series in closed form. More specifically, we uncovered a new
expression of the current-to-rate gain function $\Phi(\mu,\sigma)$ in Eq. (14)
as series of the adjoint eigenfunctions evaluated at the reset potential
$\psi_{n}({v_{\mathrm{res}}})$. We also found a tight relationship between the
moments of the ISI distribution and a suited combination of the
$\psi_{n}({v_{\mathrm{res}}})$ and the eigenvalues $\lambda_{n}$ of the
Fokker-Planck operator. We remark that all these findings are model-
independent highlighting the descriptive power of the spectral expansion and
further extending recent results reported in [35].

From the same side of our Rosetta stone, we found an alternative way to
compute both $\psi_{n}({v_{\mathrm{res}}})$ and $\lambda_{n}$. Indeed, relying
on the residue theorem, Eqs. (18) and (19) show how to obtain from the Laplace
transform $\rho(s)$ of the ISI density both these coefficients. Notably, these
equations provide a great advantage in practical terms to estimate numerically
the eigenvalues $\lambda_{n}$ even when their analytic expressions are not
available. This is due to the fact that usually eigenvalues results from the
numerical search of the roots of the spectral equation (21) [26, 35], or of
the zeros of numerically derived functions [34, 4]. Evaluating residues
instead requires computation of a limit or of a line integral around the
border of a domain centered around some rough guess for $\lambda_{n}$. The
expected computational advantages resulting from this residue-based approach
can pave the way to effective implementations of the spectral expansion
formalism in the numerical integration of the network dynamics of spiking
neurons.

For the specific case of LIF neurons, we also derived additional analytic
expressions for the coupling coefficients $c_{n0}$ involved in the evaluation
of both the out-of-equilibrium dynamics and the critical points of coupled
networks via the firing rate equation (27). Here we resorted to the linear
response to a small perturbation of the input as the other side of our Rosetta
stone. Pairing the perturbative response carried out from the spectral
expansion approach and directly from the Fokker-Planck equation, we finally
obtained Eqs. (38) and (39). These expressions for $c_{n0}$ result to involve
only the eigenfunctions at the reset potential $\psi_{n}({v_{\mathrm{res}}})$
and the eigenvalues $\lambda_{n}$. These coefficients in turn appear to be
tightly linked in Eq. (16) to the ISI density $\rho(t)$ of isolated cells. In
other words the coupling coefficients are the expression of a single-neuron
feature valid under stationary condition, rather than being a direct function
of the synaptic efficacy $J$. Although this may appear as a contradiction, we
remark that in Eq. (41) these single-neuron features in $c_{n0}^{(\epsilon)}$
mainly play a modulatory role of the synaptic coupling $KJ$, which is fully
taken into account in the mean-field expressions for the current moments
$\mu(\nu)$ and $\sigma^{2}(\nu)$.

In the framework of the spectral expansion of $p(v,t)$, approximated
expressions for $c_{n0}^{(\mu)}$ has been previously derived in [21], finding

$c_{n0}(\mu)\simeq\frac{\nu_{0}}{\lambda_{n}+1}C(\mu)$

for networks of LIF neurons working in drift-dominated regime
($\mu{\tau_{\mathrm{m}}}>{v_{\mathrm{thr}}}$) and with small synaptic noise
(${x_{\mathrm{t}}}\gg 1$). This is perfectly compatible with what we derived
here and reported in Eq. (38), provided that
$C(\mu)=\left(\psi_{n}^{\prime}({x_{\mathrm{t}}})-\psi_{n}^{\prime}({x_{\mathrm{r}}})\right)/(\sigma_{0}\lambda_{n})$.
Note that our result has a more general applicability being valid also under
subthreshold noise-dominated regime
($\mu{\tau_{\mathrm{m}}}<{v_{\mathrm{thr}}}$). Another expression for
$c_{n0}^{(\mu)}$ has been recently presented also in [35]. In this case the
spectral expansion approach targeted the refractory density, leading to obtain
for a network of renewal neurons with Poissonian ISIs and a refractory period
$\tau_{0}$, the following expression

$c_{n0}(\mu)=\frac{\partial_{\mu}\Phi}{\lambda_{n}+(1-\tau_{0}\lambda_{n})\nu_{0}}\xrightarrow[\tau_{0}\to
0]{}\frac{\partial_{\mu}\Phi}{\lambda_{n}+\nu_{0}}\,.$

This equation appears to be qualitatively different from our Eq. (38). A
possible explanation of such disagreement on one hand may reside in the fact
that LIF neurons are not Poissonian processes. On the other hand, this might
highlight a qualitative difference in the population dynamics described by the
refractory density approach compared to the classical one focused on the
evolution of the membrane potential density.

Note that, our derivation of the coupling terms $c_{n0}$ for LIF neurons can
in principle be applied to other IF neuron models. We expect in this case
qualitatively similar results, confirming that the ISI density-based
coefficients $c_{n0}^{(\epsilon)}$ eventually modulate the system response
with strength proportional to $KJ$. We also draw the reader attention to the
fact that here we did not report any expression for the coupling coefficients
$C_{nm}$ between nonstationary modes ($n,m\neq 0$). They are involved in the
nonlinear response of the neuronal network playing a role when the time
derivative $\dot{\nu}$ of the firing rate is relatively large and the density
$p(v,t)$ is significantly different from the stationary one ($\phi_{0}$) [26].
Provided that such constraints are not violated, the network dynamics can be
fully described by taking into account the only coefficients of the spectral
expansion we studied, even when the firing rate is driven outside equilibrium
[4, 28].

Regarding the applicability to the out-of-equilibrium dynamics, we further
point out that all the mentioned coefficients depend on the infinitesimal
moments $\mu$ and $\sigma$ of the synaptic current. This means that under the
extended mean-field approximation that leads to Eq. (3), all the coefficients
depend from time to time on the activity $\nu(t)$ of the network. Here $\mu$
and $\sigma$ must be interpreted as parameters. Indeed, it is true that our
Rosetta stone approach allowed us to carry out these spectral expansion
coefficients from the quasi-equilibrium dynamics (linear response and
relaxation in uncoupled sets), i.e., considering $\mu$ and $\sigma$ as fixed.
Nonetheless, for a given $\mu$ and $\sigma$ they contribute to faithful
represent a snapshot of the moving basis
$\left\\{|\phi_{n}\rangle\right\\}_{n\in\mathbb{Z}}$ onto which the density
$p(v,t)$ is projected. The basis moves as the drift and diffusion terms in the
Fokker-Planck equation (5) vary describing a sequence of stochastic processes
(i.e., the membrane potentials of the neurons) locally homogeneous in time. We
believe this is an important point to stress. Starting from the linearizable
(quasi-equilibrium) dynamics of spiking neuron networks, the spectral
expansion coefficients can be derived paving the way to describe the same
system outside equilibrium.

##  Appendix A Derivation of Eq. (35)

To derive the formula for the perturbation of the density
$p_{1}({x_{\mathrm{t}}})$ we need first to find the following coefficients
$a$, $b$ and $d$:

$p_{1}(x)=\left\\{\begin{matrix}a\,f_{1}(x)+b\,f_{2}(x)+Q(x)&{x_{\mathrm{r}}}<x<{x_{\mathrm{t}}}\\\
d\,f_{2}(x)+Q(x)&x<{x_{\mathrm{r}}}\end{matrix}\right.\,.$

They can be evaluated by taking into account the boundary conditions of the
Fokker-Planck equation:

$\begin{split}af_{1}({x_{\mathrm{t}}})+bf_{2}({x_{\mathrm{t}}})&=-Q({x_{\mathrm{t}}})\\\
a\left[f^{\prime}_{1}({x_{\mathrm{t}}})-f^{\prime}_{1}({x_{\mathrm{r}}})\right]+&\\\
b\left[f^{\prime}_{2}({x_{\mathrm{t}}})-f^{\prime}_{2}({x_{\mathrm{r}}})\right]+d\,f^{\prime}_{2}({x_{\mathrm{r}}})&=\Delta
Q^{\prime}({x_{\mathrm{r}}})-Q^{\prime}({x_{\mathrm{t}}})\\\
a\,f_{1}({x_{\mathrm{r}}})+b\,f_{2}({x_{\mathrm{r}}})-d\,f_{2}({x_{\mathrm{r}}})&=-\Delta
Q({x_{\mathrm{r}}})\end{split}\,.$

We note that the homogeneous system is equivalent to the eigenfunction
problem. Indeed, the setting of the associated determinant to zero gives back
the spectral equation, i.e. the boundary condition
$\psi_{n}({x_{\mathrm{r}}})=\psi_{n}({x_{\mathrm{r}}})$, for the eigenvalues
$\lambda_{n}$ [14]. Solving the above system we have

$a=\frac{Q({x_{\mathrm{t}}})f_{2}({x_{\mathrm{r}}})f^{\prime}_{2}({x_{\mathrm{t}}})-f_{2}({x_{\mathrm{t}}})\left[\Delta
Q({x_{\mathrm{r}}})f^{\prime}_{2}({x_{\mathrm{r}}})+f_{2}({x_{\mathrm{r}}})\left(Q^{\prime}({x_{\mathrm{t}}})-\Delta
Q^{\prime}({x_{\mathrm{t}}})\right)\right]}{\mathcal{W}({x_{\mathrm{t}}})f_{2}({x_{\mathrm{r}}})-\mathcal{W}({x_{\mathrm{r}}})f_{2}({x_{\mathrm{t}}})}$

and

$b=-\frac{f_{2}({x_{\mathrm{r}}})\left[Q_{(}{x_{\mathrm{r}}})\left(f^{\prime}_{1}({x_{\mathrm{r}}})-f^{\prime}_{1}({x_{\mathrm{t}}})\right)+f_{1}({x_{\mathrm{t}}})\left(Q^{\prime}({x_{\mathrm{t}}})-\Delta
Q^{\prime}({x_{\mathrm{r}}})\right)\right]+f^{\prime}_{2}({x_{\mathrm{t}}})\left[\Delta
Q({x_{\mathrm{r}}})f_{1}({x_{\mathrm{t}}})-Q({x_{\mathrm{t}}})f_{1}({x_{\mathrm{r}}})\right]}{\mathcal{W}({x_{\mathrm{t}}})f_{2}({x_{\mathrm{r}}})-\mathcal{W}({x_{\mathrm{r}}})f_{2}({x_{\mathrm{t}}})}\,.$

These coefficients determine $p_{1}(x)$ in the domain of interest
($x>{x_{\mathrm{r}}}$), this is why we do not report the expression for $d$.
Evaluating the derivative $p_{1}^{\prime}({x_{\mathrm{t}}})$ of the above
density perturbation at the threshold potential eventually leads to Eq. (35).

##  Appendix B Derivation of the transfer function $H^{(P)}_{\epsilon}$

From Eq. (33) we know that the generic transfer function $H^{(P)}_{\epsilon}$
is proportional to $p_{1}^{\prime}({x_{\mathrm{t}}})$ given by Eq. (35), whose
derivation has been detailed in the previous Appendix. Starting from this,
here we provide the key steps needed to obtain the explicit expressions for
$H^{(P)}_{\epsilon}$ reported in Eqs. (36) and (37) for the specific case of
the standard LIF neuron model. If the input modulation involves only the mean
$\mu$, $Q_{\mu}(x)\propto\phi^{\prime}_{0}(x)$. As in the main text, here we
assume an implicit dependence on $s$. Remembering that
$\mathcal{L}_{x0}\phi_{0}=\partial_{x}(x\phi_{0})+\frac{1}{2}\partial_{x}^{2}\phi_{0}=0$
we can write

$Q^{\prime}_{\mu}(x)=-2x\,Q_{\mu}(x)-2\phi_{0}(x)=\frac{\mathcal{W}^{\prime}(x)}{\mathcal{W}(x)}\,Q_{\mu}(x)-2\phi_{0}(x)\,,$

where we taken into account that from Eq. (31) the Wronskian for the LIF
neuron reduces to the differential equation
$\mathcal{W}^{\prime}(x)=-2x\,\mathcal{W}(x)$. The first term of the numerator
in Eq. (35) then reduces to

$\begin{split}\frac{\mathcal{W}_{Q,f_{2}}({x_{\mathrm{t}}})}{\mathcal{W}({x_{\mathrm{t}}})Z}&=\left.\frac{Q^{\prime}_{\mu}f_{2}-Q_{\mu}f_{2}}{\mathcal{W}Z}\right|_{x={x_{\mathrm{t}}}}\\\
&=\left.-\frac{Q_{\mu}}{Z}\left(-\frac{\mathcal{W}^{\prime}f_{2}}{\mathcal{W}^{2}}+\frac{f^{\prime}_{2}}{\mathcal{W}}\right)-2\phi_{0}\frac{f_{2}}{\mathcal{W}Z}\right|_{x={x_{\mathrm{t}}}}\\\
&=\left.-\left(Q_{\mu}+2\phi_{0}\right)\psi^{\prime}\right|_{x={x_{\mathrm{t}}}}\\\
&=-Q_{\mu}({x_{\mathrm{t}}})\psi^{\prime}({x_{\mathrm{t}}})\end{split}$

where we used Eq. (32) replacing $\lambda_{n}$ with $s$ and defining
$Z(s)=Z_{n}|_{\lambda_{n}=s}$, and we took into account the boundary condition
$\phi_{0}({x_{\mathrm{t}}})=0$. The same can be done with the other term
$\mathcal{W}_{\Delta
Q,f_{2}}({x_{\mathrm{r}}})/\mathcal{W}({x_{\mathrm{r}}})Z$ by taking into
account the other boundary condition about the flux conservation,
$Q(x_{t})=\Delta Q(x_{r})$, eventually leading to rewrite Eq. (35) as

$p^{\prime}_{1}({x_{\mathrm{t}}})=-Q_{\mu}(x_{t})\frac{\psi^{\prime}(x_{t})-\psi^{\prime}(x_{r})}{\psi(x_{t})-\psi(x_{r})}\,.$

Analogous results hold whenever the particular solution is proportional to the
derivative of the stationary eigenfunction.

In the case of a changing variance $\sigma^{2}$ of the input current, the
forcing term is
$F_{0}(t)=\frac{\nu_{1}(t)}{\sigma_{0}^{2}}\phi^{\prime\prime}_{0}$. The
particular solution is then

$Q^{\sigma}(x)=K(s)\phi^{\prime\prime}_{0}(x)=\frac{\gamma(s)}{2(s+2)}\phi^{\prime\prime}_{0}(x)\,.$

Using the stationary equation $\mathcal{L}_{x0}\phi_{0}=0$, we have that

$\partial^{3}_{x}\phi_{0}=\left(-2x\partial^{2}_{x}-4\partial_{x}\right)\phi_{0}=(4x^{2}-4)\phi^{\prime}_{0}\,.$

By making use of these expressions in Eq. (35), we eventually obtain for the
first term of the numerator

$\frac{\mathcal{W}_{Q,f_{2}}(x)}{\mathcal{W}(x)Z}=\frac{K\phi^{\prime}_{0}}{\mathcal{W}Z}\left(-4f_{2}+2x(2xf_{2}+f^{\prime}_{2})\right)\,.$

As above, the differential definition of $\mathcal{W}$ and Eq. (32) lead us to
derive the results presented in the main text. Note that, in order to obtain
the transfer function $H^{(P)}_{\epsilon}$ the complete definition of the
variation of the firing rate is $\nu_{1}=-\frac{1}{2}p_{1}^{\prime}(x)$.

## References

  * [1] Larry F Abbott and Carl van Vreeswijk.  Asynchronous states in networks of pulse-coupled oscillators.  Phys. Rev. E, 48(2):1483–1490, 1993.
  * [2] Daniel J Amit and Nicolas Brunel.  Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex.  Cereb. Cortex, 7(3):237–52, 1997.
  * [3] Daniel J Amit and Misha Tsodyks.  Quantitative study of attractor neural network retrieving at low spike rates: I. substrate—spikes, rates and neuronal gain.  Network, 2(3):259–273, 1991.
  * [4] Moritz Augustin, Josef Ladenbauer, Fabian Baumann, and Klaus Obermayer.  Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: Comparison and implementation.  PLoS Comput. Biol., 13(6):e1005545, 2017.
  * [5] Nicolas Brunel.  Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons.  J. Comput. Neurosci., 8(3):183–208, 2000.
  * [6] Nicolas Brunel, Frances S Chance, Nicolas Fourcaud, and Larry F Abbott.  Effects of synaptic noise and filtering on the frequency response of spiking neurons.  Phys. Rev. Lett., 86(10):2186–9, 2001.
  * [7] Nicolas Brunel and Vincent Hakim.  Fast global oscillations in networks of integrate-and-fire neurons with low firing rates.  Neural Comput., 11(7):1621–1671, 1999.
  * [8] David Cai, Louis Tao, Michael Shelley, and David W McLaughlin.  An effective kinetic representation of fluctuation-driven neuronal networks with application to simple and complex cells in visual cortex.  Proc. Natl. Acad. Sci. USA, 101(20):7757–62, 2004.
  * [9] R M Capocelli and Luigi M Ricciardi.  Diffusion approximation and first passage time problem for a model neuron.  Kybernetik, 8(6):214–23, 1971.
  * [10] D R Cox and H D Miller.  The theory of stochastic processes.  Chapman and Hall/CRC, 1977.
  * [11] D A Darling and Arnold Siegert.  The first passage problem for a continuous Markov process.  Ann. Math. Stat., 24(4):624–639, 1953.
  * [12] Marc de Kamps.  A simple and stable numerical solution for the population density equation.  Neural Comput., 15(9):2129–46, 2003.
  * [13] Javier DeFelipe, Lidia Alonso-Nanclares, and Jon I Arellano.  Microstructure of the neocortex: comparative aspects.  J. Neurocytol., 31(3-5):299–316, 2002.
  * [14] Taşkın Deniz and Stefan Rotter.  Solving the two-dimensional Fokker-Planck equation for strongly correlated neurons.  Phys. Rev. E, 95(1-1):012412, 2017.
  * [15] Ronald F Fox.  Stochastic resonance in a double well.  Phys. Rev. A, 39(8):4148–4153, 1989.
  * [16] Stefano Fusi and Maurizio Mattia.  Collective behavior of networks with linear (VLSI) integrate-and-fire neurons.  Neural Comput., 11(3):633–52, 1999.
  * [17] Luca Gammaitoni, Peter Hänggi, Peter Jung, and Fabio Marchesoni.  Stochastic resonance.  Rev. Mod. Phys., 70(1):223–287, 1998.
  * [18] George L Gerstein and Benoit Mandelbrot.  Random walk models for the spike activity of a single neuron.  Biophys. J., 4(1):41–68, 1964.
  * [19] Wulfram Gerstner, Werner M. Kistler, Richard Naud, and Liam Paninski.  Neuronal dynamics: From single neurons to networks and models of cognition.  Cambridge University Press, 2014.
  * [20] Bruce W Knight.  Dynamics of encoding in neuron populations: some general mathematical features.  Neural Comput., 12(3):473–518, 2000.
  * [21] Bruce W Knight, Dimitri Manin, and Lawrence Sirovich.  Dynamical models of interacting neuron populations in visual cortex.  In E.C. Gerf, editor, Symposium on Robotics and Cybernetics: Computational Engineering in Systems Applications, pages 1–5, Cite Scientifique, Lille, France, 1996. Cite Scientifique.
  * [22] L D Landau and E M Lifshits.  Statistical physics.  Pergamon Press, Oxford, 2nd edition, 1969.
  * [23] Benjamin Lindner and Lutz Schimansky-Geier.  Transmission of noise coded versus additive signals through a neuronal ensemble.  Phys. Rev. Lett., 86(14):2934–2937, 2001.
  * [24] Henry Markram, Eilif Muller, Srikanth Ramaswamy, Michael W Reimann, Marwan Abdellah, and et al.  Reconstruction and simulation of neocortical microcircuitry.  Cell, 163(2):456–492, 2015.
  * [25] Maurizio Mattia, Matteo Biggio, Andrea Galluzzi, and Marco Storace.  Dimensional reduction in networks of non-Markovian spiking neurons: Equivalence of synaptic filtering and heterogeneous propagation delays.  PLoS Comput. Biol., 15(10):e1007404, 2019.
  * [26] Maurizio Mattia and Paolo Del Giudice.  Population dynamics of interacting spiking neurons.  Phys. Rev. E, 66(5 Pt 1):051917, 2002.
  * [27] Maurizio Mattia and Paolo Del Giudice.  Finite-size dynamics of inhibitory and excitatory interacting spiking neurons.  Phys. Rev. E, 70:052903, 2004.
  * [28] Maurizio Mattia and Gianni V Vinci.  Low dimensional dynamics of spiking neuron networks.  Zenodo, 2021.  DOI: 10.5281/zenodo.5518215.
  * [29] A D McLachlan and M A Ball.  Time-dependent Hartree-Fock theory for molecules.  Rev. Mod. Phys., 36(3):844–855, 1964.
  * [30] Carver Mead.  Analog VLSI and neural system.  Addison-Wesley, Reading, MA, 1989.
  * [31] Duane Q Nykamp and Daniel Tranchina.  A population density approach that facilitates large-scale modeling of neural networks: analysis and an application to orientation tuning.  J. Comput. Neurosci., 8(1):19–50, 2000.
  * [32] Frank W J Olver, Daniel W Lozier, Ronald F Boisvert, and Charles W Clark.  NIST Handbook of Mathematical Functions.  Cambridge University Press, New York, NY, 2010.
  * [33] A Omurtag, Bruce W Knight, and Lawrence Sirovich.  On the simulation of large populations of neurons.  J. Comput. Neurosci, 8(1):51–63, 2000.
  * [34] Srdjan Ostojic.  Interspike interval distributions of spiking neurons driven by fluctuating inputs.  J. Neurophysiol., 106(1):361–73, 2011.
  * [35] Bastian Pietras, Noé Gallice, and Tilo Schwalger.  Low-dimensional firing-rate dynamics for populations of renewal-type spiking neurons.  Phys. Rev. E, 102(2):34–38, 2020.
  * [36] Christian Pozzorini, Skander Mensi, Olivier Hagens, Richard Naud, Christof Koch, and Wulfram Gerstner.  Automated high-throughput characterization of single neurons by means of simplified spiking models.  PLoS Comput. Biol., 11(6):1–29, 2015.
  * [37] Magnus J E Richardson.  Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive.  Phys Rev E, 76(2 Pt 1):021919, 2007.
  * [38] Magnus J E Richardson and Rupert Swarbrick.  Firing-rate response of a neuron receiving excitatory and inhibitory synaptic shot noise.  Phys. Rev. Lett., 105(17):1–4, 2010.
  * [39] Hannes Risken.  The Fokker-Planck Equation.  Springer, Berlin, Heidelberg, 1984.
  * [40] Evan S Schaffer, Srdjan Ostojic, and Larry F Abbott.  A Complex-Valued Firing-Rate Model That Approximates the Dynamics of Spiking Networks.  PLoS Comput. Biol., 9(10):e1003301, 2013.
  * [41] Arnold Siegert.  On the first passage time probability problem.  Phys. Rev., 81(4):617–623, 1951.
  * [42] Corinne Teeter, Ramakrishnan Iyer, Vilas Menon, Nathan Gouwens, David Feng, Jim Berg, Aaron Szafer, Nicholas Cain, Hongkui Zeng, Michael Hawrylycz, Christof Koch, and Stefan Mihalas.  Generalized leaky integrate-and-fire models classify multiple neuron types.  Nat. Commun., 9(1):709, 2018.
  * [43] Alessandro Treves.  Mean-field analysis of neuronal spike dynamics.  Network, 4(3):259–84, 1993.
  * [44] Henry C Tuckwell.  Introduction to theoretical neurobiology: volume 2, nonlinear and stochastic theories, volume 8.  Cambridge University Press, 1988.
  * [45] Gianni V Vinci, Lorenzo Fasano, and Maurizio Mattia.  Universal population dynamics with simplified spiking neurons.  Zenodo, 2021.  DOI: 10.5281/zenodo.5646222.