Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS Comput Biol**|**v.13(6); 2017 June**|**PMC5507472

Formats

Article sections

Authors

Related links

PLoS Comput Biol. 2017 June; 13(6): e1005545.

Published online 2017 June 23. doi: 10.1371/journal.pcbi.1005545

PMCID: PMC5507472

Moritz Augustin,^{#}^{1,}^{2,}^{*} Josef Ladenbauer,^{#}^{1,}^{2,}^{3,}^{*} Fabian Baumann,^{1,}^{2} and Klaus Obermayer^{1,}^{2}

Ralf Haefner, Editor^{}

University of Rochester, UNITED STATES,

The authors have declared that no competing interests exist.

**Conceptualization:**MA JL.**Data curation:**MA JL FB.**Formal analysis:**MA JL.**Funding acquisition:**MA JL KO.**Investigation:**MA JL FB.**Methodology:**MA JL FB.**Project administration:**MA JL.**Resources:**KO.**Software:**MA JL FB.**Supervision:**MA JL.**Validation:**MA JL FB KO.**Visualization:**MA JL FB.**Writing – original draft:**MA JL.**Writing – review & editing:**MA JL FB KO.

Received 2016 September 30; Accepted 2017 May 1.

Copyright © 2017 Augustin et al

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

The spiking activity of single neurons can be well described by a nonlinear integrate-and-fire model that includes somatic adaptation. When exposed to fluctuating inputs sparsely coupled populations of these model neurons exhibit stochastic collective dynamics that can be effectively characterized using the Fokker-Planck equation. This approach, however, leads to a model with an infinite-dimensional state space and non-standard boundary conditions. Here we derive from that description four simple models for the spike rate dynamics in terms of low-dimensional ordinary differential equations using two different reduction techniques: one uses the spectral decomposition of the Fokker-Planck operator, the other is based on a cascade of two linear filters and a nonlinearity, which are determined from the Fokker-Planck equation and semi-analytically approximated. We evaluate the reduced models for a wide range of biologically plausible input statistics and find that both approximation approaches lead to spike rate models that accurately reproduce the spiking behavior of the underlying adaptive integrate-and-fire population. Particularly the cascade-based models are overall most accurate and robust, especially in the sensitive region of rapidly changing input. For the mean-driven regime, when input fluctuations are not too strong and fast, however, the best performing model is based on the spectral decomposition. The low-dimensional models also well reproduce stable oscillatory spike rate dynamics that are generated either by recurrent synaptic excitation and neuronal adaptation or through delayed inhibitory synaptic feedback. The computational demands of the reduced models are very low but the implementation complexity differs between the different model variants. Therefore we have made available implementations that allow to numerically integrate the low-dimensional spike rate models as well as the Fokker-Planck partial differential equation in efficient ways for arbitrary model parametrizations as open source software. The derived spike rate descriptions retain a direct link to the properties of single neurons, allow for convenient mathematical analyses of network states, and are well suited for application in neural mass/mean-field based brain network models.

Characterizing the dynamics of biophysically modeled, large neuronal networks usually involves extensive numerical simulations. As an alternative to this expensive procedure we propose efficient models that describe the network activity in terms of a few ordinary differential equations. These systems are simple to solve and allow for convenient investigations of asynchronous, oscillatory or chaotic network states because linear stability analyses and powerful related methods are readily applicable. We build upon two research lines on which substantial efforts have been exerted in the last two decades: (i) the development of single neuron models of reduced complexity that can accurately reproduce a large repertoire of observed neuronal behavior, and (ii) different approaches to approximate the Fokker-Planck equation that represents the collective dynamics of large neuronal networks. We combine these advances and extend recent approximation methods of the latter kind to obtain spike rate models that surprisingly well reproduce the macroscopic dynamics of the underlying neuronal network. At the same time the microscopic properties are retained through the single neuron model parameters. To enable a fast adoption we have released an efficient Python implementation as open source software under a free license.

There is prominent evidence that information in the brain, about a particular stimulus for example, is contained in the collective neuronal spiking activity averaged over populations of neurons with similar properties (population spike rate code) [1, 2]. Although these populations can comprise a large number of neurons [3], they often exhibit low-dimensional collective spiking dynamics [4] that can be measured using neural mass signals such as the local field potential or electroencephalography.

The behavior of cortical networks at that level is often studied computationally by employing simulations of multiple (realistically large or subsampled) populations of synaptically coupled individual spiking model neurons. A popular choice of single cell description for this purpose are two-variable integrate-and-fire models [5, 6] which describe the evolution of the fast (somatic) membrane voltage and an adaptation variable that represents a slowly-decaying potassium current. These models are computationally efficient and can be successfully calibrated using electrophysiological recordings of real cortical neurons and standard stimulation protocols [5, 7–10] to accurately reproduce their subthreshold and spiking activity. The choice of such (simple) neuron models, however, does not imply reasonable (short enough) simulation durations for a recurrent network, especially when large numbers of neurons and synaptic connections between them are considered.

A fast and mathematically tractable alternative to simulations of large networks are population activity models in terms of low-dimensional ordinary differential equations (i.e., which consist of only a few variables) that typically describe the evolution of the spike rate. These reduced models can be rapidly solved and allow for convenient analyses of the dynamical network states using well-known methods that are simple to implement. A popular example are the Wilson-Cowan equations [11], which were also extended to account for (slow) neuronal adaptation [12] and short-term synaptic depression [13]. Models of this type have been successfully applied to qualitatively characterize the possible dynamical states of coupled neuronal populations using phase space analyses [11–13], yet a direct link to more biophysically described networks of (calibrated) spiking neurons in terms of model parameters is missing.

Recently, derived population activity models have been proposed that bridge the gap between single neuron properties and mesoscopic network dynamics. These models are described by integral equations [14, 15] or partial differential equations [16, 17]

Here we derive simple models in terms of low-dimensional ordinary differential equations (ODEs) for the spike rate dynamics of sparsely coupled adaptive nonlinear integrate-and-fire neurons that are exposed to noisy synaptic input. The derivations are based on a Fokker-Planck equation that describes the neuronal population activity in the mean-field limit of large networks. We develop reduced models using recent methodological advances on two different approaches: the first is based on a spectral decomposition of the Fokker-Planck operator under two different slowness assumptions [18–20]. In the second approach we consider a cascade of linear temporal filters and a nonlinear function which are determined from the Fokker-Planck equation and semi-analytically approximated, building upon [21]. Both approaches are extended for an adaptation current, a nonlinear spike generating current and recurrent coupling with distributed synaptic delays.

We evaluate the developed low-dimensional spike rate models quantitatively in terms of reproduction accuracy in a systematic manner over a wide range of biologically plausible parameter values. In addition, we provide numerical implementations for the different reduction methods as well as the Fokker-Planck equation under a free license as open source project.

For the derived models in this contribution we use the adaptive exponential integrate-and-fire (aEIF) model [5] to describe individual neurons, which is similar to the model proposed by Izhikevich [6] but includes biophysically meaningful parameters and a refined description of spike initiation. However, the presented derivations are equally applicable when using the Izhikevich model instead (requiring only a small number of simple substitutions in the code).

Through their parameters the derived models retain a direct, quantitative link to the underlying spiking model neurons, and they are described in a well-established, convenient form (ODEs) that can be rapidly solved and analyzed. Therefore, these models are well suited (i) for mathematical analyses of dynamical states at the population level, e.g., linear stability analyses of attractors, and (ii) for application in multi-population brain network models. Apart from a specific network setting, the derived models are also appropriate as a spike rate description of individual neurons under noisy input conditions.

The structure of this article contains mildly redundant model specifications allowing the readers who are not interested in the methodological foundation to directly read the self-contained Sect. Results.

The quantity of our interest is the population-averaged number of spikes emitted by a large homogeneous network of *N* sparsely coupled aEIF model neurons per small time interval, i.e., the spike rate *r*_{N}(*t*). The state of neuron *i* at time *t* is described by the membrane voltage *V*_{i}(*t*) and adaptation current *w*_{i}(*t*), which evolve piecewise continuously in response to overall synaptic current *I*_{syn,i} = *I*_{ext,i}(*t*) + *I*_{rec,i}(*t*). This input current consists of fluctuating network-external drive *I*_{ext,i} = *C*[*μ*_{ext}(*t*) + *σ*_{ext}(*t*)*ξ*_{ext,i}(*t*)] with membrane capacitance *C*, time-varying moments *μ*_{ext}, ${\sigma}_{\text{ext}}^{2}$ and unit Gaussian white noise process *ξ*_{ext,i} as well as recurrent input *I*_{rec,i}. The latter causes delayed postsynaptic potentials (i.e., deflections of *V*_{i}) of small amplitude *J* triggered by the spikes of *K* presynaptic neurons (see Sect. Methods for details).

Here we present two approaches of how the spike rate dynamics of the large, stochastic delay-differential equation system for the 2*N* states (*V*_{i}, *w*_{i}) can be described by simple models in terms of low-dimensional ODEs. Both approaches (i) take into account adaptation current dynamics that are sufficiently slow, allowing to replace the individual adaptation current *w*_{i} by its population-average *w*, governed by

$$\frac{d\u27e8w\u27e9}{dt}=\frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{{\tau}_{w}}+b\phantom{\rule{0.166667em}{0ex}}r\left(t\right),$$

(1)

where *a*, *E*_{w}, *b*, *τ*_{w} are the adaptation current model parameters (subthreshold conductance, reversal potential, spike-triggered increment, time constant, respectively), *V*_{∞} is the steady-state membrane voltage averaged across the population (which can vary over time, see below), and *r* is the spike rate of the respective low-dimensional model. Furthermore, both approaches (ii) are based on the observation that the collective dynamics of a large, sparsely coupled (and noise driven) network of integrate-and-fire type neurons can be well described by a Fokker-Planck equation. In this *intermediate* Fokker-Planck (FP) model the overall synaptic input is approximated by a mean part with additive white Gaussian fluctuations, *I*_{syn,i}/*C* ≈ *μ*_{syn}(*t*, *r*_{d}) + *σ*_{syn}(*t*, *r*_{d})*ξ*_{i}(*t*), that are uncorrelated between neurons. The moments of the overall synaptic input,

$${\mu}_{\text{syn}}={\mu}_{\text{ext}}\left(t\right)+JK{r}_{d}\left(t\right),\phantom{\rule{2.em}{0ex}}{\sigma}_{\text{syn}}^{2}={\sigma}_{\text{ext}}^{2}\left(t\right)+{J}^{2}K{r}_{d}\left(t\right),$$

(2)

depend on time via the moments of the external input and, due to recurrent coupling, on the delayed spike rate *r*_{d}. The latter is governed by

$$\frac{d{r}_{d}}{dt}=\frac{r-{r}_{d}}{{\tau}_{d}},$$

(3)

which corresponds to individual propagation delays drawn from an exponentially distributed random variable with mean *τ*_{d}. The FP model involves solving a partial differential equation (PDE) to obtain the time-varying membrane voltage distribution *p*(*V*, *t*) and the spike rate *r*(*t*).

The first reduction approach is based on the spectral decomposition of the Fokker-Planck operator ℒ and leads to the following two low-dimensional models: the “basic” model variant (spec_{1}) is given by a complex-valued differential equation describing the spike rate evolution in its real part,

$$\frac{d\tilde{r}}{dt}={\lambda}_{1}(\tilde{r}-{r}_{\infty}),\phantom{\rule{2.em}{0ex}}r\left(t\right)=\text{Re}\left\{\tilde{r}\right\},$$

(4)

where λ_{1}(*μ*_{tot}, *σ*_{tot}) is the dominant eigenvalue of ℒ and *r*_{∞}(*μ*_{tot}, *σ*_{tot}) is the steady-state spike rate. Its parameters λ_{1}, *r*_{∞}, and *V*_{∞} (cf. Eq (1)) depend on the *total* input moments given by *μ*_{tot}(*t*) = *μ*_{syn} − *w*/*C* and ${\sigma}_{\text{tot}}^{2}\left(t\right)={\sigma}_{\text{syn}}^{2}$ which closes the model (Eqs (1)–(4)). The other, “advanced” spectral model variant (spec_{2}) is given by a real-valued second order differential equation for the spike rate,

$${\beta}_{2}\phantom{\rule{0.166667em}{0ex}}\ddot{r}+{\beta}_{1}\phantom{\rule{0.166667em}{0ex}}\dot{r}+{\beta}_{0}\phantom{\rule{0.166667em}{0ex}}r={r}_{\infty}-r-{\beta}_{c},$$

(5)

where the dots denote time derivatives. Its parameters *β*_{2}, *β*_{1}, *β*_{0}, *β*_{c}, *r*_{∞} and *V*_{∞} depend on the total input moments (*μ*_{tot}, ${\sigma}_{\text{tot}}^{2}$) as follows: the latter two parameters explicitly as in the basic model above, the former four indirectly via the first two dominant eigenvalues λ_{1}, λ_{2} and via additional quantities obtained from the (stationary and the first two nonstationary) eigenfunctions of ℒ and its adjoint ℒ^{*}. Furthermore, the parameter *β*_{c} depends explicitly on the population-averaged adaptation current *w*, the delayed spike rate *r*_{d}, and on the first and second order time derivatives of the external moments *μ*_{ext} and ${\sigma}_{\text{ext}}^{2}$.

The second approach is based on a Linear-Nonlinear (LN) cascade, in which the population spike rate is generated by applying to the time-varying mean and standard deviation of the overall synaptic input, *μ*_{syn} and *σ*_{syn}, separately a linear temporal filter, followed by a common nonlinear function. These three components–two linear filters and a nonlinearity–are extracted from the Fokker-Planck equation. Approximating the linear filters using exponentials and damped oscillating functions yields two model variants: In the basic “exponential” (LN_{exp}) model the filtered mean *μ*_{f} and standard deviation *σ*_{f} of the overall synaptic input are given by

$$\frac{d{\mu}_{\mathrm{f}}}{dt}=\frac{{\mu}_{\text{syn}}-{\mu}_{\mathrm{f}}}{{\tau}_{\mu}},\phantom{\rule{2.em}{0ex}}\frac{d{\sigma}_{\mathrm{f}}}{dt}=\frac{{\sigma}_{\text{syn}}-{\sigma}_{\mathrm{f}}}{{\tau}_{\sigma}},$$

(6)

where the time constants *τ*_{μ}(*μ*_{eff}, *σ*_{eff}), *τ*_{σ}(*μ*_{eff}, *σ*_{eff}) depend on the *effective* (filtered) input mean *μ*_{eff}(*t*) = *μ*_{f} − *w*/*C* and standard deviation *σ*_{eff}(*t*) = *σ*_{f}. The “damped oscillator” (LN_{dos}) model variant, on the other hand, describes the filtered input moments by

$${\ddot{\mu}}_{\mathrm{f}}+\frac{2}{\tau}{\dot{\mu}}_{\mathrm{f}}+(\frac{2}{{\tau}^{2}}+{\omega}^{2})\phantom{\rule{1pt}{0ex}}{\mu}_{\mathrm{f}}=\frac{1+{\tau}^{2}{\omega}^{2}}{\tau}\phantom{\rule{1pt}{0ex}}(\frac{{\mu}_{\text{syn}}}{\tau}+{\dot{\mu}}_{\text{syn}}),$$

(7)

$$\frac{d{\sigma}_{\mathrm{f}}}{dt}=\frac{{\sigma}_{\text{syn}}-{\sigma}_{\mathrm{f}}}{{\tau}_{\sigma}},$$

(8)

where the time constants *τ*(*μ*_{tot}, *σ*_{tot}), *τ*_{σ}(*μ*_{tot}, *σ*_{tot}) and the angular frequency *ω*(*μ*_{tot}, *σ*_{tot}) depend on the *total* input moments defined above. In both LN model variants the spike rate is obtained by the nonlinear transformation of the *effective* input moments through the steady-state spike rate,

(9)

and the steady-state mean membrane voltage *V*_{∞} (cf. Eq (1)) is also evaluated at (*μ*_{eff}, *σ*_{eff}).

These four models (spec_{1}, spec_{2}, LN_{exp}, LN_{dos}) from both reduction approaches involve a number of parameters that depend on the strengths of synaptic input and adaptation current only via the *total* or *effective* input moments. We refer to these parameters as *quantities* below to distinguish them from fixed (independent) parameters. The computational complexity when numerically solving the models forward in time (for different parametrizations) can be greatly reduced by precomputing those quantities for a range of values for the *total*/*effective* input moments and using look-up tables during time integration. Changing any parameter value of the external input, the recurrent coupling or the adaptation current does not require renewed precomputations, enabling rapid explorations of parameter space and efficient (linear) stability analyses of network states.

The full specification of the “ground truth” system (network of aEIF neurons), the derivations of the intermediate description (FP model) and the low-dimensional spike rate models complemented by concrete numerical implementations are provided in Sect. Methods (that is complemented by the supporting material S1 Text). In Fig 1 we visualize the outputs of the different models using an example excitatory aEIF network exposed to external input with varying mean *μ*_{ext}(*t*) and standard deviation *σ*_{ext}(*t*).

Here, and in the subsequent two sections, we assess the accuracy of the four low-dimensional models to reproduce the spike rate dynamics of the underlying aEIF population. The *intermediate* FP model is included for reference. The derived models generate population activity in response to overall synaptic input moments *μ*_{syn} and ${\sigma}_{\text{syn}}^{2}$. These depend on time via the external moments *μ*_{ext}(*t*) and ${\sigma}_{\text{ext}}^{2}\left(t\right)$, and the delayed spike rate *r*_{d}(*t*). Therefore, it is instrumental to first consider an uncoupled population and suitable variations of external input moments that effectively mimic a range of biologically plausible presynaptic spike rate dynamics. This allows us to systematically compare the reproduction performance of the different models over a manageable parameter space (without *K*, *J*, *τ*_{d}), yet it provides useful information on the accuracy for recurrent networks.

For many network settings the dominant effect of synaptic coupling is on the mean input (cf. Eq (2)). Therefore, we consider first in detail time-varying mean but constant variance of the input. Specifically, to account for a wide range of oscillation frequencies for presynaptic spike rates, *μ*_{ext} is described by an Ornstein-Uhlenbeck (OU) process

$${\dot{\mu}}_{\text{ext}}=\frac{\overline{\mu}-{\mu}_{\text{ext}}}{{\tau}_{\text{ou}}^{\mu}}+\sqrt{\frac{2}{{\tau}_{\text{ou}}^{\mu}}}{\vartheta}_{\mu}\xi \left(t\right),$$

(10)

where ${\tau}_{\text{ou}}^{\mu}$ denotes the correlation time, $\overline{\mu}$ and _{μ} are the mean and standard deviation of the stationary normal distribution, i.e., ${lim}_{t\to \infty}{\mu}_{\text{ext}}\left(t\right)\sim \mathcal{N}(\overline{\mu},{\vartheta}_{\mu}^{2})$, and *ξ* is a unit Gaussian white noise process. Sample time series generated from the OU process are filtered using a Gaussian kernel with a small standard deviation *σ*_{t} to obtain sufficiently differentiable time series ${\tilde{\mu}}_{\text{ext}}$ (due to the requirements of the spec_{2} model and the LN_{dos} model). The filtered realization ${\tilde{\mu}}_{\text{ext}}\left(t\right)$ is then used for all models to allow for a quantitative comparison of the different spike rate responses to the same input. The value of *σ*_{t} we use in this study effectively removes very large oscillation frequencies which are rarely observed, while lower frequencies [22] are passed.

The parameter space we explore covers large and small correlation times ${\tau}_{\text{ou}}^{\mu}$, strong and weak input mean $\overline{\mu}$ and standard deviation *σ*_{ext}, and for each of these combinations we consider an interval from small to large variation magnitudes _{μ}. The values of ${\tau}_{\text{ou}}^{\mu}$ and _{μ} determine how rapid and intense *μ*_{ext}(*t*) fluctuates.

We apply two performance measures, as in [21]. One is given by Pearson’s correlation coefficient,

$$\rho ({r}_{N},r)\u2254\frac{{\sum}_{k=1}^{M}\phantom{\rule{1pt}{0ex}}({r}_{N}\left({t}_{k}\right)-{\overline{r}}_{N})\phantom{\rule{1pt}{0ex}}(r\left({t}_{k}\right)-\overline{r})}{\sqrt{{\sum}_{k=1}^{M}\phantom{\rule{1pt}{0ex}}{({r}_{N}\left({t}_{k}\right)-{\overline{r}}_{N})}^{2}\phantom{\rule{1pt}{0ex}}{\sum}_{k=1}^{M}\phantom{\rule{1pt}{0ex}}{(r\left({t}_{k}\right)-\overline{r})}^{2}}},$$

(11)

between the (discretely given) spike rates of the aEIF population and each derived model with time averages ${\overline{r}}_{N}=1/M\phantom{\rule{1pt}{0ex}}{\sum}_{k=1}^{M}\phantom{\rule{1pt}{0ex}}{r}_{N}\left({t}_{k}\right)$ and $\overline{r}=1/M\phantom{\rule{1pt}{0ex}}{\sum}_{k=1}^{M}\phantom{\rule{1pt}{0ex}}r\left({t}_{k}\right)$ over a time interval of length *t*_{M} − *t*_{1}. For comparison, we also include the correlation coefficient between the aEIF population spike rate and the time-varying mean input, *ρ*(*r*_{N}, *μ*_{ext}). In addition, to assess absolute differences we calculate the root mean square (RMS) distance,

$${\text{d}}_{\text{RMS}}\left({r}_{N},r\right)\u2254\sqrt{\frac{1}{M}{\displaystyle \sum _{k=1}^{M}{\left({r}_{N}\left({t}_{k}\right)-r\left({t}_{k}\right)\right)}^{2}}},$$

(12)

where *M* denotes the number of elements of the respective time series (*r*_{N}, *r*).

We find that three of the four low-dimensional spike rate models (spec_{2}, LN_{exp}, LN_{dos}) very well reproduce the spike rate *r*_{N} of the aEIF neurons: for the LN_{exp} model *ρ* > 0.95 and for the spec_{2} and LN_{dos} models *ρ* 0.8 (each) over the explored parameter space, see Fig 2. Only the basic spectral model (spec_{1}) is substantially less accurate. Among the best models, the simplest (LN_{exp}) overall outperforms spec_{2} and LN_{dos}, in particular for fast and strong mean input variations. However, in the strongly mean-driven regime the best performing model is spec_{2}.

We observe that the performance of any of the spike rate models decreases (with model-specific slope) with (i) increasing variation strength _{μ} larger than a certain (small) value, and with (ii) smaller ${\tau}_{\text{ou}}^{\mu}$, i.e., faster changes of *μ*_{ext}. For small values of _{μ} fluctuations of *r*_{N}, which are caused by the finite aEIF population size *N* and do not depend on the fluctuations of *μ*_{ext}, deteriorate the performance measured by *ρ* (see also [21], p.13 right). This explains why *ρ* does not increase as _{μ} decreases (towards zero) for any of the models. Naturally, the FP model is the by far most accurate spike rate description in terms of both measures, correlation coefficient *ρ* and RMS distance. This is not surprising because the four low-dimensional models are derived from that (infinite-dimensional) representation. Thus, the performance of the FP system defines an upper bound on the correlation coefficient *ρ* and a lower bound on the RMS distance for the low-dimensional models.

In detail: for moderately fast changing mean input (large ${\tau}_{\text{ou}}^{\mu}$) the three models spec_{2}, LN_{exp} and LN_{dos} exhibit excellent reproduction performance with *ρ* > 0.95, and spec_{1} shows correlation coefficients of at least *ρ* = 0.9 (Fig 2A), which is substantially better than *ρ*(*r*_{N}, *μ*_{ext}). The small differences between the three top models can be better assessed from the RMS distance measure. For large input variance ${\sigma}_{\text{ext}}^{2}$ the two LN models perform best (cf. Fig 2A, top, and for an example, 2C). For weak input variance and large mean (small *σ*_{ext}, large $\overline{\mu}$) the spec_{2} model outperforms the LN models, unless the variation magnitude _{μ} is very large. For small mean $\overline{\mu}$, where transient activity is interleaved with periods of quiescence, the LN_{exp} model performs best, except for weak variations _{μ}, where LN_{dos} is slightly better (see Fig 2A, bottom).

Stronger differences in performance emerge when considering faster changes of the mean input *μ*_{ext}(*t*) (i.e., for small ${\tau}_{\text{ou}}^{\mu}$), see Fig 2B, and for examples, Fig 2C. The spec_{1} model again performs worst with *ρ* values even below the input/output correlation baseline *ρ*(*r*_{N}, *μ*_{ext}) for large mean input $\overline{\mu}$ (cf. Fig 2B, left). The spec_{1} spike rate typically decays too slowly (cf. Fig 2C). The three better performing models differ as follows: for large input variance and mean (large *σ*_{ext} and $\overline{\mu}$), where the spike rate response to the input is rather fast (cf. increased *ρ*(*r*_{N}, *μ*_{ext})), the performance of all three models in terms of *ρ* is very high, but the RMS distance measure indicates that LN_{exp} is the most accurate model (cf. Fig 2B, top). For weak mean input LN_{exp} is once again the top model while LN_{dos} and, especially noticeable, spec_{2} show a performance decline (see example in Fig 2C). For weak input variance (Fig 2A, bottom), where significant (oscillatory) excursions of the spike rates in response to changes in the mean input can be observed (see also Fig 1), we obtain the following benchmark contrast: for large mean drive $\overline{\mu}$ the spec_{2} model performs best, except for large variation amplitudes _{μ}, at which LN_{exp} is more accurate. Smaller mean input on the other hand corresponds to the most sensitive regime where periods of quiescence alternate with rapidly increasing and decaying spike rates. The LN_{exp} model shows the most robust and accurate spike rate reproduction in this setting, while LN_{dos} and spec_{2} each exhibit decreased correlation and larger RMS distances–spec_{2} even for moderate input variation intensities _{μ}. The slowness approximation underlying the spec_{2} model likely induces an error due to the fast external input changes in comparison with the rather slow intrinsic time scale by the dominant eigenvalue, ${\tau}_{\text{ou}}^{\mu}=5\phantom{\rule{1pt}{0ex}}\text{ms}$ vs. 1/|Re{λ_{1}}| ≈ 15 ms (cf. visualization of the spectrum in Sect. *Spectral models*). Note that for these weak inputs the distribution of the spike rate is rather asymmetric (cf. Fig 2B). Interestingly the LN_{dos} model performs worse than LN_{exp} for large mean input variations (i.e., large _{μ}) in general, and only slightly better for small input variance and mean input variations that are not too large and fast.

We would like to note that decreasing the Gaussian filter width *σ*_{t} to smaller values, e.g., fractions of a millisecond, can lead to a strong performance decline for the spec_{2} model because of its explicit dependence on first and second order time derivatives of the mean input.

Furthermore, we show how the adaptation parameters affect the reproduction performance of the different models in Fig 3. The adaptation time constant *τ*_{w} and spike-triggered adaptation increment *b* are varied simultaneously (keeping their product constant) such that the average spike rate and adaptation current, and thus the spiking regime, remain comparable for all parametrizations. As expected, the accuracy of the derived models decreases for faster adaptation current dynamics, due to the adiabatic approximation that relies on sufficiently slow adaptation (cf. Sect. Methods). Interestingly however, the performance of all reduced models (except spec_{1}) declines only slightly as the adaptation time constant decreases to the value of the membrane time constant (which means the assumption of separated time scales underlying the adiabatic approximation is clearly violated). This kind of robustness is particularly pronounced for input with large baseline mean $\overline{\mu}$ and small noise amplitude *σ*_{ext}, cf. Fig 3B.

For perfectly balanced excitatory and inhibitory synaptic coupling the contribution of presynaptic activity to the mean input *μ*_{syn} is zero by definition, but the input variance ${\sigma}_{\text{syn}}^{2}$ is always positively (linearly) affected by a presynaptic spike rate–even for a negative synaptic efficacy *J* (cf. Eq (2)). To assess the performance of the derived models in this scenario, but within the reference setting of an uncoupled population, we consider constant external mean drive *μ*_{ext} and let the variance ${\sigma}_{\text{ext}}^{2}\left(t\right)$ evolve according to a filtered OU process (such as that used for the mean input *μ*_{ext} in the previous section) with parameters $\overline{{\sigma}^{2}}$ and _{σ2} of the stationary normal distribution $\mathcal{N}(\overline{{\sigma}^{2}},{\vartheta}_{{\sigma}^{2}}^{2})$, correlation time ${\tau}_{\text{ou}}^{{\sigma}^{2}}$ and Gaussian filter standard deviation *σ*_{t} as before.

The results of two input parametrizations are shown in Fig 4. For large input mean *μ*_{ext} and rapidly varying variance ${\sigma}_{\text{ext}}^{2}\left(t\right)$ the spike rate response of the aEIF population is very well reproduced by the FP model and, to a large extent, by the spec_{2} model (cf. Fig 4A). This may be attributed to the fact that the latter model depends on the first two time derivatives of the input variance ${\sigma}_{\text{ext}}^{2}$. The LN models cannot well reproduce the rapid spike rate excursions in this setting, and the spec_{1} model performs worst, exhibiting time-lagged spike rate dynamics compared to *r*_{N}(*t*) which leads to a very small value of correlation coefficient *ρ* (below the input/output correlation baseline $\rho ({r}_{N},{\sigma}_{\text{ext}}^{2})$). For smaller mean input *μ*_{ext} and moderately fast varying variance ${\sigma}_{\text{ext}}^{2}\left(t\right)$ (larger correlation time ${\tau}_{\text{ou}}^{{\sigma}^{2}}$) the fluctuating aEIF population spike rate is again nicely reproduced by the FP model while the rate response of the spec_{2} model exhibits over-sensitive behavior to changes in the input variance, as indicated by the large RMS distance (see Fig 4B). This effect is even stronger for faster variations, i.e., smaller ${\tau}_{\text{ou}}^{{\sigma}^{2}}$ (cf. supplementary visualization S1 Fig). The LN models perform better in this setting, and the spec_{1} model (again) performs worst in terms of correlation coefficient *ρ* due to its time-lagged spike rate response.

It should be noted that the lowest possible value of the input standard deviation, i.e., *σ*_{ext} (plus a nonnegative number in case of recurrent input) cannot be chosen completely freely but must be large enough ($\gtrsim \text{\hspace{0.17em}}0.5\phantom{\rule{1pt}{0ex}}\text{mV}/\sqrt{\text{ms}}$) for our parametrization. This is due to theoretical reasons (Fokker-Planck formalism) and practical reasons (numerics for Fokker-Planck solution and for calculation of the derived quantities, such as *r*_{∞}).

To demonstrate the applicability of the low-dimensional models for network analyses we consider a recurrently coupled population of aEIF neurons that produces self-sustained network oscillations by the interplay of strong excitatory feedback and spike-triggered adaptation or, alternatively, by delayed recurrent synaptic inhibition [16, 23]. The former oscillation type is quite sensitive to changes in input, adaptation and especially coupling parameters for the current-based type of synaptic coupling considered here and due to lack of (synaptic) inhibition and refractoriness. For example, a small increase in coupling strength can lead to a dramatic (unphysiologic) increase in oscillation amplitude because of strong recurrent excitation. Hence we consider a difficult setting here to evaluate the reduced spike rate models–in particular, when the network operates close to a bifurcation.

In Fig 5A we present two example parametrizations from a region (in parameter space) that is characterized by stable oscillations. This means the network exhibits oscillatory spike rate dynamics for constant external input moments *μ*_{ext} and ${\sigma}_{\text{ext}}^{2}$. The derived models reproduce the limit cycle behavior of the aEIF network surprisingly well, except for small frequency and amplitude deviations (FP, spec_{2}, LN_{dos}, LN_{exp}) and larger frequency mismatch (spec_{1}), see Fig 5A, top. For weaker input moments and increased spike-triggered adaptation strength the network is closer to a Hopf bifurcation [16, 23]. It is, therefore, not surprising that the differences in oscillation period and amplitude are more prominent (cf. Fig 5A, bottom). The bifurcation point of the LN_{exp} model is slightly shifted, shown by the slowly damped oscillatory convergence to a fixed point. This suggests that the bifurcation parameter value of each of the derived models is not far from the true critical parameter value of the aEIF network but can quantitatively differ (slightly) in a model-dependent way.

The second type of oscillation is generated by delayed synaptic inhibition [22] and does not depend on the (neuronal) inhibition that is provided by an adaptation current. To demonstrate this independence the adaptation current was disabled (by setting the parameters *a* = *b* = 0) for the two respective examples that are shown in Fig 5B. Similarly as for the previous oscillation type, the low-dimensional models (except spec_{1}) reproduce the spike rate limit cycle of the aEIF network surprisingly well, in particular for weak external input (see Fig 5B, top). For larger external input and stronger inhibition with shorter delay the network operates close to a Hopf bifurcation, leading to larger differences in oscillation amplitude and frequency in a model-dependent way (Fig 5B, bottom). Note that the intermediate (Fokker-Planck) model very well reproduces the inhibition-based type of oscillation which demonstrates the applicability of the underlying mean-field approximation. We would also like to note that enabling the adaptation current dynamics (only) leads to decreased average spike rates but does not affect the reproduction accuracy.

We would like to emphasize that the previous comprehensive evaluations for an uncoupled population provide a deeper insight on the reproduction performance–also for a recurrent network–than the four examples shown here, as explained in the Sect. *Performance for variations of the mean input*. For example, the (improved) reproduction performance for increased input variance in the uncoupled setting (cf. Fig 2) informs about the reproduction performance for networks of excitatory and inhibitory neurons that are roughly balanced, i.e., where the overall input mean is rather small compared to the input standard deviation.

We have developed efficient implementations of the derived models using the Python programming language and by employing the library Numba for low-level machine acceleration [24]. These include: (i) the numerical integration of the Fokker-Planck model using an accurate finite volume scheme with implicit time discretization (cf. Sect. Methods), (ii) the parallelized precalculation of the quantities required by the low-dimensional spike rate models and (iii) the time integration of the latter models, as well as example scripts demonstrating (i)–(iii). The code is available as open source software under a free license at GitHub: https://github.com/neuromethods/fokker-planck-based-spike-rate-models

With regards to computational cost, summarizing the results of several aEIF network parametrizations, the duration to generate population activity time series for the low-dimensional spike rate models is usually several orders of magnitude smaller compared to numerical simulation of the original aEIF network and a few orders of magnitude smaller in comparison to the numerical solution of the FP model. For example, considering a population of 50,000 coupled neurons with 2% connection probability, a single simulation run of 5 s and the same integration time step across the models, the computation times amounted to 1.1–3.6 s for the low-dimensional models (with order–fast to slow–LN_{exp}, spec_{1}, LN_{dos}, spec_{2}), about 100 s for the FP model and roughly 1500 s for the aEIF network simulation on a dual-core laptop computer. The time difference to the network simulation substantially increases with the numbers of neurons and connections, and with spiking activity within the network due to the extensive propagation of synaptic events. Note that the speedup becomes even more pronounced with increasing number of populations, where the runtimes of the FP model and the aEIF network simulation scale linearly and the low-dimensional models show a sublinear runtime increase due to vectorization of the state variables representing the different populations.

The derived low-dimensional (ODE) spike rate models are very efficient to integrate given that the required input-dependent parameters are available as precalulated look-up quantities. For the grids used in this contribution, the precomputation time was 40 min. for the cascade (LN_{exp}, LN_{dos}) models and 120 min. for the spectral (spec_{1}, spec_{2}) models, both on a hexa-core desktop computer. The longer calculation time for the spectral models was due to the finer internal grid for the mean input (see S1 Text).

Note that while the time integration of the spec_{2} model is on the same order as for the other low-dimensional models its implementation complexity is larger because of the many quantities it depends on, cf. Eqs (63)–(66).

In this contribution we have developed four low-dimensional models that approximate the spike rate dynamics of coupled aEIF neurons and retain all parameters of the underlying model neurons. These simple spike rate models were derived in two different ways from a Fokker-Planck PDE that describes the evolving membrane voltage distribution in the mean-field limit of large networks, and is complemented by an ODE for the population-averaged slow adaptation current. Two of the reduced spike rate models (spec_{1} and spec_{2}) were obtained by a truncated spectral decomposition of the Fokker-Planck operator assuming vanishingly slow (for spec_{1}) or moderately slow (for spec_{2}) changes of the input moments. The other two reduced models (LN_{exp} and LN_{dos}) are described by a cascade of linear filters (one for the input mean and another for its standard deviation) and a nonlinearity which were derived from the Fokker-Planck equation, and subsequently the filters were semi-analytically approximated. Our approaches build upon [18–20] as well as [21], and extend those methods for adaptive nonlinear integrate-and-fire neurons that are sparsely coupled with distributed delays (cf. Sect. Methods).

We have compared the different spike rate representations for a range of biologically plausible input statistics and found that three of the reduced models (spec_{2}, LN_{exp} and LN_{dos}) accurately reproduce the spiking activity of the underlying aEIF population while one model (spec_{1}) shows the least accuracy. Among the best models, the simplest (LN_{exp}) was the most robust and (somewhat surprisingly) overall outperformed spec_{2} and LN_{dos}–especially in the sensitive regime of rapidly changing sub- and suprathreshold mean drive and in general for rapid and strong input variations. The LN_{exp} model did not exhibit exaggerated deflections in that regime as compared to the other two models. This result is likely due to the importance of the quantitatively correct decay time of the filter for the mean input in the LN_{exp} model, while the violations of the slowness assumptions for the spec_{2} and LN_{dos} models seem more harmful in this regime. In the strongly mean-driven regime, however, the best performing model was spec_{2} for variations both in the mean drive (as long as those variations are not too strong and fast) and for variations of the input variance.

We have also demonstrated that the low-dimensional models well reproduce the dynamics of recurrently coupled aEIF populations in terms of asynchronous states (see Fig 1) and spike rate oscillations (cf. Fig 5), where mild deviations at critical (bifurcation) parameter values are expected due to the approximative nature of the model reduction.

The computational demands of the low-dimensional models are very modest in comparison to the aEIF network and also to the integration of the Fokker-Planck PDE, for which we have developed a novel finite volume discretization scheme. We would like to emphasize that any change of a parameter value for input, coupling or adaptation current does not require renewed precomputations. To facilitate the application of the presented models we have made available implementations that precompute all required quantities and numerically integrate the derived low-dimensional spike rate models as well as the Fokker-Planck equation, together with example (Python) scripts, as open source software.

Since the derived models are formulated in terms of simple ODEs, they allow to conveniently perform linear stability analyses, e.g., based on the eigenvalues of the Jacobian matrix of the respective vector field. In this way network states can be rapidly characterized by quantifying the bifurcation structure of the population dynamics–including regions of the parameter space where multiple fixed points and/or limit cycle attractors co-exist. For a characterization of stable network states by numerical continuation and an assessment of their controllability through neuromodulators using the LN_{exp} model see [23] ch. 4.2 and [25]. Furthermore, the low-dimensional models are well suited to be employed in large neuronal networks of multiple populations for efficient simulations of population-averaged activity time series. Overall, the LN_{exp} model seems a good candidate for that purpose considering its accuracy and robustness, as well as its computational and implementational simplicity.

We considered a homogeneous population of neurons in the sense that the parameter values across model neurons are identical except those for synaptic input. Thereby we assume that neurons with similar dynamical properties can be grouped into populations [3]. Heterogeneity is incorporated by distributed synaptic delays, by sparse random coupling, and by fluctuating external inputs for each neuron. The (reduced) population models further allow for heterogeneous synaptic strengths that are sampled from a Gaussian distribution and can be included in a straightforward way [16, 26] (see also Sect. Methods). Distributed values for other parameters (of the isolated model neurons within the same population) are currently not supported.

The presented mean-field network model can be easily adjusted for multiple populations. In this case we obtain a low-dimensional ODE for each population and the overall synaptic moments for population *k* become

$${\mu}_{\text{syn},k}={\mu}_{\text{ext},k}\left(t\right)+\sum _{l}\phantom{\rule{1pt}{0ex}}{J}_{kl}{K}_{kl}{r}_{d,kl}\left(t\right),\phantom{\rule{2.em}{0ex}}{\sigma}_{\text{syn},k}^{2}={\sigma}_{\text{ext},k}^{2}\left(t\right)+\sum _{l}\phantom{\rule{1pt}{0ex}}{J}_{kl}^{2}{K}_{kl}{r}_{d,kl}\left(t\right),$$

(13)

where *J*_{kl} is the synaptic strength for the *K*_{kl} neurons from population *l* targeting neurons from population *k* and *r*_{d,kl} is the delayed spike rate of population *l* affecting population *k* (cf. Eq (2)). For each pair of coupled populations we may consider identical or distributed delays (using distributions from the exponential family) as well as identical or distributed synaptic strengths (sampled from a Gaussian distribution).

Here we described synaptic interaction by delayed (delta) current pulses with delays sampled from an exponential distribution. This description leads to a fluctuating overall synaptic input current with white noise characteristics. Interestingly, for the mean-field dynamics this setting is very similar to considering exponentially decaying synaptic currents with a decay constant that matches that of the delay distribution, although the overall synaptic input current is a colored noise process in that case, see [27] and, for an intuitive explanation [28].

A conductance-based model of synaptic coupling can also be considered in principle [16, 29], which results in a multiplicative noise process for the overall synaptic input. This, however, would in general impede the beneficial concept of precalculated “look-up” quantities that are unaffected by the input and coupling parameters.

It should be noted that most current- or conductance-based models of synaptic coupling (including the one considered here) can produce unphysiologically large amounts of synaptic current in case of high presynaptic activity, unless the coupling parameters are carefully tuned. This problem can be solved, for example, by considering a (more realistic) model of synaptic coupling based on [30], from which activity-dependent coupling terms can be derived for the mean-field and reduced population models [23] ch. 4.2. Using that description ensures robust simulation of population activity time series without having to fine-tune the coupling parameter values, which is particularly useful for multi-population network models. In this contribution though we used for simplicity a basic synaptic coupling model that has frequently been applied in the mean-field literature.

The Gaussian stochastic process driving the individual neurons could also be substituted by colored noise, which would lead to a Fokker-Planck model with increased dimensionality [31]. However, this would require more complex and computationally expensive numerical schemes not only to solve that model but also for the different dimension reduction approaches.

To derive low-dimensional models of population activity we approximated the adaptation current by its population average, justified by its slow dynamics compared to the other time scales of the system. This approximation is equivalent to a first order moment closure method [17]. In case of a faster adaptation time scale the approximation may be improved by considering second and higher order moments [17, 32].

The mean-field models presented here can well reproduce the dynamics of population-averaged state variables (that is, spike rate, mean membrane voltage, and mean adaptation current) for large populations (*N* → ∞ in the derivation). Fluctuations of those average variables due to the finite size of neuronal populations, however, are not captured. Hence, it would be interesting to extend the mean-field models so as to reproduce these (so-called) finite size effects, for example, by incorporating an appropriate stochastic process [18] or using concepts from [15].

For uncoupled EIF populations (without an adaptation current) and constant input standard deviation it has been shown that the LN cascade approximation performs well for physiological ranges of amplitude and time scale for mean input variations [21]. Our results for the cascade models are consistent with [21], but the performance is substantially improved for the sensitive low (baseline) input regime (LN_{exp} and LN_{dos}, also in absence of adaptation), and damped oscillatory behavior (including over- and undershoots) is accounted for by the LN_{dos} model.

To achieve these improvements we semi-analytically fit the linear filters derived from the Fokker-Planck equation using exponential and damped oscillator functions considering a range of input frequencies. The approximation can be further improved by using more complex functions, such as a damped oscillator with two time scales. That, however, can lead to less robustness (i.e., undesired model behavior) for rapid and strong changes of the input moments (cf. Sect. Methods).

LN cascade models are frequently applied in neuroscience to describe population activity, and the model components are often determined from electrophysiological recordings using established techniques. The methodology presented here contributes to establishing quantitative links between networks of spiking neurons, a mesoscopic description of population activity and recordings at the population level.

Here we provide a new numerical solver for the eigenvalue problem of the Fokker-Planck operator and its adjoint. This allows to compute the full spectrum together with associated eigenfunctions and is applicable to nonlinear integrate-and-fire models, extending [18, 19, 33].

Using that solver the spec_{2} model, which is based on two eigenvalues, can be further improved by interpolating its coefficients, Eqs (63)–(66), around the double eigenvalues at the spectrum’s real-to-complex transition. This interpolation would effectively smooth the quantities—e.g., preventing the jumps and kinks that are present in the visualization of Sect. *Spectral models*—and is expected to increase the spike rate reproduction accuracy (particularly for weak mean input) beyond what was reported in this contribution.

The spec_{2} model can also be extended to yield a third order ODE with everywhere smooth coefficients by considering an additional eigenvalue (cf. Sect. *Remarks on the spectrum*).

Moreover, the spec_{2} model, and more generally the whole spectral decomposition approach, can be extended to account for a refractory period in the presence of time-varying total input moments, e.g., by building upon previous attempts [18, 34, 35].

Furthermore, it could be beneficial to explicitly quantify the approximation error due to the slowness assumption that underlies the spec_{2} model by integration of the (original) spectral representation of the Fokker-Planck model.

Both reduced spectral models allow for a refined description of the mean adaptation current dynamics, cf. Eq (1), by replacing the mean membrane voltage *V* with its steady-state value *V*_{∞}, using that the membrane voltage distribution is available through the eigenfunctions of the Fokker-Planck operator.

The numerical eigenvalue solver can be extended in a straightforward way to yield quantities that are required by the original spectral representation of the Fokker-Planck model and by the corresponding stochastic equation for finite population size [18].

In addition to the work we build upon [18–21] (cf. Sect. Methods) there are a few other approaches to derive spike rate models from populations of spiking neurons. Some methods also result in an ODE system, taking into account (slow) neuronal adaptation [17, 26, 36–38] or disregarding it [39]. The settings differ from the work presented here in that (i) the intrinsic neuronal dynamics are adiabatically neglected [17, 26, 36, 37], (ii) only uncoupled populations [38] or all-to-all connected networks [17, 36, 39] are assumed in contrast to sparse connectivity, and (iii) (fixed) heterogeneous instead of fluctuating input is considered [39]. Notably, these previous methods yield rather qualitative agreements with the underlying spiking neuron population activity except for [39] where an excellent quantitative reproduction for (non-adaptive) quadratic integrate-and-fire oscillators with quenched input randomness is reported.

Other approaches yield mesoscopic representations of population activity in terms of model classes that are substantially less efficient to simulate and more complicated to analyze than low-dimensional ODEs [14–17, 40–42]. The spike rate dynamics in these models has been described (i) by a rather complex ODE system that depends on a stochastic jump process derived for integrate-and-fire neurons without adaptation [40], (ii) by PDEs for recurrently connected aEIF [16] or Izhikevich [17] neurons, (iii) by an integro-PDE with displacement for non-adaptive neurons [42] or (iv) by integral equations that represent the (mean) activity of coupled phenomenological spiking neurons without [41] and with adaptation [14, 15].

Furthermore, the stationary condition of a noise-driven population of adaptive EIF neurons [32, 43, 44] and the first order spike rate response to weak input modulations [43, 44] have been analyzed using the Fokker-Planck equation. Ref. [32] also considered a refined approximation of the (purely spike-triggered) adaptation current including higher order moments.

It may be interesting for future studies to explore ways to extend the presented methods and relax some of the underlying assumptions, in particular, considering (i) the diffusion approximation (via shot noise input, e.g., [45, 46]), (ii) the Poisson assumption (e.g., using the concept from [47] in combination with results from [48]) and (iii) (noise) correlations (see, e.g., [49]).

Here we present all models in detail—the aEIF network (*ground truth*), the mean-field FP system (*intermediate model*) and the low-dimensional models: spec_{1}, spec_{2}, LN_{exp}, LN_{dos}—including step-by-step derivations and essential information on the respective numerical solution methods. An implementation of these models using Python is made available at GitHub: https://github.com/neuromethods/fokker-planck-based-spike-rate-models

We consider a large (homogeneous) population of *N* synaptically coupled aEIF model neurons [5]. Specifically, for each neuron (*i* = 1, …, *N*), the dynamics of the membrane voltage *V*_{i} is described by

$$C\frac{d{V}_{i}}{dt}={I}_{\mathrm{L}}\left({V}_{i}\right)+{I}_{\text{exp}}\left({V}_{i}\right)-{w}_{i}+{I}_{\text{syn},i}\left(t\right),$$

(14)

where the capacitive current through the membrane with capacitance *C* equals the sum of three ionic currents and the synaptic current *I*_{syn,i}. The ionic currents consist of a linear leak current *I*_{L}(*V*_{i}) = −*g*_{L}(*V*_{i} − *E*_{L}) with conductance *g*_{L} and reversal potential *E*_{L}, a nonlinear term *I*_{exp}(*V*_{i}) = *g*_{L} Δ_{T} exp((*V*_{i} − *V*_{T})/Δ_{T}) that approximates the rapidly increasing Na^{+} current at spike initiation with threshold slope factor Δ_{T} and effective threshold voltage *V*_{T}, and the adaptation current *w*_{i} which reflects a slowly deactivating K^{+} current. The adaptation current evolves according to

$${\tau}_{w}\frac{d{w}_{i}}{dt}=a({V}_{i}-{E}_{w})-{w}_{i},$$

(15)

with adaptation time constant *τ*_{w}. Its strength depends on the subthreshold membrane voltage via conductance *a*. *E*_{w} denotes its reversal potential. When *V*_{i} increases beyond *V*_{T}, it diverges to infinity in finite time due to the exponentially increasing current *I*_{exp}(*V*_{i}), which defines a spike. In practice, however, the spike is said to occur when *V*_{i} reaches a given value *V*_{s}—the spike voltage. The downswing of the spike is not explicitly modeled; instead, when *V*_{i} ≥ *V*_{s}, the membrane voltage *V*_{i} is instantaneously reset to a lower value *V*_{r}. At the same time, the adaptation current *w*_{i} is incremented by a value of parameter *b*, which implements suprathreshold (spike-dependent) activation of the adaptation current.

Immediately after the reset, *V*_{i} and *w*_{i} are clamped (i.e., remain constant) for a short refractory period *T*_{ref}, and subsequently governed again by Eqs (14) and (15). At the end of the Methods section we describe how (optionally) a spike shape can be included in the aEIF model, together with the associated small changes for the models derived from it.

To complete the network model the synaptic current in Eq (14) needs to be specified: for each cell it is given by the sum of recurrent and external input, *I*_{syn,i} = *I*_{rec,i}(*t*) + *I*_{ext,i}(*t*). Recurrent synaptic input is received from *K* other neurons of the network, that are connected in a sparse (*K* *N*) and uniformly random way, and is modeled by

$${I}_{\text{rec},i}=C\phantom{\rule{1pt}{0ex}}{\displaystyle \sum _{j}\phantom{\rule{1pt}{0ex}}{J}_{ij}}\phantom{\rule{1pt}{0ex}}{\displaystyle \sum _{{t}_{j}}\phantom{\rule{1pt}{0ex}}\delta}(t-{t}_{j}-{d}_{ij}),$$

(16)

where *δ* denotes the Dirac delta function. Every spike by one of the *K* presynaptic neurons with indices *j* and spike times *t*_{j} causes a postsynaptic membrane voltage jump of size *J*_{ij}. The coupling strength is positive (negative) for excitation (inhibition) and of small magnitude. Here it is chosen to be constant, i.e., *J*_{ij} = *J*. Each of these membrane voltage deflections occur after a time delay *d*_{ij} that takes into account (axonal and dendritic) spike propagation times and is sampled (independently) from a probability distribution *p*_{d}. In this work we use exponentially distributed delays, i.e., *p*_{d}(*τ*) = exp(−*τ*/*τ*_{d})/*τ*_{d} (for *τ* ≥ 0) with mean delay *τ*_{d}.

The second type of synaptic input is a fluctuating current generated from network-external neurons,

(17)

with time-varying moments *μ*_{ext} and ${\sigma}_{\text{ext}}^{2}$, and unit Gaussian white noise process *ξ*_{ext,i}. The latter is uncorrelated with that of other neurons *j* ≠ *i*, i.e., *ξ*_{ext,i}(*t*)*ξ*_{ext,j}(*t* + *τ*) = *δ*(*τ*)*δ*_{ij}, where · denotes expectation (w.r.t. the joint ensemble of noise realizations at times *t* and *t* + τ) and *δ*_{ij} is the Kronecker delta. This external current, for example, accurately approximates the input generated from a large number of independent Poisson neurons that produce instantaneous postsynaptic potentials of small magnitude, cf. [48].

The spike rate *r*_{N} of the network is defined as the population-averaged number of emitted spikes per time interval [*t*, *t* + Δ*T*],

$${r}_{N}(t)=\frac{1}{N}\phantom{\rule{1pt}{0ex}}{\displaystyle \sum _{i=1}^{N}\phantom{\rule{1pt}{0ex}}\frac{1}{\Delta T}\phantom{\rule{1pt}{0ex}}}{\displaystyle {\int}_{t}^{t+\Delta T}\phantom{\rule{1pt}{0ex}}{\displaystyle \sum _{{t}_{i}}\phantom{\rule{1pt}{0ex}}\delta}}(s-{t}_{i})ds,$$

(18)

where the interval size Δ*T* is practically chosen small enough to capture the dynamical structure and large enough to yield a comparably smooth time evolution for a finite network, i.e., *N* < ∞.

We chose values for the neuron model parameters to describe cortical pyramidal cells, which exhibit “regular spiking” behavior and spike frequency adaptation [7, 50, 51]. For the complete parameter specification see Table 1.

All network simulations were performed using the Python software BRIAN2 [52, 53] with C++ code generation enabled for efficiency. The aEIF model Eqs (14) and (15) were discretized using the Euler-Maruyama method with equidistant time step Δ*t* and initialized with *w*_{i}(0) = 0 and *V*_{i}(0) that is (independently) sampled from a Gaussian initial distribution *p*_{0}(*V*) with mean *V*_{r} − *δV* and standard deviation *δV*/2 where *δV* = *V*_{T} − *V*_{r}. Note that the models derived in the following Sects. do not depend on this particular initial density shape but allow for an arbitrary (density) function *p*_{0}.

The time scales of (slow) K^{+} channel kinetics which are effectively described by the adaptation current *w*_{i}, cf. Eq (15), are typically much larger than the faster membrane voltage dynamics modeled by Eq (14), i.e., *τ*_{w} *C*/*g*_{L} [54–57]. This observation justifies to replace the individual adaptation current *w*_{i} in Eq (14) by its average across the population, ${\langle w\rangle}_{N}=1/N\phantom{\rule{1pt}{0ex}}{\sum}_{i=1}^{N}{w}_{i}\left(t\right)$, in order to reduce computational demands and enable further analysis. The mean adaptation current is then governed by [16, 26, 48, 58]

$$\frac{d{\u27e8w\u27e9}_{N}}{dt}=\frac{a({\u27e8V\u27e9}_{N}-{E}_{w})-{\u27e8w\u27e9}_{N}}{{\tau}_{w}}+b\phantom{\rule{0.166667em}{0ex}}{r}_{N}\left(t\right),$$

(19)

where *V*_{N} denotes the time-varying population average of the membrane voltage of non-refractory neurons.

The dynamics of the population-averaged adaptation current reflecting the non-refractory proportion of neurons are well captured by Eq (19) as long as *T*_{ref} is small compared to *τ*_{w}. In this (physiologically plausible) case *w*_{N} from Eq (19) can be considered equal to the average adaptation current over the refractory proportion of neurons [16, 48].

For large networks (*N* → ∞) the recurrent input can be approximated by a mean part with additive fluctuations, ${I}_{\text{rec}}{}_{,\text{i}}/C\approx JK{r}_{d}\left(t\right)+J\sqrt{K{r}_{d}\left(t\right)}{\xi}_{\text{rec},i}\left(t\right)$ with delayed spike rate

(20)

i.e., the spike rate convolved with the delay distribution, and unit white Gaussian noise process *ξ*_{rec,i} that is uncorrelated to that of any other neuron [16, 18, 26, 34].

The step is valid under the assumptions of (i) sufficiently many incoming synaptic connections (*K* 1) with small enough weights |*J*_{ij}| in comparison with *V*_{T} − *V*_{r} and sufficient presynaptic activity (diffusion approximation) (ii) that neuronal spike trains can be approximated by independent Poisson processes (Poisson assumption) and (iii) that the correlations between the fluctuations of synaptic inputs for different neurons vanish (mean-field limit). The latter assumption is fulfilled by sparse and uniformly random synaptic connectivity, but also when synaptic strengths *J*_{ij} and delays *d*_{ij} are independently distributed (in case of less sparse or random connections) [18].

This approximation of the recurrent input allows to replace the overall synaptic current in Eq (14) by *I*_{syn,i} = *C*[*μ*_{syn}(*t*, *r*_{d}) + *σ*_{syn}(*t*, *r*_{d})*ξ*_{i}(*t*)] with overall synaptic moments

$${\mu}_{\text{syn}}={\mu}_{\text{ext}}\left(t\right)+JK{r}_{d}\left(t\right),\phantom{\rule{2.em}{0ex}}{\sigma}_{\text{syn}}^{2}={\sigma}_{\text{ext}}^{2}\left(t\right)+{J}^{2}K{r}_{d}\left(t\right),$$

(21)

and (overall) unit Gaussian white noise *ξ*_{i} that is uncorrelated to that of any other neuron. Here we have used that external *I*_{ext,i} and recurrent synaptic current *I*_{rec,i} are independent from each other.

The resulting mean-field dynamics of the membrane voltage is given by

$$\frac{d{V}_{i}}{dt}=\frac{{I}_{\mathrm{L}}\left({V}_{i}\right)+{I}_{\text{exp}}\left({V}_{i}\right)-\u27e8w\u27e9}{C}+{\mu}_{\text{syn}}(t,{r}_{d})+{\sigma}_{\text{syn}}(t,{r}_{d}){\xi}_{i}\left(t\right),$$

(22)

and corresponds to a McKean-Vlasov type of equation with distributed delays [59] and discontinuity due to the reset mechanism [60] that complements the dynamics of *V*_{i} as before. The population-averaged adaptation current *w* = lim_{N → ∞}*w*_{N} is governed by

$$\frac{d\u27e8w\u27e9}{dt}=\frac{a(\u27e8V\u27e9-{E}_{w}\phantom{{}^{2}})-\u27e8w\u27e9}{{\tau}_{w}}+b\phantom{\rule{0.166667em}{0ex}}r\left(t\right),$$

(23)

with mean membrane voltage (of non-refractory neurons), *V* = lim_{N→∞}*V*_{N}, and spike rate *r* = lim_{N→∞,Δt→0}
*r*_{N}(*t*).

Remarks: Instead of exponentially distributed synaptic delays we may also consider other continuous densities *p*_{d}, identical delays, *p*_{d}(*τ*) = *δ*(*τ* − *d*) with *d* > 0, or no delays at all, *p*_{d}(*τ*) = *δ*(*τ*). Instead of identical synaptic strengths one may also consider strengths *J*_{ij} that are drawn independently from a normal distribution with mean *J*_{m} and variance *J*_{v} instead, in which case the overall synaptic moments become *μ*_{syn} = *μ*_{ext}(*t*) + *J*_{m}*Kr*_{d}(*t*) and ${\sigma}_{\text{syn}}^{2}={\sigma}_{\text{ext}}^{2}\left(t\right)+({J}_{m}^{2}+{J}_{v})K{r}_{d}\left(t\right)$, cf. [16, 26].

In the membrane voltage evolution, Eq (22), individual neurons are exchangeable as they are described by the same stochastic equations and are coupled to each other exclusively through the (delayed) spike rate via the overall synaptic moments *μ*_{syn} and ${\sigma}_{\text{syn}}^{2}$. Therefore, the adiabatic and mean-field approximations allow us to represent the collective dynamics of a large network by a (1+1-dimensional) Fokker-Planck equation [16, 18, 26, 34],

$$\frac{\partial}{\partial t}p(V,t)+\frac{\partial}{\partial V}{q}_{p}(V,t)=0\phantom{\rule{2.em}{0ex}}\text{for}\phantom{\rule{1.em}{0ex}}V\in (-\infty ,{V}_{\mathrm{s}}],\phantom{\rule{0.277778em}{0ex}}t>0,$$

(24)

which describes the evolution of the probability density *p*(*V*, *t*) to find a neuron in state *V* at time *t* (in continuity form). The probability flux is given by

$${q}_{p}(V,t)=(\frac{{I}_{\mathrm{L}}\left(V\right)+{I}_{\text{exp}}\left(V\right)}{C}+{\mu}_{\text{tot}}\left(t\right))\phantom{\rule{1pt}{0ex}}p(V,t)-\frac{{\sigma}_{\text{tot}}^{2}\left(t\right)}{2}\frac{\partial}{\partial V}p(V,t),$$

(25)

with *total* input mean and standard deviation,

(26)

(27)

Note that the mean adaptation current (simply) subtracts from the synaptic mean in the drift term, cf. Eq (22).

The mean adaptation current evolves according to Eq (23) with time-dependent mean membrane voltage (of the non-refractory neurons)

$$\u27e8V\u27e9=\frac{{\int}_{-\infty}^{{V}_{\mathrm{s}}}vp(v,t)dv}{{\int}_{-\infty}^{{V}_{\mathrm{s}}}p(v,t)dv}.$$

(28)

The spike rate *r* is obtained by the probability flux through *V*_{s},

(29)

To account for the reset condition of the aEIF neuron dynamics and ensuring that probability mass is conserved, Eq (24) is complemented by the reinjection condition,

$${q}_{p}({V}_{\mathrm{r}}^{+},t)-{q}_{p}({V}_{\mathrm{r}}^{-},t)={q}_{p}({V}_{\mathrm{s}},t-{T}_{\text{ref}}),$$

(30)

where ${q}_{p}\left({V}_{\mathrm{r}}^{+}\right)\u2254{lim}_{V\searrow {V}_{\mathrm{r}}}{q}_{p}\left(V\right)$ and ${q}_{p}\left({V}_{\mathrm{r}}^{-}\right)\u2254{lim}_{V\nearrow {V}_{\mathrm{r}}}{q}_{p}\left(V\right)$, an absorbing boundary at *V*_{s},

(31)

and a natural (reflecting) boundary condition,

$$\underset{V\to -\infty}{lim}{q}_{p}(V,t)=0.$$

(32)

Together with the initial membrane voltage distribution *p*(*V*, 0) = *p*_{0}(*V*) and mean adaptation current *w*(0) = 0 the Fokker-Planck mean-field model is now completely specified.

Note that *p*(*V*, *t*) only reflects the proportion of neurons which are not refractory at time *t*, given by $P\left(t\right)={\int}_{-\infty}^{{V}_{\mathrm{s}}}p(v,t)dv\phantom{\rule{4pt}{0ex}}=1-{\int}_{t-{T}_{\text{ref}}}^{t}r\left(s\right)ds$ (<1 for *T*_{ref} > 0 and *r*(*t*) > 0). The total probability density that the membrane voltage is *V* at time *t* is given by *p*(*V*, *t*) + *p*_{ref}(*V*, *t*) with refractory density *p*_{ref}(*V*, *t*) = (1 − *P*(*t*)) *δ*(*V* − *V*_{r}). At the end of the Methods section we describe how an (optional) spike shape extension for the aEIF model changes the calculation of *p*_{ref} and *V*.

In practice we consider a finite reflecting lower barrier *V*_{lb} instead of negative infinite for the numerical solution (next section) and for the low-dimensional approximations of the Fokker-Planck PDE (cf. sections below). *V*_{lb} is chosen sufficiently small in order to not distort the free diffusion of the membrane voltage for values below the reset, i.e., *V*_{lb} *V*_{r}. The density *p*(*V*, *t*) is then supported on [*V*_{lb}, *V*_{s}] for each time *t*, and in all expressions above *V* → −∞ is replaced by *V*_{lb}.

In this work we focus on low-dimensional approximations of the FP model. To obtain a reference for the reduced models it is, however, valuable to solve the (full) FP system, Eqs (24)–(32). Here we outline an accurate and robust method of solution that exploits the linear form of the FP model in contrast to previously described numerical schemes [61, 62] which both require rather small time steps due to the steeply increasing exponential current *I*_{exp} in the flux *q*_{p} close to the spike voltage *V*_{s}.

We first discretize the (finite) domain [*V*_{lb}, *V*_{s}] into *N*_{V} equidistant grid cells $[{V}_{m-\frac{1}{2}},{V}_{m+\frac{1}{2}}]$ with centers *V*_{m} (*m* = 1, …, *N*_{V}) that satisfy *V*_{1} < *V*_{2} < < *V*_{NV}, where *V*_{½} = *V*_{lb} and *V*_{NV+½} = *V*_{s} are the outmost cell borders. Within each cell the numerical approximation of *p*(*V*, *t*) is assumed to be constant and corresponds to the average value denoted by *p*(*V*_{m}, *t*). Integrating Eq (24) combined with Eq (30) over the volume of cell *m*, and applying the divergence theorem, yields

$$\frac{\partial}{\partial t}p({V}_{m},t)=\frac{{q}_{p}({V}_{m-\frac{1}{2}},t)-{q}_{p}({V}_{m+\frac{1}{2}},t)}{\Delta V}+{\delta}_{m{m}_{\mathrm{r}}}\frac{1}{\Delta V}{q}_{p}({V}_{{N}_{V}+\frac{1}{2}},t-{T}_{\text{ref}}),$$

(33)

where Δ*V* is the grid spacing and *m*_{r} corresponds to the index of the cell that contains the reset voltage *V*_{r}. To solve Eq (33) forward in time the fluxes at the borders of each cell need to be approximated. Since the Fokker-Planck PDE belongs to the class of drift-diffusion equations this can be accurately achieved by the first order Scharfetter-Gummel flux [63, 64],

$${q}_{p}({V}_{m+\frac{1}{2}},t)={v}_{m+\frac{1}{2}}\phantom{\rule{0.166667em}{0ex}}\frac{p({V}_{m},t)-p({V}_{m+1},t)\phantom{\rule{0.166667em}{0ex}}exp(-{v}_{m+\frac{1}{2}}\Delta V/D)}{1-exp(-{v}_{m+\frac{1}{2}}\Delta V/D)}\phantom{\rule{0.166667em}{0ex}},$$

(34)

where ${v}_{m+\frac{1}{2}}\left(t\right)=[{I}_{\mathrm{L}}\left({V}_{m+\frac{1}{2}}\right)+{I}_{\text{exp}}\left({V}_{m+\frac{1}{2}}\right)]/C+{\mu}_{\text{tot}}(t,{r}_{d}\left(t\right),\langle w\rangle \left(t\right))$ and $D\left(t\right)=\frac{1}{2}{\sigma}_{\text{tot}}^{2}(t,{r}_{d}\left(t\right))$ denote the drift and diffusion coefficients, respectively (cf. Eq (25)). This exponentially fitted scheme [64] is globally first order convergent [65] and yields for large drifts, |*v*_{m+½}|Δ*V* ≫ *D*, the upwind flux, sharing its stability properties. For vanishing drifts, on the other hand, the centered difference method is recovered [64], leading to more accurate solutions than the upwind scheme in regimes of strong diffusion.

For the time discretization we rewrite Eq (33) (with Eq (34)) in vectorized form and approximate the involved time derivative as first order backward difference to ensure numerical stability. This yields in each time step of length Δ*t* a linear system for the values **p**^{n+1} of the (discretized) probability density at *t*_{n+1}, given the values **p**^{n} at the previous time step *t*_{n}, and the spike rate at the time *t*_{n+1−nref} for which the refractory period has just passed,

$$(\mathbf{I}-\frac{\Delta t}{\Delta V}\phantom{\rule{0.166667em}{0ex}}{\mathbf{G}}^{n}){\mathbf{p}}^{n+1}={\mathbf{p}}^{n}+{\mathit{g}}^{n+1-{n}_{\text{ref}}},$$

(35)

with vector elements ${p}_{m}^{n}=p({V}_{m},{t}_{n})$, *m* = 1, …, *N*_{V}, and ${g}_{m}^{n+1-{n}_{\text{ref}}}={\delta}_{m{m}_{\mathrm{r}}}\frac{\Delta t}{\Delta V}r\left({t}_{n+1-{n}_{\text{ref}}}\right)$. The refractory period in time steps is given by *n*_{ref} = *T*_{ref}/Δ*t*, where the brackets denote the ceiling function, and **I** is the identity matrix. This linear equation can be efficiently solved with runtime complexity 𝒪(*N*_{V}) due to the tridiagonal structure of **G**^{n} ∈ ℝ^{NV×NV} which contains the discretization of the membrane voltage (cf. Eqs (33) and (34)), including the absorbing and reflecting boundary conditions (Eqs (31) and (32)). For details we refer to S1 Text.

The spike rate, Eq (29), in this representation is obtained by evaluating the Scharfetter-Gummel flux, Eq (34), at the spike voltage *V*_{s}, taking into account the absorbing boundary condition, Eq (31), and introducing an auxiliary ghost cell [66], with center *V*_{NV+1}, which yields

$$r\left({t}_{n+1}\right)={q}_{p}({V}_{{N}_{V}+\frac{1}{2}},{t}_{n+1})={v}_{{N}_{V}+\frac{1}{2}}\phantom{\rule{0.166667em}{0ex}}\frac{1+exp(-{v}_{{N}_{V}+\frac{1}{2}}\Delta V/D)}{1-exp(-{v}_{{N}_{V}+\frac{1}{2}}\Delta V/D)}\phantom{\rule{0.166667em}{0ex}}{p}_{{N}_{V}}^{n+1},$$

(36)

where the drift and diffusion coefficients, *v*_{NV+½} and *D*, are evaluated at *t*_{n}. The mean membrane voltage (of non-refractory neurons), Eq (28), used for the dynamics of the mean adaptation current, Eq (23), is calculated by $\langle V\rangle \left({t}_{n}\right)={\sum}_{m=1}^{{N}_{V}}{V}_{m}{p}_{m}^{n}/{\sum}_{m=1}^{{N}_{V}}{p}_{m}^{n}$.

Practically, we use the initialization ${p}_{m}^{0}={p}_{0}\left({V}_{m}\right)$ and solve in each time step the linear system, Eq (35), using the function `banded_solve` from the Python library SciPy [67]. Note that (for a recurrent network or time-varying external input) the tridiagonal matrix **G**^{n} has to be constructed in each time step *t*_{n}, which can be time consuming–especially for small Δ*V* and/or small Δ*t*. Therefore, we employ low-level virtual machine acceleration for this task through the Python package Numba [24] which yields an efficient implementation.

Remark: for a vanishing refractory period *T*_{ref} = 0 the matrix **G**^{n} would lose its tridiagonal structure due to the instantaneous reinjection, cf. Eq (36). In this case we enforce a minimal refractory period of one time step, *T*_{ref} = Δ*t*, which is an excellent approximation if the time step is chosen sufficiently small and the spike rate does not exceed biologically plausible values.

In the following sections we present two approaches of how simple spike rate models can be derived from the Fokker-Planck mean-field model described in the previous section, cf. Eqs (20), (21) and (23)–(32).

The derived models are described by low-dimensional ordinary differential equations (ODEs) which depend on a number of quantities defined in the plane of (generic) input mean and standard deviation (*μ*, *σ*). To explain this concept more clearly we consider, as an example, the steady-state spike rate, which is a quantity required by all reduced models. The steady-state spike rate as a function of *μ* and *σ*,

$${r}_{\infty}(\mu ,\sigma )\u2254\underset{t\to \infty}{lim}r(t;{\mu}_{\text{tot}}\phantom{\rule{-0.166667em}{0ex}}=\phantom{\rule{-0.166667em}{0ex}}\mu ,{\sigma}_{\text{tot}}\phantom{\rule{-0.166667em}{0ex}}=\phantom{\rule{-0.166667em}{0ex}}\sigma ),$$

(37)

denotes the stationary value of Eq (29) under replacement of the (time-varying) total moments *μ*_{tot} and ${\sigma}_{\text{tot}}^{2}$ in the probability flux *q*_{p}, Eq (25), by (constants) *μ* and *σ*^{2}, respectively. Thus the steady-state spike rate *r*_{∞} effectively corresponds to that of an uncoupled EIF population whose membrane voltage is governed by *dV*_{i}/*dt* = [*I*_{L}(*V*_{i}) + *I*_{exp}(*V*_{i})]/*C* + *μ* + *σξ*_{i}(*t*) plus reset condition, i.e., adaptation and synaptic current dynamics are detached. For a visualization of *r*_{∞}(*μ*, *σ*) see Fig 6.

When simulating the reduced models these quantities need to be evaluated for each discrete time point *t* at a certain value of (*μ*, *σ*) which depends on the overall synaptic moments *μ*_{syn}(*t*), ${\sigma}_{\text{syn}}^{2}\left(t\right)$ and on the mean adaptation current *w*(*t*) in a model-specific way (as described in the following Sects.). An example trajectory of *r*_{∞} in the (*μ*, *σ*) space for a network showing stable spike rate oscillations is shown in Fig 5.

Importantly, these quantities depend on the parameters of synaptic input (*J*, *K*, *τ*_{d}, *μ*_{ext}, *σ*_{ext}) and adaptation current (*a*, *b*, *τ*_{w}, *E*_{w}) only through their arguments (*μ*, *σ*). Therefore, for given parameter values of the EIF model (*C*, *g*_{L}, *E*_{L}, Δ_{T}, *V*_{T}, *V*_{r}, *T*_{ref}) we precalculate those quantities on a (reasonably large and sufficiently dense) grid of *μ* and *σ* values, and access them during time integration by interpolating the quantity values stored in a table. This greatly reduces the computational complexity and enables rapid numerical simulations.

The derived low-dimensional models describe the spike rate dynamics and generally do not express the evolution of the entire membrane voltage distribution. Therefore, the mean adaptation dynamics, which depends on the density *p*(*V*, *t*) (via *V*, cf. Eq (23)) is adjusted through approximating the mean membrane voltage *V* by the expectation over the steady-state distribution,

$${\u27e8V\u27e9}_{\infty}=\frac{{\int}_{-\infty}^{{V}_{\mathrm{s}}}v{p}_{\infty}\left(v\right)dv}{{\int}_{-\infty}^{{V}_{\mathrm{s}}}{p}_{\infty}\left(v\right)dv},$$

(38)

which is valid for sufficiently slow adaptation current dynamics [48, 58]. The steady-state distribution is defined as *p*_{∞}(*V*) = lim_{t → ∞}
*p*(*V*, *t*; *μ*_{tot} = *μ*, *σ*_{tot} = *σ*), representing the stationary membrane voltages of an uncoupled EIF population for generic input mean *μ* and standard deviation *σ*. The mean adaptation current in all reduced models is thus governed by

$$\frac{d\u27e8w\u27e9}{dt}=\frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{{\tau}_{w}}+b\phantom{\rule{0.166667em}{0ex}}r\left(t\right),$$

(39)

where the evaluation of quantity *V*_{∞} in terms of particular values for *μ* and *σ* at a given time *t* is model-specific (cf. following Sects.). Note again that the calculation of *V*_{∞} slightly changes when considering an (optional) spike shape extension for the aEIF model, as described at the end of the Methods section.

The Fokker-Planck model does not restrict the form of the delay distribution *p*_{d}, except that the convolution with the spike rate *r*, Eq (20), has to be well defined. Here, however, we aim at specifying the complete network dynamics in terms of a low-dimensional ODE system. Exploiting the exponential form of the delay distribution *p*_{d} we obtain a simple ordinary differential equation for the delayed spike rate,

$$\frac{d{r}_{d}}{dt}=\frac{r-{r}_{d}}{{\tau}_{d}},$$

(40)

which is equivalent to the convolution *r*_{d} = *r* * *p*_{d}.

Note that more generally any delay distribution from the exponential family allows to represent the delayed spike rate *r*_{d} by an equivalent ODE instead of a convolution integral [68]. Identical delays, *r*_{d}(*t*) = *r*(*t* − *d*), are also possible but lead to delay differential equations. Naturally, in case of no delays, we simply have *r*_{d}(*t*) = *r*(*t*).

To simulate the reduced models standard explicit time discretization schemes can be applied–directly to the first order equations of the LN_{exp} model, and for the other models (LN_{dos}, spec_{1}, spec_{2})–to the respective equivalent (real) first order systems. We would like to note that when using the explicit Euler method to integrate any of the latter three low-dimensional models a sufficiently small integration time step Δ*t* is required to prevent oscillatory artifacts. Although the explicit Euler method works well for the parameter values used in this contribution, we have additionally implemented the method of Heun, i.e., the explicit trapezoidal rule, which is second order accurate.

Following and extending [18] we can specify the Fokker-Planck operator ℒ

$$\mathcal{L}(\mu ,\sigma )\left[p\right]=-\frac{\partial}{\partial V}\phantom{\rule{1pt}{0ex}}\left[\right(\frac{{I}_{\mathrm{L}}\left(V\right)+{I}_{\text{exp}}\left(V\right)}{C}+\mu \left)\phantom{\rule{1pt}{0ex}}p\right]+\frac{{\sigma}^{2}}{2}\frac{{\partial}^{2}p}{\partial {V}^{2}},$$

(41)

for an uncoupled EIF population receiving (constant) input (*μ*, *σ*), cf. Sect. *Low-dimensional approximations*. This operator allows to rearrange the FP dynamics of the recurrent aEIF network, Eq (24), as

$$\frac{\partial p}{\partial t}=\mathcal{L}({\mu}_{\text{tot}}\left(t\right),{\sigma}_{\text{tot}}\left(t\right))\left[p\right]$$

(42)

which depends on the (time-varying) total input moments *μ*_{tot}(*t*, *r*_{d}, *w*) and ${\sigma}_{\text{tot}}^{2}(t,{r}_{d})$, cf. Eqs (26) and (27), in the drift and diffusion coefficients, respectively.

For each value of (*μ*, *σ*) the operator ℒ possesses an infinite, discrete set of eigenvalues λ_{n} in the left complex half-plane including zero [18], i.e., Re{λ_{n}} ≤ 0, and associated eigenfunctions *ϕ*_{n}(*V*) (*n* = 0, 1, 2, …) satisfying

ℒ[*ϕ*_{n}] = λ_{n}*ϕ*_{n}.

(43)

Furthermore, the boundary conditions, Eqs (30)–(32) have to be fulfilled for each eigenfunction *ϕ*_{n} separately, i.e., the absorbing boundary at the spike voltage,

(44)

and the reflecting barrier at the (finite) lower bound voltage,

(45)

must hold. The eigenflux is given by

$${q}_{{\varphi}_{n}}\left(V\right)=(\frac{{I}_{\mathrm{L}}\left(V\right)+{I}_{\text{exp}}\left(V\right)}{C}+\mu )\phantom{\rule{1pt}{0ex}}{\varphi}_{n}\left(V\right)-\frac{{\sigma}^{2}}{2}\frac{\partial}{\partial V}{\varphi}_{n}\left(V\right),$$

(46)

i.e., the flux *q*_{p} of Eq (25) with eigenfunction *ϕ*_{n} and (constant) generic input moments (*μ*, *σ*^{2}) instead of density *p* and (time-varying) total input moments, respectively. Moreover, the eigenflux *q*_{ϕn} has to be reinjected into the reset voltage, cf. Eq (30),

$${q}_{{\varphi}_{n}}\left({V}_{\mathrm{r}}^{+}\right)-{q}_{{\varphi}_{n}}\left({V}_{\mathrm{r}}^{-}\right)={q}_{{\varphi}_{n}}\left({V}_{\mathrm{s}}\right),$$

(47)

where we have neglected the refractory period, i.e., *T*_{ref} = 0. Note that incorporating a refractory period *T*_{ref} > 0 is straightforward only for the simplified case of vanishing total input moment variations, ${\dot{\mu}}_{\text{tot}}\approx {\dot{\sigma}}_{\text{tot}}^{2}\approx 0$, which is described in the following section and is not captured here in general.

The spectrum of ℒ is shown in Fig 7A and further discussed in Sect. *Remarks on the spectrum*.

Defining the non-conjugated [69] inner product $\langle \psi ,\varphi \rangle ={\int}_{{V}_{\text{lb}}}^{{V}_{\mathrm{s}}}\psi \left(v\right)\varphi \left(v\right)dv$ yields the corresponding adjoint operator ℒ^{*}(*μ*, *σ*) given by [18]

$${\mathcal{L}}^{*}=(\frac{{I}_{\mathrm{L}}\left(V\right)+{I}_{\text{exp}}\left(V\right)}{C}+\mu )\phantom{\rule{1pt}{0ex}}\frac{\partial}{\partial V}+\frac{{\sigma}^{2}}{2}\frac{{\partial}^{2}}{\partial {V}^{2}},$$

(48)

which satisfies 〈*ψ*, ℒ*ϕ*〉 = 〈ℒ^{*}*ψ*, *ϕ*〉 for any complex-valued functions *ψ* and *ϕ* that are sufficiently smooth on [*V*_{lb}, *V*_{s}]. ℒ^{*} has the same set of eigenvalues λ_{n} as ℒ but distinct associated eigenfunctions *ψ*_{n}(*V*), i.e.,

ℒ^{*}[*ψ*_{n}] = λ_{n}*ψ*_{n},

(49)

which have to satisfy three boundary conditions,

(50)

$$\frac{\partial {\psi}_{n}}{\partial V}\left({V}_{\text{lb}}\right)=0,$$

(51)

$$\frac{\partial {\psi}_{n}}{\partial V}\left({V}_{\mathrm{r}}^{+}\right)=\frac{\partial {\psi}_{n}}{\partial V}\left({V}_{\mathrm{r}}^{-}\right)$$

(52)

that are determined by integrating 〈*ψ*_{n}, ℒ[*ϕ*_{n}]〉 by parts and equalling with 〈ℒ^{*}[*ψ*_{n}], *ϕ*_{n}〉 using the conditions of ℒ, Eqs (44)–(47). Note that the last condition, Eq (52), ensures a continuous derivative of *ψ*_{n} at *V*_{r} in contrast to the eigenfunctions *ϕ*_{n} of ℒ, that have a kink at the reset due to reinjection condition, Eq (47), as shown in Fig 7B.

The eigenfunctions of ℒ and ℒ^{*} are pairwise orthogonal and in the following (without loss of generality) assumed to be scaled according to the biorthonormality condition,

⟨*ψ*_{n}, *ϕ*_{m}⟩ = *δ*_{nm}.

(53)

The membrane voltage probability density can now be expanded onto the (moving) eigenbasis of ℒ [18, 70],

$$p(V,t)=\sum _{n=0}^{\infty}{\alpha}_{n}\left(t\right){\varphi}_{n}\left(V\right),$$

(54)

where each eigenfunction *ϕ*_{n} depends on time via the total input moments *μ* = *μ*_{tot}(*t*, *r*_{d}, *w*), ${\sigma}^{2}={\sigma}_{\text{tot}}^{2}(t,{r}_{d})$, and the projection coefficients are given by *α*_{n} = *ψ*_{n}, *p* and particularly *α*_{0} = 1 [18]. Deriving *α*_{n} with respect to time (for *n* = 1, 2, …) and using the expansion, Eq (54), the Fokker-Planck Eq (42) as well as the definition of the adjoint operator yields an infinite-dimensional equation for the complex-valued projection coefficients * α*(

$$\dot{\mathit{\alpha}}=(\mathbf{\Lambda}+{\mathbf{C}}_{\mu}\phantom{\rule{0.166667em}{0ex}}{\dot{\mu}}_{\text{tot}}+{\mathbf{C}}_{{\sigma}^{2}}\phantom{\rule{0.166667em}{0ex}}{\dot{\sigma}}_{\text{tot}}^{2})\phantom{\rule{1pt}{0ex}}\mathit{\alpha}+{\mathbf{c}}_{\mu}\phantom{\rule{0.166667em}{0ex}}{\dot{\mu}}_{\text{tot}}+{\mathbf{c}}_{{\sigma}^{2}}\phantom{\rule{0.166667em}{0ex}}{\dot{\sigma}}_{\text{tot}}^{2}.$$

(55)

This dynamics is initialized by *α*_{n}(0) = *ψ*_{n}, *p*_{0}, and is complemented by (i) an expression for the spike rate,

(56)

that is obtained from Eq (29) (using Eq (54)), and (ii) by the mean adaptation and delayed spike rate dynamics, Eqs (39) and (40). The dots in Eqs (55) and (56) denote time derivatives (e.g., ${\dot{\sigma}}_{\text{tot}}^{2}\phantom{\rule{7.77771pt}{0ex}}=d{\sigma}_{\text{tot}}^{2}/dt$) and a non-conjugated scalar product of complex vectors (i.e., $\mathbf{f}\xb7\mathit{\alpha}={\sum}_{n=1}^{\infty}{f}_{n}{\alpha}_{n}$), respectively. The matrix **Λ** = diag(λ_{1}, λ_{2}, …) contains the eigenvalues of ℒ, the matrices **C**_{μ} and **C**_{σ2} have elements (**C**_{x})_{n,m} = _{x}*ψ*_{n}, *ϕ*_{m} for *x* {*μ*, *σ*^{2}} and *n*, *m* ∈ ℕ with partial derivative _{x} = /*x*, the vectors **c**_{μ} and **c**_{σ2} consist of components ${c}_{n}^{x}=\langle {\partial}_{x}{\psi}_{n},\phantom{\rule{0.166667em}{0ex}}{\varphi}_{0}\rangle $. The steady-state spike rate is given by *r*_{∞} = *q*_{ϕ0}(*V*_{s}), i.e., the flux of the eigenfunction *ϕ*_{0} = *p*_{∞} that represents the stationary membrane voltage distribution *p*_{∞} and corresponds to the (stationary) eigenvalue λ_{0} = 0 [18]. The vector **f** contains the (nonstationary) eigenfluxes evaluated at the spike voltage, *f*_{n} = *q*_{ϕn}(*V*_{s}).

Note that the quantities (**Λ**, **C**_{μ}, **C**_{σ2}, **c**_{μ}, **c**_{σ2}, *r*_{∞}, **f**, *V*_{∞}) all depend on time in Eqs (55), (56) and (39) via the total input moments (*μ*_{tot}, ${\sigma}_{\text{tot}}^{2}$). Particularly, at time *t* the (biorthonormal) solution of the eigenvalue problems for ℒ and its adjoint ℒ^{*}, Eqs (43)–(47) and (49)–(52) with *n* ∈ ℕ_{0}, is required for *μ* = *μ*_{tot}(*t*, *r*_{d}, *w*), *σ* = *σ*_{tot}(*t*, *r*_{d}). Because ℒ is a real operator its spectrum contains only real eigenvalues and/or complex conjugated pairs depending on (*μ*, *σ*). This property carries over to the eigenfunctions and therefore also to the components of all quantities above implying for example that the scalar product **f** · * α* is always real-valued.

Although the spectral representation, Eqs (55) and (56), is fully equivalent to the original (partial differential) Fokker-Planck equation without refractory period and while it contains only time derivatives it is still an infinite-dimensional and furthermore generally an implicit model. Therefore, to derive an explicit low-dimensional ordinary differential model for the spike rate *r*(*t*), it is not sufficient to truncate the expansion in Eq (54) after, e.g., two terms but additional assumptions have to be considered.

The first and simpler, derived spectral model is based on [19] and requires the strong assumption of vanishing changes of the total input moments, ${\dot{\mu}}_{\text{tot}}\approx 0$, ${\dot{\sigma}}_{\text{tot}}^{2}\phantom{\rule{7.77771pt}{0ex}}\approx 0$. Under this approximation the projection coefficient dynamics, Eq (55), simplifies to

$${\dot{\alpha}}_{n}={\lambda}_{n}{\alpha}_{n}$$

(57)

for *n* ∈ ℕ. Considering only the dominant (nonzero) eigenvalue

λ_{1} = argmin{|Re{λ_{n}}|:λ_{n} ≠ 0},

(58)

i.e., the “slowest mode”, we obtain from Eqs (56) and (57) that *r*(*t*) = *r*_{∞} + *m*_{1} Re{*f*_{1}*α*_{1}} with *m*_{1} = 1 if λ_{1} ∈ ℝ and *m*_{1} = 2 if λ_{1} ∈ ℂ∖ℝ. Here we have included for a complex eigenvalue λ_{1} also its complex conjugate ${\overline{\lambda}}_{1}$ that has the projection coefficient $\overline{{\alpha}_{1}}\left(t\right)$. They jointly yield a zero imaginary part in the scalar product of Eq (56). Defining $\tilde{r}\left(t\right)={r}_{\infty}+{m}_{1}{f}_{1}{\alpha}_{1}$ yields the complex first order equation for the spike rate [19],

$$\dot{\tilde{r}}={\lambda}_{1}(\tilde{r}-{r}_{\infty}),\phantom{\rule{2.em}{0ex}}r\left(t\right)=\text{Re}\left\{\tilde{r}\right\}.$$

(59)

While this derivation is based on neglecting changes of the total input moments, their time-variation is effectively reintroduced in the two quantities (dominant eigenvalue and steady-state rate), i.e., λ_{1}(*μ*_{tot}, *σ*_{tot}) and *r*_{∞}(*μ*_{tot}, *σ*_{tot}) with *μ*_{tot}(*t*) and *σ*_{tot}(*t*) according to Eqs (26) and (27). Therefore, the spike rate evolution, Eq (59), is complemented with the dynamics of mean adaptation current *w* and delayed spike rate *r*_{d}, Eqs (39) and (40), where the former involves the third (*μ*_{tot}, *σ*_{tot})-dependent quantity *V*_{∞}. See Figs Figs66 and and7C7C for the involved quantities depending on (generic) input *μ*, *σ*.

We call this first derived low-dimensional spike rate model, i.e., Eqs (59), (39) and (40), spec_{1}. It is very simple in comparison with the full Fokker-Planck system in the spectral representation, Eqs (55)–(56), (39) and (40), in the sense that it does not depend on “nonstationary” quantities of ℒ or its adjoint ℒ^{*} except for the dominant eigenvalue λ_{1}.

Note that under the assumption of vanishing input moment variations the dynamics of the expansion coefficients *α*_{n}, Eq (55), is simply exponentially decaying in time, cf. Eq (57). This allows to incorporate a refractory period *T*_{ref} > 0 into the spectral decomposition framework by inserting the eigenbasis expansion of the membrane voltage distribution, Eq (54), into the reinjection condition of the Fokker-Planck model, Eq (30), using *α*_{n}(*t*) = *α*_{n}(0) exp(λ_{n}*t*) and the absorbing boundary, Eq (44). This generalizes the reinjection condition for the eigenfunctions *ϕ*_{n} of ℒ from Eq (47) to

$${q}_{{\varphi}_{n}}\left({V}_{\mathrm{r}}^{+}\right)-{q}_{{\varphi}_{n}}\left({V}_{\mathrm{r}}^{-}\right)={q}_{{\varphi}_{n}}\left({V}_{\mathrm{s}}\right)exp(-{\lambda}_{n}{T}_{\text{ref}}),$$

(60)

which was applied in [19], and the corresponding boundary condition for the adjoint operator’s eigenfunctions *ψ*_{n} from Eq (50) to

(61)

Another possibility to derive a low-dimensional spike rate model is based on ongoing work of Maurizio Mattia. He has recognized that under the weaker (compared to the basic spec_{1} model above) assumption of small but not vanishing input moment changes a real-valued second order ordinary differential equation for the spike rate *r*(*t*) can be consistently derived [20]. Here we extend this approach to account for neuronal adaptation, time-varying external input moments and delay distributions.

The steps in a nutshell (for the detailed derivation see S1 Text) are (i) taking the derivative of Eq (55) (once) and of (56) (twice) w.r.t. time, (ii) considering only the first two dominant eigenvalues λ_{1} and λ_{2}, i.e., neglecting all (faster) eigenmodes that correspond to eigenvalues with larger absolute real part (“modal approximation”), (iii) assuming slowly changing input moments, i.e., small ${\dot{\mu}}_{\text{tot}}$ and ${\dot{\sigma}}_{\text{tot}}^{2}$ that allow to consider projections coefficients of that order, i.e., ${\alpha}_{n}=\mathcal{O}\left({\dot{\mu}}_{\text{tot}}\right)$ and ${\alpha}_{n}=\mathcal{O}\left({\dot{\sigma}}_{\text{tot}}^{2}\right)$ for *n* = 1, 2, and therefore to consider only linear occurences of ${\dot{\mu}}_{\text{tot}}$, ${\dot{\sigma}}_{\text{tot}}^{2}$, *α*_{1}, *α*_{2}, and neglect terms of second and higher order. The slowness approximation implies that neither the external moments *μ*_{ext}(*t*), ${\sigma}_{\text{ext}}^{2}\left(t\right)$ nor the (delayed) spike rate *r*_{d}(*t*) nor the population-averaged adaptation current *w*(*t*) should change very fast (cf. Eqs (26), (27) and (21). Note that the dynamics of *w* is assumed to be slow already, cf. Sect. *Fokker-Planck system*.

Under these approximations we obtain for the spike rate dynamics the following real second order ODE,

$${\beta}_{2}\phantom{\rule{0.166667em}{0ex}}\ddot{r}+{\beta}_{1}\phantom{\rule{0.166667em}{0ex}}\dot{r}+{\beta}_{0}\phantom{\rule{0.166667em}{0ex}}r={r}_{\infty}-r-{\beta}_{c},$$

(62)

that is complemented with the mean adaptation current and delayed spike rate dynamics, Eqs (39) and (40), and we call the model spec_{2}. The coefficients

(63)

$${\beta}_{1}=\phantom{\rule{0.166667em}{0ex}}-T+DM\frac{b}{C}-\frac{R}{{\tau}_{d}},$$

(64)

$${\beta}_{0}=-DM\frac{b}{{\tau}_{w}C}-\frac{b}{C}{H}_{\mu}+\frac{1}{{\tau}_{d}}\phantom{\rule{1pt}{0ex}}(KJ{H}_{\mu}+K{J}^{2}{H}_{{\sigma}^{2}})+\frac{R}{{\tau}_{d}^{2}},$$

(65)

$$\begin{array}{ccc}\hfill {\beta}_{\text{c}}& \hfill =\hfill & \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}{r}_{d}\phantom{\rule{1pt}{0ex}}(-\frac{1}{{\tau}_{d}}(KJ{H}_{\mu}+\phantom{\rule{0.166667em}{0ex}}K{J}^{2}{H}_{{\sigma}^{2}})-\frac{R}{{\tau}_{d}^{2}})\hfill \\ & -& ({\ddot{\mu}}_{\text{ext}}+\frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{{\tau}_{w}^{2}C})\phantom{\rule{1pt}{0ex}}DM-\ddot{{\sigma}_{\text{ext}}^{2}\phantom{\rule{-7.77771pt}{0ex}}}\phantom{\rule{7.77771pt}{0ex}}DS\hfill \\ & +& ({\dot{\mu}}_{\text{ext}}-\frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{{\tau}_{w}C})\phantom{\rule{1pt}{0ex}}{H}_{\mu}+\dot{{\sigma}_{\text{ext}}^{2}\phantom{\rule{-7.77771pt}{0ex}}}\phantom{\rule{7.77771pt}{0ex}}{H}_{{\sigma}^{2}},\hfill \end{array}$$

(66)

depend on the (lumped) quantities *D* = 1/λ_{1} 1/λ_{2}, *T* = 1/λ_{1} + 1/λ_{2}, *M* = _{μ}*r*_{∞} + **f** **c**_{μ}, *S* = _{σ2}*r*_{∞} + **f** **c**_{σ2}, *R* = *DMKJ* + *DSKJ*^{2}, *H*_{μ} = *TM* + *DMa*_{μ}*V*_{∞}/(*τ*_{w}*C*) − *DF*_{μ}, *H*_{σ2} = *TS* + *DMa*_{σ2}*V*_{∞}/(*τ*_{w}*C*) − *DF*_{σ2}, *F*_{μ} = **f** **Λ**
**c**_{μ} and *F*_{σ2} = **f** **Λ**
**c**_{σ2}. Here the diagonal eigenvalue matrix **Λ** = diag(λ_{1}, λ_{2}) and the vectors **f** = (*f*_{1}, *f*_{2})^{T}, ${\mathbf{c}}_{\mu}={({c}_{1}^{\mu},{c}_{2}^{\mu})}^{T}$, ${\mathbf{c}}_{{\sigma}^{2}}={({c}_{1}^{{\sigma}^{2}},{c}_{2}^{{\sigma}^{2}})}^{T}$ are two-dimensional in contrast to infinite as in the original dynamics, Eqs (55) and (56). These individual quantities, that also include *V*_{∞} (cf. Eq (39)), *r*_{∞} and derivatives of both w.r.t. generic mean *μ* and variance *σ*^{2}, are evaluated at the total input moments *μ* = *μ*_{tot}(*t*, *r*_{d}, *w*), ${\sigma}^{2}={\sigma}_{\text{tot}}^{2}(t,{r}_{d})$. A relevant subset of individual and lumped quantities is shown in Figs Figs66 and and7C7C.

The four coefficients, Eqs (63)–(66), are real-valued because we define λ_{2} as the second dominant eigenvalue conditioned that (λ_{1}, λ_{2}) compose either a real or a complex conjugated eigenvalue pair, i.e.,

λ_{2} = argmin{|Re{λ_{n}}| : λ_{n} ≠ λ_{1} s.t. λ_{n} ≠ 0, λ_{1} + λ_{n} ∈ ℝ},

(67)

where λ_{1} is obtained as for the (basic) spec_{1} model, cf. Eq (58). This condition ensures that all related complex quantities occur in vectors of two complex conjugate components (for example, ${f}_{1}=\overline{{f}_{2}}$) implying that the scalar products above are real-valued (e.g., **f** · **c**_{μ} ∈ ℝ) and therefore all (nine) lumped quantities, too. Note that this specific definition of λ_{2} is required only for integrate-and-fire neuron models that have a lower bound different from the reset, i.e., *V*_{lb} < *V*_{r}, as discussed in the following section.

The coefficients of the spec_{2} model require–in addition to eigenvalues λ_{n} and steady-state rate *r*_{∞} (cf. basic spec_{1} model, Eq (59))–quantities that involve the first and second eigenfunctions of the Fokker-Planck operator ℒ and its adjoint ℒ^{*}. Additionally, *β*_{1}, *β*_{0} and *β*_{c} contain the parameters of membrane voltage (*C*) and mean adaptation current dynamics (*a*, *b*, *τ*_{w}, *E*_{w}) as well as of the recurrent coupling (*K*, *J*, *τ*_{d}) and, importantly, explicit dependencies on the population-averaged adaptation current *w* and the delayed spike rate *r*_{d} (that is in addition to implicitly via *μ*_{tot} and *σ*_{tot}). Furthermore *β*_{c} depends on the first two time derivatives of the external input moments *μ*_{ext}(*t*) and ${\sigma}_{\text{ext}}^{2}\left(t\right)$. This explicit occurence of neuronal and coupling parameters, state variables and input moment derivatives in the coefficients is not expressed in the basic spectral model (spec_{1}), Eq (59). A consequence is that for the (advanced) model spec_{2} the external moments have to be provided twice differentiable or in case of non-smooth time series to be filtered (e.g., see Sect. *Performance for variations of the mean input*).

The particular coefficients, Eqs (63)–(66), are specific for the choice of an exponential delay distribution, as indicated by the occurrence of the mean delay *τ*_{d} in *β*_{1}, *β*_{0} and *β*_{c}. Other choices such as identical delays or no delays are described in the supplementary material S1 Text.

Note that the original (infinite-dimensional) spectral dynamics, Eqs (55) and (56), assumes for the refractory period a value of *T*_{ref} = 0 which carries over to the same choice for the spec_{2} model, whereas *T*_{ref} > 0 is not supported (yet).

The spike rate is by definition nonnegative, however, the model spec_{2}, Eq (62), can yield negative rates *r*(*t*), especially for small total mean input *μ*_{tot} and fast (external) input changes, e.g., when ${\dot{\mu}}_{\text{ext}}$ is large. We explicitly avoid that behaviour by setting both the spike rate of this model, *r*, and its derivative, $\dot{r}$, to zero whenever *r*(*t*)<0 and continue the integration of the differential equation afterwards.

In the previous two sections we have developed two spike rate models based on approximations of the Fokker-Planck system’s spectral representation, Eqs (55) and (56) under different slowness assumptions. In the derivation of the (simple) model spec_{1}, cf. Eq (59), temporal variations of the total input moments are completely neglected while the (advanced) model spec_{2}, cf. Eq (62), incorporates (slow) changes of the total input moments through linear terms (proportional to ${\dot{\mu}}_{\text{tot}}$ or ${\dot{\sigma}}_{\text{tot}}^{2}$) and neglects (faster) quadratic and higher order ones.

Both slowness approximations imply that the eigenvalue matrix Λ is the dominant term in the homogeneous part of Eq (55). Therefore the eigenvalues approximately correspond to decay time constants 1/|Re{λ_{n}}| in their real part and in case of complex eigenvalues they contribute damped oscillatory components with frequency |Im{λ_{n}}|/(2*π*) through their imaginary part to the dynamics. How the spectrum of the Fokker-Planck operator ℒ depends on the input moments therefore gives insights into the behavior of the derived spike rate models spec_{1} and spec_{2}, and of the FP model with slow adaptation in general.

Here we summarize the main properties of the eigenvalues λ_{n}(*μ*, *σ*) for the (uncoupled, nonadaptive) EIF neuron model as a function of generic input mean *μ* and standard deviation *σ* (cf. Sect. *Low-dimensional approximations*) that are shown in Fig 7A. Note that at time *t* the total input determines the particularly effective eigenvalue, λ_{n}(*μ*_{tot}(*t*), *σ*_{tot}(*t*)), in dependence of external input moments, (delayed) spike rate and adaptation current, cf. Eqs (26), (27) and (21. We thereby extend the findings of [18] concerning the perfect integrate-and-fire neuron with reflecting lower barrier at the reset voltage (PIFb), i.e., *dV*_{i}/*dt* = *μ* + *σξ*_{i}(*t*) with *V*_{lb} = *V*_{r}.

- The eigenvalue λ
_{0}= 0 exists for all (generic) input moments*μ*,*σ*^{2}with stationary membrane voltage distribution as corresponding eigenfunction,*ϕ*_{0}=*p*_{∞}(see Fig 7A and 7B), the respective adjoint eigenfunction is constant,*ψ*_{0}= 1, cf. [18]. The other, nonstationary eigenvalues λ_{n}with*n*≥ 1 have negative real parts for all (*μ*,*σ*) which yields stability of the stationary distribution for constant total input moments. For sufficiently small mean input*μ*_{min}the eigenvalues λ_{n}(*μ*_{min},*σ*) are real-valued for all*n*and*σ*, i.e., no (damped) oscillatory spike rate transient in that regime are present (consistent with [18]), cf. Fig 7A. - Two classes of modes can be distinguished: eigenvalues of the first kind occur in couples which are real-valued at
*μ*_{min}but merge for increasing mean input*μ*to become a complex conjugated pair of decreasing absolute real part and almost linearly increasing imaginary part (see Fig 7A). Thus they correspond in this situation to damped oscillatory dynamics with increased frequency |Im{λ_{n}}|/(2*π*) and decay time constant 1/|Re{λ_{n}}| for stronger mean input. Note that there is no single critical parameter for the real-to-complex transition, instead that depends on the input mean*μ*and standard deviation*σ*as well as on the eigenvalue index*n*in contrast to the PIFb neuron where*μ*= 0 (alone) induces the transition [18]. We call this type of eigenvalue*regular*because it is observed also for very simple integrate-and-fire models (such as PIFb).The second type of eigenvalue is real for the whole input parameter space of (*μ*,*σ*). The corresponding decay time constant is large if noise dominates while for stronger mean input the respective dynamics is negligibly fast. Note that by setting the lower bound*V*_{lb}equal to the reset voltage, this eigenvalue class is completely removed for the EIFb model, which is the EIF membrane voltage description with*V*_{lb}=*V*_{r}(not used here). This explains why this new type of eigenvalue has not been described in [18], and suggests a relationship to the diffusion of the membrane voltage for hyperpolarized neuronal states. The latter correspondence is further supported from the significant values of the respective eigenfunction*ϕ*_{n}below the reset voltage*V*_{r}in contrast to the eigenfunctions of regular eigenvalues (cf. Fig 7B). Therefore we label this second type of eigenvalue*diffusive*. - The input noise intensity
*σ*controls the spectrum’s mixture of the two eigenvalue classes as follows: weak noise favors regular modes, i.e., the dominant two eigenvalues are pairs of real (for smaller mean input*μ*) or complex conjugated eigenvalues (for larger*μ*) while the diffusive modes are irrelevantly fast in this regime (cf. Fig 7A, left). Increased input fluctuations, i.e., a larger*σ*, on the other hand leads to a spectrum with dominant (“slowest” eigenvalue λ_{1}) of the diffusive kind for small mean input*μ*while for larger*μ*the dominant mode is from the regular class (being real or complex depending on*μ*,*σ*,*n*), see Fig 7A, right. Furthermore an increased noise strength*σ*leads to a smaller decay time constant 1/|Re{λ_{n}}| for regular modes, i.e., their respective contribution to the spike rate dynamics is faster in the fluctuation-dominated regime than in the drift-dominated one, whereas for diffusive modes the opposite holds true.

The specific definition of the second dominant eigenvalue λ_{2}, Eq (67), is necessary to ensure real coefficients in the spec_{2} model, Eqs (63)–(66), for regions in (*μ*, *σ*)-space where the first dominant eigenvalue λ_{1} is of the (real) diffusive type and the second is part of a complex conjugate couple (for an example see Fig 7A, right).

Extracting the required dominant eigenvalues λ_{1} (for both spectral models: spec_{1}, spec_{2}) and λ_{2} (for spec_{2} model only) according to Eqs (58) and (67) in the (*μ*, *σ*)-plane leads to points of instantaneous changes in the imaginary part (λ_{1} and λ_{2}) and the real part (only λ_{2}) due to transitions from a dominant diffusive mode for small mean input *μ* to a regular eigenvalue (or pair) becoming dominant for larger *μ* (see Fig 7A and 7C and the last property above). These discontinuities (which lie on a one-dimensional curve *μ*(*σ*)) could be avoided by either restricting to *V*_{lb} = *V*_{r} (no diffusive modes) or by deriving a third spectral model based on three eigenvalues: the dominant regular pair together with the dominant diffusive eigenvalue. Making the latter extension is straightforward by using the same steps and slowness approximation as for the model spec_{2} and would yield a 3rd order ODE for the spike rate *r*(*t*) with smooth coefficients and is expected to show increased reproduction accuracy (especially for small mean input) compared to the model spec_{2}.

The properties above imply for the low-dimensional models spec_{1} and spec_{2}, that both enable damped oscillatory spike rates for mean-dominated input (for large *μ*) since then the first two dominant eigenvalues (λ_{1} and λ_{2}) have nonzero imaginary parts and are complex conjugates of each other (see Fig 7C). Furthermore in this case the effective time constant, 1/|Re{λ_{1}}| = 1/|Re{λ_{2}}|, is large (especially for small input standard deviation *σ*). For noise-dominated input, i.e., when *σ* (and not *μ*) is large, on the other hand, the corresponding spike rate dynamics is fast and does not contain an oscillatory component.

Here we present a numerical solution method of the Fokker-Planck boundary eigenvalue problem (BEVP) for the operator ℒ(*μ*, *σ*), Eqs (43)–(47), and its adjoint ℒ^{*}(*μ*, *σ*), Eqs (49)–(52). The solution of the two BEVPs in terms of eigenvalues λ_{n} and associated (biorthonormal) eigenfunctions (*ϕ*_{n}(*V*), *ψ*_{n}(*V*)) as well as quantities that are derived from those and are required by the models spec_{1} and spec_{2}, is obtained for a rectangle of the input mean *μ* and standard deviation *σ* (see Sect. *Low-dimensional Approximations*). Note that the numerical method is not restricted to the EIF neuron model and supports other integrate-and-fire models as well (e.g., perfect, leaky or quadratic).

The eigenequation ℒ(*μ*, *σ*)[*ϕ*] = λ*ϕ*, i.e., Eq (43) (omitting the eigenvalue index *n*), represents a second order ODE (cf. Eq (41)), that is equivalent to the following first order system for the eigenflux *q*_{ϕ}(*V*) and the eigenfunction *ϕ*(*V*),

$$-\frac{d}{dV}\left(\begin{array}{c}{q}_{\varphi}\\ \varphi \end{array}\right)=\underset{=\text{A}}{\underbrace{\left(\begin{array}{cc}0& \lambda \\ \frac{2}{{\sigma}^{2}}& -2\frac{g\left(V\right)+\mu}{{\sigma}^{2}}\end{array}\right)}}\phantom{\rule{1pt}{0ex}}\left(\begin{array}{c}{q}_{\varphi}\\ \varphi \end{array}\right)$$

(68)

with coefficient matrix **A**(*V*) that has a nonlinear component through *g*(*V*) = [*I*_{L}(*V*) + *I*_{exp}(*V*)]/*C* that contains leak and exponential (membrane) currents. Here it was used the form of the eigenflux *q*_{ϕ} (cf. Eq (46)), and that ℒ(*μ*, *σ*)[*ϕ*] = - ∂_{V}*q*_{ϕ} for generic input moments *μ* and *σ*^{2} (cf. Eq (41)).

Basically a direct discretization, e.g., by a finite difference approximation, of the membrane voltage derivatives in the System (68) can be applied. In combination with the boundary conds., Eqs (44)–(47), this leads to a (sparse) matrix eigenvalue problem that allows for application of standard (Arnoldi iteration based) numerical solvers. However, the convergence properties of this approach are very poor in the sense that extremely small voltage steps Δ*V* have to be chosen for the finite differences. Thus, the technique is inefficient as huge systems appear or even inaccurate due to amplified round-off errors by ill-conditioned matrices.

Here we propose an alternative solution procedure which is based on a reformulation of the System (68) with corresponding boundary conditions, Eqs (44)–(47), as a complex-valued algebraic root finding problem

$$\lambda \text{\hspace{0.17em}}\mapsto \text{\hspace{0.17em}}{q}_{\varphi}({V}_{\text{lb}};\lambda )\stackrel{!}{=}0$$

(69)

whose solutions are the eigenvalues λ_{n}. To evaluate the nonlinear function (the left hand side of this equation) for an arbitrary λ ∈ ℂ (not necessary an eigenvalue) Eq (68) is integrated backward starting from the spike voltage *V*_{s} (initializing one component to satisfy the absorbing boundary cond., Eq (44), i.e., *ϕ*(*V*_{s}) = 0, and another component that can be chosen arbitrarily, *q*_{ϕ}(*V*_{s}) ∈ ℂ∖{0}, due to the linearity of the problem) via the reset *V*_{r} (where the reinjection cond., Eq (47), is enforced, i.e., ${q}_{\varphi}\left({V}_{\mathrm{r}}^{-}\right)={q}_{\varphi}\left({V}_{\mathrm{r}}^{+}\right)-{q}_{\varphi}\left({V}_{\mathrm{s}}\right)$, which induces a discontinuity in *q*_{ϕ} at the reset voltage) finally to the lower bound voltage *V*_{lb}. There, a nonzero value of *q*_{ϕ}(*V*_{lb}) ∈ ℂ indicates that λ is not an eigenvalue since in this case *ϕ*(*V*) violates the reflecting boundary condition, Eq (45). *q*_{ϕ}(*V*_{lb}) = 0 on the other hand shows that λ is an eigenvalue with corresponding eigenfunction *ϕ*(*V*) and eigenflux *q*_{ϕ}(*V*). Note that when considering a nonzero refractory period the generalized version of the reinjection cond., Eq (60), i.e., ${q}_{\varphi}\left({V}_{\mathrm{r}}^{-}\right)={q}_{\varphi}\left({V}_{\mathrm{r}}^{+}\right)-{q}_{\varphi}\left({V}_{\mathrm{s}}\right){e}^{-\lambda {T}_{\text{ref}}}$, is enforced instead of the expression above. The latter is only valid for the spec_{1} model, Eq (59) (see Sect. *Basic model: one eigenvalue, negligible input variations*) and makes the Fokker-Planck eigenvalue problem nonlinear due to the exponentiation of λ in Eq (60).

The (complex-valued) root finding problem, Eq (69) can be solved numerically to yield a target eigenvalue λ_{n} using an iterative procedure (for example a variant of Newton’s method) given that a sufficiently close initial approximation ${\tilde{\lambda}}_{n}\in \mathbb{C}$ is available. In our Python implementation we apply Powell’s hybrid method as implemented in MINPACK wrapped through SciPy [67] to the equivalent real system of the two variables Re{λ} and Im{λ}.

Appropriate initial approximations ${\tilde{\lambda}}_{n}$ can be achieved (i) by exploiting that for sufficiently small generic mean input *μ*_{min} all eigenvalues have zero imaginary part (see Sect. *Remarks on the spectrum*). In that case the eigenvalues are given by the roots {λ_{0}, λ_{1}, …} of *q*(*V*_{lb}; λ)–the one-dimensional function of the real-valued eigenvalue candidate λ (−∞, 0]–which are obtained, for example, by (dense) evaluation of that function in a sufficiently large (sub)interval below zero. Furthermore (ii) all eigenvalues depend continuously on the (input) parameters *μ* and *σ* (cf. Fig 7A), i.e., for a small step in the respective parameter space the solution λ_{n} of Eq (69) for the last parametrization is a very good initial value ${\tilde{\lambda}}_{n}$ for the current parametrization. With this initialization the solution of Eq (69) is typically found in a few steps of the chosen Newton-like method (see S1 Text for details).

To efficiently and accurately evaluate *q*_{ϕ}(*V*_{lb}; λ) in each iteration (of the root finding algorithm), we perform an exponential (backward) integration of the System (68). The resulting scheme is based on truncating the Magnus expansion of the exact solution after one term [71]. Here matrix exponential function evaluations of the form exp[**A**(*V*)Δ*V*] occur that are calculated using an analytic expression [72]. Note that ODE solvers that either do not exploit the linear structure of Eq (68) at all (such as the explicit Euler method but also higher order Runge-Kutta methods) or utilize linearity only in one variable (e.g., [73]) have poor convergence behaviour and thus require very small step sizes Δ*V* due to the strong nonlinearity *g*(*V*) as a consequence of the large value of the current *I*_{exp} close to *V*_{s} and significantly large absolute eigenvalues |λ_{n}|, respectively.

For more information on the numerical solver we refer to S1 Text, where details regarding the adjoint operator, the exponential integration, the initial eigenvalue determination and the (parallelized) treatment of the input parameters are included.

Linear-Nonlinear (LN) cascade models of neuronal activity are often applied in neuroscience, because they are simple and efficient, and the model components can be estimated using established experimental procedures [21, 74, 75]. Here we use the LN cascade as an ansatz to develop a low-dimensional model and we determine its components from the underlying Fokker-Planck model. This section builds upon [21] and extends that approach for recurrently coupled aEIF neurons; specifically, by taking into account an adaptation current and variations of the input variance. Furthermore, we consider an improved approximation of the derived linear filters and include an (optional) explicit description of the spike shape, cf. [23] (ch 4.2).

The cascade models considered here produce spike rate output by applying to the time-varying mean *μ*_{syn} and standard deviation *σ*_{syn} of the (overall) synaptic input, cf. Eq (21), separately a linear temporal filter, *D*_{μ} and *D*_{σ}, followed by a common nonlinear function *F*. That is,

(70)

(71)

(72)

where *μ*_{f} and *σ*_{f} denote the filtered mean and filtered standard deviation of the input, respectively. ${D}_{\mu}*{\mu}_{\text{syn}}\left(t\right)={\int}_{0}^{\infty}{D}_{\mu}\left(\tau \right){\mu}_{\text{syn}}(t-\tau )d\tau $ is the convolution between *D*_{μ} and *μ*_{syn}. The filters *D*_{μ}, *D*_{σ} are adaptive in the sense that they depend on the mean adaptation current *w* and on the (arbitrary) baseline input in terms of baseline mean ${\mu}_{\text{syn}}^{0}$ and standard deviation ${\sigma}_{\text{syn}}^{0}$. For improved readability these dependencies are not explicitly indicated in Eqs (71) and (72). Note, that the nonlinearity *F* also depends on *w*, which is governed by Eq (23). Since the mean adaptation current depends on the mean membrane voltage *V* we also consider a nonlinear mapping *H* for that population output quantity,

⟨*V*⟩(*t*) = *H* (*μ*_{f}, *σ*_{f}, ⟨*w*⟩).

(73)

For the derivation below it is instructive to first consider an uncoupled population, i.e., the input moments do not depend on *r*_{d} for now. In particular, the input statistics are described by ${\mu}_{\text{syn}}\left(t\right)={\mu}_{\text{syn}}^{0}+{\mu}_{\text{syn}}^{1}\left(t\right)$ and ${\sigma}_{\text{syn}}\left(t\right)={\sigma}_{\text{syn}}^{0}+{\sigma}_{\text{syn}}^{1}\left(t\right)$. In the following, we derive the components *F*, *D*_{μ} and *D*_{σ} from the Fokker-Planck model for small amplitude variations ${\mu}_{\text{syn}}^{1}$, ${\sigma}_{\text{syn}}^{1}$ and for a slowly varying adaptation current (as already assumed). We then approximate the derived linear filter components using suitable functions such that the convolutions can be expressed in terms of simple ODEs. Finally, we account for time-varying baseline input (${\mu}_{\text{syn}}^{0}\left(t\right)$, ${\sigma}_{\text{syn}}^{0}\left(t\right)$) and for recurrent coupling in the resulting low-dimensional spike rate models.

We first expand *F* in Eq (70) around the baseline ${\mu}_{\mathrm{f}}={\mu}_{\text{syn}}^{0}$, ${\sigma}_{\mathrm{f}}={\sigma}_{\text{syn}}^{0}$, *w* = *w*_{0} to linear order, assuming that the amplitudes of ${\mu}_{\text{syn}}^{1}$ and ${\sigma}_{\text{syn}}^{1}$ are small, and the mean adaptation current varies slowly compared to the input moments, to obtain the approximation for Eq (70),

$$\begin{array}{ccc}\hfill r\left(t\right)\approx F({\mu}_{\text{syn}}^{0},{\sigma}_{\text{syn}}^{0},\u27e8w\u27e9)& \hfill +\hfill & {D}_{\mu}*{\mu}_{\text{syn}}^{1}\left(t\right)\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial \mu}F({\mu}_{\text{syn}}^{0},{\sigma}_{\text{syn}}^{0},\u27e8w\u27e9)\hfill \\ & +& {D}_{\sigma}*{\sigma}_{\text{syn}}^{1}\left(t\right)\phantom{\rule{0.166667em}{0ex}}\frac{\partial}{\partial \sigma}F({\mu}_{\text{syn}}^{0},{\sigma}_{\text{syn}}^{0},\u27e8w\u27e9).\hfill \end{array}$$

(74)

Due to slow adaptation (*w*(*t*) = *w*_{0} + *w*_{1}(*t*) with vanishing *w*_{1}) we have neglected the expansion term in the direction of *w* and replaced *w*_{0} = *w* in the approximation above. Note also that, without loss of generality, we have assumed normalized filters, ${\int}_{0}^{\infty}{D}_{\mu}\left(\tau \right)d\tau ={\int}_{0}^{\infty}{D}_{\sigma}\left(\tau \right)d\tau =1$. Under the same assumptions the output from the Fokker-Planck model (Eqs (23)–(32)) can be approximated as

$$r\left(t\right)\approx {r}_{\infty}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})+{R}_{\mu}*{\mu}_{\text{syn}}^{1}\left(t\right)+{R}_{\sigma}*{\sigma}_{\text{syn}}^{1}\left(t\right),$$

(75)

$$\frac{d\u27e8w\u27e9}{dt}\approx \frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{{\tau}_{w}}+b\phantom{\rule{0.166667em}{0ex}}r\left(t\right),$$

(76)

$${\mu}_{\text{tot}}^{0}\left(t\right)={\mu}_{\text{syn}}^{0}-\u27e8w\u27e9/C,\phantom{\rule{2.em}{0ex}}{\sigma}_{\text{tot}}^{0}={\sigma}_{\text{syn}}^{0},$$

(77)

where *r*_{∞} and *V*_{∞} are the steady-state spike rate and mean membrane voltage of a population of EIF neurons in response to an input of *total* mean ${\mu}_{\text{tot}}^{0}$ plus Gaussian white noise with standard deviation ${\sigma}_{\text{tot}}^{0}$. In particular, *V*_{∞} reflects the mean over all nonrefractory neurons, cf. Eq (38). *R*_{μ} and *R*_{σ} are the so-called linear rate response functions of the population for weak modulations of the input mean and standard deviation around ${\mu}_{\text{tot}}^{0}$ and ${\sigma}_{\text{tot}}^{0}$, respectively [21, 73, 76]. Comparing Eqs (74) and (75) we obtain for the nonlinearity *F* and the linear filters *D*_{μ}, *D*_{σ},

(78)

$${D}_{\mu}\left(t\right)=\frac{{R}_{\mu}\left(t\right)}{\frac{\partial}{\partial \mu}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})},$$

(79)

$${D}_{\sigma}\left(t\right)=\frac{{R}_{\sigma}\left(t\right)}{\frac{\partial}{\partial \sigma}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})}.$$

(80)

Furthermore, the function *H* in Eq (73) is given by

(81)

ensuring that the mean membrane voltage corresponds with the instantaneous spike rate estimate of the model (in every time step). Note that *R*_{μ} and *R*_{σ} depend on ${\mu}_{\text{tot}}^{0}$ and ${\sigma}_{\text{tot}}^{0}$ (which is again not explicitly indicated for improved readability).

Fortunately, the quantities *r*_{∞}, *V*_{∞}, *R*_{μ}, and *R*_{σ} can be calculated from the Fokker-Planck equation using an efficient numerical method [73]. In particular, for *r*_{∞} and *V*_{∞} we need to solve a linear boundary value problem (BVP), and *R*_{μ}(*t*), *R*_{σ}(*t*) are calculated in the Fourier domain, where we need to solve two linear BVPs to obtain ${\widehat{R}}_{\mu}\left(f\right)$, ${\widehat{R}}_{\sigma}\left(f\right)$ for each frequency *f*. It is worth noting that the refractory period is included in a straightforward way and does not increase the complexity of the BVPs to be solved, see [23, 73] (and our provided code).

To express the LN model with adaptation, Eqs (70)–(73), (23) and (78)–(81), in terms of a low-dimensional ODE system we next approximate the linear filters *D*_{μ} and *D*_{σ} using suitable functions. The shapes of the true filters (proportional to *R*_{μ} and *R*_{σ}, cf. Eqs (79) and (80)) for different input parameter values ($\mu ={\mu}_{\text{tot}}^{0}$, $\sigma ={\sigma}_{\text{tot}}^{0}$) are shown in Fig 8A and 8B.

We first consider the linear filter *D*_{μ} (Eq (79)) and apply the approximation

(82)

motivated by the exponential decay exhibited by *R*_{μ}, particularly for large input variance compared to its mean. Note that *D*_{μ} depends on ${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$ and therefore the scaling parameter *A*_{μ} and the time constant *τ*_{μ} both depend on ${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$, which is not explicitly indicated. *A*_{μ} and *τ*_{μ} may be determined analytically using asymptotic results for the Fourier transform ${\widehat{R}}_{\mu}$ of the linear rate response for vanishing and very large frequencies, respectively [21, 76],

$$\underset{f\to 0}{lim}{\widehat{R}}_{\mu}\left(f\right)=\frac{\partial}{\partial \mu}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0}),\phantom{\rule{2.em}{0ex}}\underset{f\to \infty}{lim}\phantom{\rule{1pt}{0ex}}{\widehat{R}}_{\mu}\left(f\right)=\frac{{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})}{i2\pi f{\Delta}_{\mathrm{T}}}.$$

(83)

To guarantee that these asymptotics are matched by the Fourier transform *A*_{μ}*τ*_{μ}/(1 + *i*2*πfτ*_{μ}) of the exponential, taking into account the scaling factor in Eq (79), we obtain

$${A}_{\mu}=\frac{1}{{\tau}_{\mu}},\phantom{\rule{2.em}{0ex}}{\tau}_{\mu}=\frac{{\Delta}_{\mathrm{T}}}{{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})}\frac{\partial}{\partial \mu}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0}).$$

(84)

Note that matching the zero frequency limit in the Fourier domain is equivalent to the natural requirement that the time integral ${\int}_{0}^{\infty}{D}_{\mu}\left(\tau \right)d\tau $ of the linear filter is reproduced exactly, that is, the approximation is normalized appropriately. An advantage of this approximation is that it is no longer required to calculate the linear rate response function *R*_{μ} explicitly. On the other hand, as only the limit *f* → ∞ is used for fitting in addition to the normalization constraint, the approximation of the linear filter can be poor for a range of intermediate frequencies (in the Fourier domain), in particular for small input mean and standard deviation (see Fig 8A here, and Fig. 4B,C in [21]). To improve the approximation for intermediate frequencies we use the same normalization condition, which fixes the parameter *A*_{μ} = 1/*τ*_{μ}, and we determine *τ*_{μ} by a least-squares fit of ${\widehat{D}}_{\mu}$ over the range of frequencies *f* [0, 1] kHz. In both cases, using the approximation Eq (82) the filtered mean input *μ*_{f}(*t*) = *D*_{μ} * *μ*_{syn}(*t*) can be equivalently obtained by solving the simple scalar ODE,

$$\frac{d{\mu}_{\mathrm{f}}}{dt}=\frac{{\mu}_{\text{syn}}-{\mu}_{\mathrm{f}}}{{\tau}_{\mu}}.$$

(85)

Recall that *τ*_{μ} depends on ${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$. This exponentially decaying filter is part of the LN_{exp} cascade model variant.

A shortcoming of the approximation Eqs (82) and (85) above is that it cannot reproduce damped oscillations exhibited by the true linear filter, in particular, for large input mean and small variance (see Fig 8A). Therefore, we introduce an alternative approximation using a damped oscillator function,

(86)

Note that here *B*_{μ}, *τ* and *ω* depend on ${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$, which is not explicitly indicated. We fix the scaling parameter *B*_{μ} by (again) requiring that the approximation is normalized to reproduce the time integral of the true linear filter *D*_{μ}, which yields *B*_{μ} = (1 + (*τ*^{2}
*ω*^{2}))/*τ*. The remaining two parameters *τ* and *ω* are determined such that the dominant oscillation frequency is reproduced. Specifically, the approximation should match ${\widehat{D}}_{\mu}$ at the frequencies ${f}_{\mathrm{R}}={\text{argmax}}_{f}\phantom{\rule{0.166667em}{0ex}}\mid \text{Re}\left\{{\widehat{D}}_{\mu}\left(f\right)\right\}\mid $ and ${f}_{\mathrm{I}}={\text{argmax}}_{f}\phantom{\rule{0.166667em}{0ex}}\mid \text{Im}\left\{{\widehat{D}}_{\mu}\left(f\right)\right\}\mid $ in the Fourier domain as close as possible. We would like to note that using the method of least-squares over a range of frequencies instead can generate approximated filters which decay to zero instantly, particularly for large input mean and small variance (not shown). For such inputs a damped oscillator with a single decay time constant is too simple to fit the complete, rather complex linear filter shape. With Eq (86) the filtered mean input can be obtained by solving the second order ODE

$${\ddot{\mu}}_{\mathrm{f}}+\frac{2}{\tau}{\dot{\mu}}_{\mathrm{f}}+(\frac{2}{{\tau}^{2}}+{\omega}^{2})\phantom{\rule{1pt}{0ex}}{\mu}_{\mathrm{f}}=\frac{1+{\tau}^{2}{\omega}^{2}}{\tau}\phantom{\rule{1pt}{0ex}}(\frac{{\mu}_{\text{syn}}}{\tau}+{\dot{\mu}}_{\text{syn}}).$$

(87)

Using this damped oscillator filter gives rise to the LN_{dos} cascade model variant.

Considering the linear filter *D*_{σ} (Eq (80)) we use the approximation *D*_{σ}(*t*) ≈ *A*_{σ} exp(−*t*/*τ*_{σ}), with *A*_{σ} = 1/*τ*_{σ} and *τ*_{σ} determined by a least-squares fit of ${\widehat{D}}_{\sigma}$ for frequencies *f* [0, 1] kHz, as long as $\frac{\partial}{\partial \sigma}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})>0$. When this condition is not fulfilled, which occurs for large input mean and small variance, the full linear filter cannot be properly fit by an exponential function. This may be seen by the asymptotic behavior [77]

$$\underset{f\to \infty}{lim}{\widehat{R}}_{\sigma}\left(f\right)=\frac{{\sigma}_{\text{tot}}^{0}{r}_{\infty}\phantom{\rule{1pt}{0ex}}({\mu}_{\text{tot}}^{0},{\sigma}_{\text{tot}}^{0})}{i2\pi f{\Delta}_{\mathrm{T}}^{2}},$$

(88)

which implies negative ${lim}_{f\to \infty}{\widehat{D}}_{\sigma}\left(f\right)$ (cf. Eq (80)) that cannot be approximated for nonnegative *τ*_{σ}. In this case we use *τ*_{σ} → 0 which effectively yields *D*_{σ}(*t*) ≈ *δ*(*t*), justified by the observation that the full filter rapidly relaxes to zero (see Fig 8B). The linear filter application can be implemented by solving

$$\frac{d{\sigma}_{\mathrm{f}}}{dt}=\frac{{\sigma}_{\text{syn}}-{\sigma}_{\mathrm{f}}}{{\tau}_{\sigma}},$$

(89)

where *σ*_{f}(*t*) = *D*_{σ} * *σ*_{syn}(*t*) is the filtered input standard deviation. This filter is used in both model variants (LN_{exp} and LN_{dos}). Note (again) that *τ*_{σ} depends on ${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$.

In the derivation above we considered that the synaptic input mean and standard deviation, *μ*_{syn}(*t*) and *σ*_{syn}(*t*), vary around ${\mu}_{\text{syn}}^{0}$ and ${\sigma}_{\text{syn}}^{0}$ with small magnitudes. To extend the LN cascade model(s) to inputs that show large deviations from their baseline values we let the linear filters adjust to a changing input baseline in the following way: using the exponentially decaying mean input filter (LN_{exp} model) the filter parameters *τ*_{μ} and *τ*_{σ} are evaluated at

(90)

in every time step, i.e., these parameters adapt to the *effective* input moments. Using the damped oscillating mean input filter (LN_{dos} model), on the other hand, the filter parameters *τ*, *ω* and *τ*_{σ} are evaluated at *μ*_{tot}(*t*), *σ*_{tot}(*t*) given by Eqs (26) and (27), i.e., these parameters adapt directly to the *total* input moments, assuming that these moments do not fluctuate too vigorously. Note that because of these adjustments we do not need to consider a (particular) input baseline. Two remarks are in place: (i) the parameters of the damped oscillator cannot be adapted to a changing input baseline using the *effective* input mean (with *μ*_{f} given by Eq (87)), because this can lead to stable oscillations (for an uncoupled population) and thus decreased reproduction performance (not shown); (ii) for input moments that change very rapidly the reproduction performance of the LN_{dos} model variant may be improved by alternatively evaluating the parameters *τ*, *ω* and *τ*_{σ} at ${\mu}_{a}\left(t\right)={\mu}_{\mathrm{f}}^{a}-\langle w\rangle /C$, *σ*_{eff}(*t*), with ${\mu}_{\mathrm{f}}^{a}$ governed by Eq (85) (combining LN_{exp} and LN_{dos}), cf. [23] (ch. 4.2).

Finally, recurrent coupling within the population is included (in both model variants) by the dependence of the synaptic input moments on the delayed spike rate, *μ*_{syn}(*t*, *r*_{d}), ${\sigma}_{\text{syn}}^{2}(t,{r}_{d})$, with *r*_{d} given by Eq (40). For the LN_{dos} model we can then replace ${\dot{\mu}}_{\text{syn}}$ in Eq (87) by

$${\dot{\mu}}_{\text{syn}}={\dot{\mu}}_{\text{ext}}+JK{\dot{r}}_{d}={\dot{\mu}}_{\text{ext}}+JK\frac{r-{r}_{d}}{{\tau}_{d}}.$$

(91)

In case of identical (constant) propagation delays within the population this term would be ${\dot{\mu}}_{\text{syn}}={\dot{\mu}}_{\text{ext}}+JK\dot{r}(t-d)$ and in case of recurrent coupling without delays we would have

$${\dot{\mu}}_{\text{syn}}={\dot{\mu}}_{\text{ext}}+JK\phantom{\rule{1pt}{0ex}}\left[\frac{\partial \phantom{\rule{0.166667em}{0ex}}{r}_{\infty}}{\partial \mu}\phantom{\rule{1pt}{0ex}}\right({\dot{\mu}}_{\mathrm{f}}-\frac{a({\u27e8V\u27e9}_{\infty}-{E}_{w})-\u27e8w\u27e9}{C{\tau}_{w}}+\frac{b\phantom{\rule{0.166667em}{0ex}}r}{C})+\frac{\partial \phantom{\rule{0.166667em}{0ex}}{r}_{\infty}}{\partial \sigma}\frac{{\sigma}_{\text{syn}}-{\sigma}_{\mathrm{f}}}{{\tau}_{\sigma}}].$$

(92)

To summarize both LN cascade models, the population spike rate and mean membrane voltage are described by Eqs (70) and (73), using Eqs (78), (81) and (23), respectively. For the LN_{exp} model input filtering is governed by Eqs (85) and (89), where the filter parameters are evaluated at *μ*_{eff}(*t*), *σ*_{eff}(*t*) (Eq (90)). For the LN_{dos} model input filtering is described by Eqs (87) and (89), where the filter parameters are evaluated at *μ*_{tot}(*t*), *σ*_{tot}(*t*) given by Eqs (26) and (27). The system for the recurrent network under consideration is closed by Eqs (21) and (40) which relate population spike rate output and overall synaptic input moments. Note that both LN systems here are fully equivalent to the respective ones specified in the Sect. *Model reduction*.

In order to efficiently simulate the derived cascade rate models it is highly recommended to precalculate the quantities which are needed in each time step on a (*μ*, *σ*)-rectangle (see Sect. *Low-Dimensional Approximations* above). For the LN_{exp} model these quantities are the filter time constants *τ*_{μ}, *τ*_{σ}, and for the LN_{dos} model we need the quantities *τ*, *ω* and *τ*_{σ} (all displayed in Fig 8C). Both models additionally require the steady-state quantities *r*_{∞} and *V*_{∞} shown in Fig 6. An efficient implementation to obtain these quantities using Python with the package Numba for low-level virtual machine acceleration is available. Recall, that changing any parameter value of the external input, the recurrent synaptic input or the adaptation current does not require renewed precomputations.

If desired, it is also possible to obtain initial values for the variables of the cascade models (LN_{exp} and LN_{dos} variants) that correspond to a given initial distribution of membrane voltage and adaptation current values {*V*_{i}(0)}, {*w*_{i}(0)} of a population of *N* aEIF neurons. We can calculate *w*(0) = 1/*N*∑_{i}
*w*_{i}(0) and determine *μ*_{f}(0), *σ*_{f}(0) by requiring that the initial membrane voltage distribution of the respective LN model *p*_{∞}(*V*; *μ*_{f}(0) − *w*(0)/*C*, *σ*_{f}(0)) matches the initial voltage distribution from the aEIF population as close as possible (e.g., using the Kolmogorov–Smirnov statistic). For the LN_{dos} model we additionally set ${\dot{\mu}}_{\mathrm{f}}\left(0\right)=0$. This means we assume vanishing changes in the input history which underlies the initial membrane voltage distribution and filter parameters in the LN models (i.e., ${\dot{\mu}}_{\text{syn}}\approx 0$, ${\dot{\sigma}}_{\text{syn}}\approx 0$ for a sufficiently long time interval prior to *t* = 0).

The components of the LN model are derived in the limit of small amplitude variations of *μ*_{syn} and *σ*_{syn}. However, the approximation also provides an exact description of the population dynamics for very slow variations of *μ*_{syn} and *σ*_{syn}, where the spike rate, mean membrane voltage and adaptation current are well approximated by their steady-state values in each time step.

Here we approximated the derived linear filters using exponential and damped oscillator functions. We would like to note that, for a given baseline input (${\mu}_{\text{tot}}^{0}$, ${\sigma}_{\text{tot}}^{0}$) the filter application using the latter function (Eq (86)) can be equivalently described by a complex-valued ODE [23] (ch. 4.2). Furthermore, the true linear filters *D*_{μ} and *D*_{σ} can be substantially better approximated by a damped oscillator function with two time scales (i.e., two exponentials) each. In these three cases, however, the ODE representation for the filter application can lead to decreased reproduction performance when the baseline input changes very rapidly (due to increased sensitivity to variations of the filter parameters).

In this contribution the membrane voltage spike shape has been neglected (typical for IF type neuron models) by clamping *V*_{i} and *w*_{i} during the refractory period, justified by the observation that it is rather stereotyped and its duration is very brief. Furthermore, the spike shape is believed to contain little information compared to the time at which the spike occurs. Nevertheless, it can be incorporated in the aEIF model in a straightforward way using the following reset condition, as suggested in [43]: When *V*_{i} reaches the spike voltage *V*_{s} from below we let *V*_{i} decrease linearly from *V*_{s} to *V*_{r} during the refractory period and increment the adaptation current *w*_{i} ← *w*_{i} + *b* at the onset of that period. That is, *V*_{i} and *w*_{i} are not clamped during the refractory period, instead, *V*_{i} has a fixed time course and *w*_{i} is incremented by *b* and then governed again by Eq (15). This modification implies that the average membrane voltage in Eq (23) needs to be calculated over all neurons (and not only the nonrefractory ones), that is, *V* is calculated with respect to *p* + *p*_{ref}, where ${p}_{\text{ref}}(V,t)={\int}_{0}^{{T}_{\text{ref}}}r(t-s)\delta (V-{V}_{\text{sp}}\left(s\right))ds$ with spike trajectory *V*_{sp}(*t*) = *V*_{s} + (*V*_{r} − *V*_{s})*t*/*T*_{ref}, cf. [43]. The same applies to the steady-state mean membrane potential in Eqs (1), (39) and (76), i.e., *V*_{∞} is then given by

$${\u27e8V\u27e9}_{\infty}={\int}_{-\infty}^{{V}_{\mathrm{s}}}v{p}_{\infty}\left(v\right)dv+(1-{\int}_{-\infty}^{{V}_{\mathrm{s}}}{p}_{\infty}\left(v\right)dv)\phantom{\rule{1pt}{0ex}}\frac{{V}_{\mathrm{r}}+{V}_{\mathrm{s}}}{2},$$

(93)

instead of Eq (38). Notably, the accuracy of the adiabatic approximation (Eq (15)) does not depend on the refractory period *T*_{ref} in this case. That type of spike shape can therefore be considered in the FP model and the low-dimensional models in a straightforward way without significant additional computational demand. Note, however, that for the spec_{2} model a nonzero refractory period is not supported (see above). For an evaluation of the spike shape extension in terms of reproduction accuracy of the LN models see [23] (Fig. 4.15 in [23]).

A) Numerical integration of the time-dependent Fokker-Planck equation. B) Derivation of the model spec_{2} based on the Fokker-Planck operator. C) Numerical solver for the nonlinear Fokker-Planck eigenvalue problem.

(PDF)

Click here for additional data file.^{(266K, pdf)}

We thank Maurizio Mattia for several fruitful discussions, sharing his unpublished manuscript on the advanced spectral model and the permission to apply it to our system. Additionally we thank Srdjan Ostojic for providing his insights and code related to the eigenvalue problem of the Fokker-Planck operator. We would finally like to acknowledge that Douglas J. Sterling wrote the initial version of the software framework which was used to integrate and compare the different models.

The author JL was funded via German Research Foundation (Collaborative Research Center no. 910), http://www.itp.tu-berlin.de/collaborative_research_center_910/. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

Numerical implementations of the considered models written in the Python programming language are available under a free license at: https://github.com/neuromethods/fokker-planck-based-spike-rate-models.

1.
Shadlen MN, Newsome WT. The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J Neurosci. 1998; 18(10):3870–3896. [PubMed]

2.
Gerstner W, Kistler WM, Naud R, Paninski L. Neuronal dynamics: from single neurons to networks and models of cognition. Cambridge University Press; 2014.

3.
Harris KD, Shepherd GMG. The neocortical circuit: themes and variations. Nat Neurosci. 2015; 18(2):170–181. doi: 10.1038/nn.3917
[PMC free article] [PubMed]

4.
Okun M, Steinmetz NA, Cossell L, Iacaruso MF, Ko H, Barthó P, et al.
Diverse coupling of neurons to populations in sensory cortex. Nature. 2015; 521(7553):511–515. doi: 10.1038/nature14273
[PMC free article] [PubMed]

5.
Brette R, Gerstner W. Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J Neurophysiol. 2005; 94(5):3637–3642. doi: 10.1152/jn.00686.2005
[PubMed]

6.
Izhikevich EM. Simple model of spiking neurons. IEEE Trans Neural Netw. 2003; 14(6):1569–1572. doi: 10.1109/TNN.2003.820440
[PubMed]

7.
Badel L, Lefort S, Brette R, Petersen CCH, Gerstner W, Richardson MJE. Dynamic I-V curves are reliable predictors of naturalistic pyramidal-neuron voltage traces. J Neurophysiol. 2008; 99(2):656–666. doi: 10.1152/jn.01107.2007
[PubMed]

8. Naud R. The dynamics of adapting neurons. PhD Thesis, École Polytechnique Fédérale de Lausanne; 2011.

9.
Harrison PM, Badel L, Wall MJ, Richardson MJE. Experimentally verified parameter sets for modelling heterogeneous neocortical pyramidal-cell populations. PLOS Comput Biol. 2015; 11(8):e1004165
doi: 10.1371/journal.pcbi.1004165
[PMC free article] [PubMed]

10.
Rauch A, La Camera G, Lüscher HR, Senn W, Fusi S, Luscher HR. Neocortical pyramidal cells respond as integrate-and-fire neurons to in vivo-like input currents. J Neurophysiol. 2003; 90(3):1598–1612. doi: 10.1152/jn.00293.2003
[PubMed]

11.
Wilson HR, Cowan JD. Excitatory and inhibitory interactions in localized populations of model neurons. Biophys J. 1972; 12(1):1–24. doi: 10.1016/S0006-3495(72)86068-5
[PubMed]

12.
Latham PE, Richmond BJ, Nirenberg S, Nelson PG. Intrinsic dynamics in neuronal networks. I. Theory. J Neurophysiol. 2000; 83(2):828–835. [PubMed]

13.
Gigante G, Deco G, Marom S, Del Giudice P. Network events on multiple space and time scales in cultured neural networks and in a stochastic rate model. PLOS Comput Biol. 2015; 11(11):e1004547
doi: 10.1371/journal.pcbi.1004547
[PMC free article] [PubMed]

14.
Naud R, Gerstner W. Coding and decoding with adapting neurons: a population approach to the peri-stimulus time histogram. PLOS Comput Biol. 2012; 8(10):e1002711
doi: 10.1371/journal.pcbi.1002711
[PMC free article] [PubMed]

15.
Schwalger T, Deger M, Gerstner W. Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size. PLOS Comput Biol. 2017; 13(4):e1005507
doi: 10.1371/journal.pcbi.1005507
[PMC free article] [PubMed]

16.
Augustin M, Ladenbauer J, Obermayer K. How adaptation shapes spike rate oscillations in recurrent neuronal networks. Front Comput Neurosci. 2013; 7(9):1–11. [PMC free article] [PubMed]

17.
Nicola W, Ly C, Campbell SA. One-dimensional population density approaches to recurrently coupled networks of neurons with noise. SIAM J Appl Math. 2015; 75(5):2333–2360. doi: 10.1137/140995738

18.
Mattia M, Del Giudice P. Population dynamics of interacting spiking neurons. Phys Rev E. 2002; 66(5):051917
doi: 10.1103/PhysRevE.66.051917 [PubMed]

19.
Schaffer ES, Ostojic S, Abbott LF. A complex-valued firing-rate model that approximates the dynamics of spiking networks. PLOS Comput Biol. 2013; 9(10):e1003301
doi: 10.1371/journal.pcbi.1003301
[PMC free article] [PubMed]

20.
Mattia M. Low-dimensional firing rate dynamics of spiking neuron networks. arXiv Prepr. 2016; 1609.08855(1):1–8.

21.
Ostojic S, Brunel N. From spiking neuron models to linear-nonlinear models. PLOS Comput Biol. 2011; 7(1):e1001056
doi: 10.1371/journal.pcbi.1001056
[PMC free article] [PubMed]

22.
Brunel N, Wang XJ. What determines the frequency of fast network oscillations with irregular neural discharges? I. synaptic dynamics and excitation-inhibition balance. J Neurophysiol. 2003; 90(1):415–430. doi: 10.1152/jn.01095.2002
[PubMed]

23. Ladenbauer J. The collective dynamics of adaptive neurons: insights from single cell and network models. PhD Thesis, Technische Universität Berlin; 2015.

24. Lam SK, Pitrou A, Seibert S. Numba: A LLVM-based python JIT compiler. Proc Second Work LLVM Compil Infrastruct HPC–LLVM’15. 2015; p. 1–6.

25.
Ladenbauer J, Augustin M, Obermayer K. Intrinsic control mechanisms of neuronal network dynamics In: Schöll E, Klapp SHL, Hövel P, editors. Control of self-organizing nonlinear systems. Springer; 2016. p. 441–460.

26.
Gigante G, Mattia M, Del Giudice P. Diverse population-bursting modes of adapting spiking neurons. Phys Rev Lett. 2007; 98(14):148101
doi: 10.1103/PhysRevLett.98.148101
[PubMed]

27. Biggio M, Storace M, Mattia M. Equivalence between synaptic current dynamics and heterogeneous propagation delays in spiking neuron networks. arXiv Prepr. 2017; 1704.02780(1):1–14.

28.
Roxin A, Montbrió E. How effective delays shape oscillatory dynamics in neuronal networks. Physica D. 2011; 240(3):323–345. doi: 10.1016/j.physd.2010.09.009

29.
Richardson M. Effects of synaptic conductance on the voltage distribution and firing rate of spiking neurons. Phys Rev E. 2004; 69(5):051918
doi: 10.1103/PhysRevE.69.051918 [PubMed]

30.
Destexhe A, Mainen ZF, Sejnowski TJ. An efficient method for computing synaptic conductances based on a kinetic model of receptor binding. Neural Comput. 1994; 6(1):14–18. doi: 10.1162/neco.1994.6.1.14

31.
Schwalger T, Droste F, Lindner B. Statistical structure of neural spiking under non-Poissonian or other non-white stimulation. J Comput Neurosci. 2015; 39(1):29–51. doi: 10.1007/s10827-015-0560-x
[PubMed]

32.
Hertäg L, Durstewitz D, Brunel N. Analytical approximations of the firing rate of an adaptive exponential integrate-and-fire neuron in the presence of synaptic noise. Front Comput Neurosci. 2014; 8(116):1–22. [PMC free article] [PubMed]

33.
Ostojic S. Inter-spike interval distributions of spiking neurons driven by fluctuating inputs. J Neurophysiol. 2011; 106(1):361–373. doi: 10.1152/jn.00830.2010
[PubMed]

34.
Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci. 2000; 8(3):183–208. doi: 10.1023/A:1008925309027
[PubMed]

35.
Apfaltrer F, Ly C, Tranchina D. Population density methods for stochastic neurons with realistic synaptic kinetics: firing rate dynamics and fast computational methods. Network. 2006; 17(4):373–418. doi: 10.1080/09548980601069787
[PubMed]

36.
Nesse WH, Borisyuk A, Bressloff PC. Fluctuation-driven rhythmogenesis in an excitatory neuronal network with slow adaptation. J Comput Neurosci. 2008; 25(2):317–333. doi: 10.1007/s10827-008-0081-y
[PubMed]

37.
Zerlaut Y, Destexhe A. A mean-field model for conductance-based networks of adaptive exponential integrate-and-fire neurons. arXiv Prepr. 2017; 1703.00698(1):1–23.

38.
Buchin A, Chizhov A. Firing-rate model of a population of adaptive neurons. Biophysics. 2010; 55(4):664–673. doi: 10.1134/S0006350910040135 [PubMed]

39.
Montbrió E, Pazó D, Roxin A. Macroscopic description for networks of spiking neurons. Phys Rev X. 2015; 5(2):021028.

40.
Zhang JW, Rangan AV. A reduction for spiking integrate-and-fire network dynamics ranging from homogeneity to synchrony. J Comput Neurosci. 2015; 38(2):355–404. doi: 10.1007/s10827-014-0543-3
[PubMed]

41.
Gerstner W. Population dynamics of spiking neurons: fast transients, ascynchronous states, and locking. Neural Comput. 2000; 12(1):43–89. doi: 10.1162/089976600300015899
[PubMed]

42.
Iyer R, Menon V, Buice M, Koch C, Mihalas S. The influence of synaptic weight distribution on neuronal population dynamics. PLOS Comput Biol. 2013; 9(10):e1003248
doi: 10.1371/journal.pcbi.1003248
[PMC free article] [PubMed]

43.
Richardson MJE. Dynamics of populations and networks of neurons with voltage-activated and calcium-activated currents. Phys Rev E. 2009; 80(2):021928
doi: 10.1103/PhysRevE.80.021928 [PubMed]

44.
Rosenbaum R. A diffusion approximation and numerical methods for adaptive neuron models with stochastic inputs. Front Comput Neurosci. 2016; 10(39):1–20. [PMC free article] [PubMed]

45.
Nykamp DQ, Tranchina D. A population density approach that facilitates large-scale modeling of neural networks: extension to slow inhibitory synapses. Neural Comput. 2001; 13:511–546. doi: 10.1162/089976601300014448
[PubMed]

46.
Richardson MJE, Swarbrick R. Firing-rate response of a neuron receiving excitatory and inhibitory synaptic shot noise. Phys Rev Lett. 2010; 105(17):178102
doi: 10.1103/PhysRevLett.105.178102
[PubMed]

47.
Renart A, Moreno-Bote R, Wang XJ, Parga N. Mean-driven and fluctuation-driven persistent activity in recurrent networks. Neural Comput. 2007; 19(1):1–46. doi: 10.1162/neco.2007.19.1.1
[PubMed]

48.
Ladenbauer J, Augustin M, Obermayer K. How adaptation currents change threshold, gain and variability of neuronal spiking. J Neurophysiol. 2014; 111(5):939–953. doi: 10.1152/jn.00586.2013
[PubMed]

49.
Rosenbaum R, Smith MA, Kohn A, Rubin JE, Doiron B. The spatial structure of correlated neuronal variability. Nat Neurosci. 2016; 20(1):1–10. doi: 10.1038/nn.4433 [PMC free article] [PubMed]

50.
Destexhe A. Self-sustained asynchronous irregular states and up-down states in thalamic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. J Comput Neurosci. 2009; 27(3):493–506. doi: 10.1007/s10827-009-0164-4
[PubMed]

51.
Wang XJ, Liu YH, Sanchez-Vives MV, Mccormick DA. Adaptation and temporal decorrelation by single neurons in the primary visual cortex. J Neurophysiol. 2003; 89(6):3279–3293. doi: 10.1152/jn.00242.2003
[PubMed]

52.
Stimberg M, Goodman DFM, Benichoux V, Brette R. Equation-oriented specification of neural models for simulations. Front Neuroinform. 2014; 8(6):1–14. [PMC free article] [PubMed]

53.
Goodman DFM, Brette R. The brian simulator. Front Neurosci. 2009; 3(2):192–197. doi: 10.3389/neuro.01.026.2009
[PMC free article] [PubMed]

54.
Brown DA, Adams PR. Muscarinic suppression of a novel voltage-sensitive K+ current in a vertebrate neurone. Nature. 1980; 283(5748):673–676. doi: 10.1038/283673a0
[PubMed]

55.
Sanchez-Vives MV, McCormick DA. Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat Neurosci. 2000; 3(10):1027–1034. doi: 10.1038/79848
[PubMed]

56.
Sanchez-Vives MV, Nowak LG, McCormick DA. Cellular mechanisms of long-lasting adaptation in visual cortical neurons in vitro. J Neurosci. 2000; 20(11):4286–4299. [PubMed]

57.
Stocker M. Ca(2+)-activated K+ channels: molecular determinants and function of the SK family. Nat Rev Neurosci. 2004; 5(10):758–770. doi: 10.1038/nrn1516
[PubMed]

58.
Brunel N, Hakim V, Richardson MJE. Firing-rate resonance in a generalized integrate-and-fire neuron with subthreshold resonance. Phys Rev E. 2003; 67(5):051916
doi: 10.1103/PhysRevE.67.051916 [PubMed]

59.
Touboul J. Limits and dynamics of stochastic neuronal networks with random heterogeneous delays. J Stat Phys. 2012; 149(4):569–597. doi: 10.1007/s10955-012-0607-6

60.
Delarue F, Inglis J, Rubenthaler S, Tanré E. Global solvability of a networked integrate-and-fire model of McKean-Vlasov type. Ann Appl Probab. 2015; 25(4):2096–2133. doi: 10.1214/14-AAP1044

61.
Marpeau F, Barua A, Josić K. A finite volume method for stochastic integrate-and-fire models. J Comput Neurosci. 2009; 26(3):445–457. doi: 10.1007/s10827-008-0121-7
[PubMed]

62.
Cáceres MJ, Carrillo JA, Tao L. A numerical solver for a nonlinear Fokker-Planck equation representation of neuronal network dynamics. J Comput Phys. 2011; 230(4):1084–1099. doi: 10.1016/j.jcp.2010.10.027

63.
Scharfetter DL, Gummel HK. Large-signal analysis of a silicon tead diode oscillator. IEEE Trans Electron Dev. 1969; 16(1):64–77. doi: 10.1109/T-ED.1969.16566

64.
Gosse L. Computing qualitatively correct approximations of balance laws. vol. 2
Springer; 2013.

65.
Farrell PA, Gartland EC Jr. On the Scharfetter-Gummel discretization for drift-diffusion continutity equations In: Miller JJH, editor. Computational methods for boundary and interior layers in several dimensions. Boole Press; 1991. p. 51–79.

66.
LeVeque RJ. Finite volume methods for hyperbolic problems. Cambridge University Press; 2002.

67.
Oliphant TE. Python for scientific computing. Comput Sci Eng. 2007; 9(3):10–20. doi: 10.1109/MCSE.2007.58

68.
MacDonald N. Time lags in biological models In: Levin SA, editor. Lecture Notes in Biomath. 27th ed
Springer; 1978. p. 1–114.

69.
Risken H. The Fokker-Planck equation: methods of solution and applications. Springer; 1996.

70.
Knight BW, Manin D, Sirovich L. Dynamical models of interacting neuron populations
Proc Symp Robot Cybern Lille-France. 1996; p. 4–8.

71.
Hochbruck M, Ostermann A. Exponential integrators. Acta Numer. 2010; 19(5):209–286. doi: 10.1017/S0962492910000048

72.
Bernstein DS, So W. Some explicit formulas for the matrix exponential. IEEE Trans Autom Control. 1993; 38(8):1228–1232. doi: 10.1109/9.233156

73.
Richardson MJE. Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive. Phys Rev E. 2007; 76(2):021919
doi: 10.1103/PhysRevE.76.021919 [PubMed]

74.
Chichilnisky EJ. A simple white noise analysis of neuronal light responses. Network. 2001; 12(2):199–213. doi: 10.1080/713663221
[PubMed]

75.
Pillow JW, Shlens J, Paninski L, Sher A, Litke AM, Chichilnisky EJ, et al.
Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature. 2008; 454(7207):995–999. doi: 10.1038/nature07140
[PMC free article] [PubMed]

76.
Fourcaud-Trocmé N, Hansel D, van Vreeswijk C, Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating inputs. J Neurosci. 2003; 23(37):11628–11640. [PubMed]

77.
Fourcaud-Trocmé N, Brunel N. Dynamics of the instantaneous firing rate in response to changes in input statistics. J Comput Neurosci. 2005; 18(3):311–321. doi: 10.1007/s10827-005-0337-8
[PubMed]

Articles from PLoS Computational Biology are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |