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

**|**J Math Neurosci**|**v.7; 2017**|**PMC5725415

Formats

Article sections

- Abstract
- Introduction
- Single Spiking Neuron Model
- The Rate-Reduced Neuron Model
- Augmenting Neural Fields
- Discussion
- References

Authors

Related links

J Math Neurosci. 2017; 7: 13.

Published online 2017 December 11. doi: 10.1186/s13408-017-0055-3

PMCID: PMC5725415

Koen Dijkstra,^{}^{1} Yuri A. Kuznetsov,^{1,}^{2} Michel J. A. M. van Putten,^{3,}^{4} and Stephan A. van Gils^{1}

Koen Dijkstra, Email: ln.etnewtu@artskjid.neok.

Received 2017 March 29; Accepted 2017 November 21.

Copyright © The Author(s) 2017

We present a simple rate-reduced neuron model that captures a wide range of complex, biologically plausible, and physiologically relevant spiking behavior. This includes spike-frequency adaptation, postinhibitory rebound, phasic spiking and accommodation, first-spike latency, and inhibition-induced spiking. Furthermore, the model can mimic different neuronal filter properties. It can be used to extend existing neural field models, adding more biological realism and yielding a richer dynamical structure. The model is based on a slight variation of the Rulkov map.

Networks of coupled neurons quickly become analytically intractable and computationally infeasible due to their large state and parameter spaces. Therefore, starting with the work of Beurle [1], a popular modeling approach has been the development of continuum models, called neural fields, that describe the average activity of large populations of neurons (Wilson and Cowan [2, 3], Nunez [4], Amari [5, 6]). In neural field models, the network architecture is represented by connectivity functions and the corresponding transmission delays, while differential operators characterize synaptic dynamics. All intrinsic properties of the underlying neuronal populations are condensed into firing rate functions, which replace individual neuronal action potentials and map the sum of all incoming synaptic currents to an outgoing firing rate. While some neural field models incorporate spike-frequency adaptation (Pinto and Ermentrout [7, 8], Coombes and Owen [9], Amari [10, 11]), more complex spiking behavior such as postinhibitory rebound, phasic spiking and accommodation, first-spike latency, and inhibition-induced spiking is mostly absent, an exception being the recent reduction of the Izhikevich neuron (Nicola and Campbell [12], Visser and van Gils [13]).

Here, we present a rate-reduced model that is based on a slight modification of the Rulkov map (Rulkov [14], Rulkov et al. [15]), a phenomenological, map-based single neuron model. Similar to Izhikevich neurons (Izhikevich [16]), the Rulkov map can mimic a wide variety of biologically realistic spiking patterns, all of which are preserved by our rate formulation. The rate-reduced model can therefore be used to incorporate all the aforementioned types of spiking behavior into existing neural field models.

This paper is organized as follows. In Sect. 2, we present the single spiking neuron model our rate-reduced model is based upon, and illustrate different spiking patterns and filter properties. In Sect. 3 we heuristically reduce the single neuron model to a rate-based formulation, and show that the rate-reduced model preserves spiking and filter properties. We give an example of a neural field that is augmented with our rate model in Sect. 4 and end with a discussion in Sect. 5.

In this section we present a phenomenological, map-based single neuron model, which is a slight modification of the Rulkov map (Rulkov [14], Rulkov et al. [15]). The Rulkov map was designed to mimic the spiking and spiking-bursting activity of many real biological neurons. It has computational advantages because the map is easier to iterate than continuous dynamical systems. Furthermore, as we will show in this paper, it is straightforward to obtain a rate-reduced version of a slightly modified version of the Rulkov model.

The Rulkov map consists of a fast variable *v*, resembling the neuronal membrane potential, and a slow adaptation variable *a*. In our modification of the original model, the adaptation only implicitly depends on the membrane potential through a binary spiking variable. As we will show in the next section, this modification allows for an easy decoupling of the membrane potential and adaptation variable, and therefore a straightforward rate reduction of the model. The cost of the modification is the loss of subthreshold oscillation dynamics. The modified Rulkov map is given by

$$\{\begin{array}{c}{v}_{n+1}=f({v}_{n},{v}_{n-1},\kappa {u}_{n}-{a}_{n}-\theta ),\hfill \\ {a}_{n+1}={a}_{n}-\epsilon ({a}_{n}+(1-\kappa ){u}_{n}-\gamma {s}_{n}),\hfill \end{array}$$

SNM

where the piecewise continuous function *f*:ℝ^{3} → ℝ is given by

$$f({x}_{1},{x}_{2},{x}_{3})=\{\begin{array}{cc}\frac{2500+150{x}_{1}}{50-{x}_{1}}+50{x}_{3}\hfill & \text{if}{x}_{1}0,\hfill \\ 50+50{x}_{3}\hfill & \text{if}0\le {x}_{1}50+50{x}_{3}\phantom{\rule{1em}{0ex}}\wedge \phantom{\rule{1em}{0ex}}{x}_{2}0,\hfill \\ -50\hfill & \text{otherwise}.\hfill \end{array}$$

1

The form of *f* is chosen to mimic the shape of neuronal action potentials. The variable *u* in (SNM) represents external (synaptic) input to the cell, which we assume to be given, and *s* is a binary indicator variable, given by

$${s}_{n}=\{\begin{array}{cc}1\hfill & \text{if the neuron spiked at iteration}n\text{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

2

A Rulkov neuron spikes at iteration *n* if its membrane potential is reset to *v*_{n+1} = −50 in the next iteration. It follows from (1) that the spiking condition in (2) is satisfied if and only if

3

The dependence of *v*_{n+1} on *v*_{n−1} in (SNM) ensures that a neuron always spikes if its membrane potential is non-negative for two consecutive iterations, independent of the external input *u*. To mimic spiking patterns of real biological neurons, one time step should correspond to approximately 0.5 ms of time.

The parameter 0 < *ε* < 1 in (SNM) sets the time scale of the adaptation variable and *γ* determines the adaptation strength. The parameter *θ* can be interpreted as a spiking threshold: for constant external input *u*_{n} = *φ*, the neuron spikes persistently if and only if *φ* > *θ*. After a change of variable *a*_{n} → *a*_{n} + (1 − *κ*)*φ* and parameter *θ* → *θ* − *φ*, constant external input vanishes. Therefore, the asymptotic response to constant input does not depend on the parameter *κ*. However, the parameter *κ* can be used to tune the transient response of the neuron to changes in external input, as it determines how input is divided between the fast and the slow subsystem of (SNM). For parameter values *κ* ∈ [0, 1], *κ* can be interpreted as the fraction of the input that is applied to the fast subsystem, and therefore determines (together with *ε*) how quickly the membrane potential dynamics react to changes in input. Since the effective drive of the system is given by *κ**u*_{n} − *a*_{n}, changes in external input are initially magnified for *κ* > 0. Asymptotically, this is then counterbalanced by additional adaption. Finally, for *κ* < 0, the initial response of the membrane potential to a change in input is reversed, i.e. an increase in external input initially has an inhibitory effect, and a decrease in external input initially has an excitatory effect.

The Rulkov map (SNM) with 0 < *ε* ≪ 1 is a slow-fast system, and we can explore the fast spiking dynamics of the model by assuming the suprathreshold drive *κ**u*_{n} − *a*_{n} − *θ* = *ς* is constant. In this case, (SNM) reduces to the fast subsystem

$${v}_{n+1}=\{\begin{array}{cc}\frac{2500+150{v}_{n}}{50-{v}_{n}}+50\varsigma \hfill & \text{if}{v}_{n}0,\hfill \\ 50+50\varsigma \hfill & \text{if}0\le {v}_{n}50+50\varsigma ,\hfill \\ -50\hfill & \text{otherwise}.\hfill \end{array}$$

FSS

The map (FSS) undergoes a saddle-node bifurcation at *ς* = 0 (Fig. 1). For *ς* < 0 there exist a stable and an unstable fixed point, given by

$${v}_{\mathrm{s}}=25(\varsigma -2-\sqrt{{\varsigma}^{2}-8\varsigma}),\phantom{\rule{2em}{0ex}}{v}_{\mathrm{u}}=25(\varsigma -2+\sqrt{{\varsigma}^{2}-8\varsigma}),$$

4

respectively (Fig. 1A), while the system will settle into a stable periodic orbit for *ς* > 0 (Fig. 1B). In the former case the unstable fixed point acts as an excitation threshold: if the value of the membrane potential exceeds this point, it will spike once and then decay back to the stable equilibrium. Since the unstable fixed point *v*_{u} always lies to the right of the ‘reset potential’ *v* = −50, a stable fixed point and a periodic orbit can never coexist. This guarantees that we can define a firing rate function *S*:ℝ → ℚ for the fast subsystem (FSS), given by

$$S(\varsigma )=\{\begin{array}{cc}0\hfill & \text{for}\varsigma \le 0,\hfill \\ \frac{1}{P(\varsigma )}\hfill & \text{for}\varsigma 0,\hfill \end{array}$$

5

where *P*:ℝ_{>0} → ℕ maps the drive to the period of the corresponding stable limit cycle of (FSS). The fast subsystem (FSS) is piecewise-defined on the ‘left’ interval (−∞, 0), the ‘middle’ interval [0, 50 + 50*ς*), and the ‘right’ interval [50 + 50*ς*, ∞). The left interval is mapped to the left and middle interval, and the middle and right interval are mapped to right and left interval, respectively. The period of a limit cycle of (FSS) therefore only depends on the number of iterations in the left interval. Note, however, that the shape of the function *f* given in (1) can easily be changed to support bistability in the fast subsystem, which allows for some additional dynamics such as ‘chattering’, a response of periodic bursts of spikes to constant input (Rulkov [14]).

Izhikevich [17] classified different features of biological spiking neurons, most of which can be mimicked by our modified Rulkov model (SNM). In the following, we discuss the role of the model parameters with the help of a few physiologically relevant examples.

Tonically spiking (also called ‘fast spiking’) neurons respond to a step input with spike trains of constant frequency. Most inhibitory neurons are fast spiking (Izhikevich [17]). In the modified Rulkov model this can be achieved by choosing a ‘large’ $(1>\epsilon >\frac{1}{10})$ value for the time scale parameter, in which case the influence of a single spike on the adaptation variable decays very fast. Therefore, the value of the adaptation variable is dominated by the timing of the last spike and the influence of older spikes is negligible (Fig. 2A). Since the time scale separation is small, the qualitative dynamics does not depend on *κ*.

Different types of spiking patterns generated by the single neuron model (SNM). Corresponding parameter values (*θ*, *κ*, *ε*, *γ*) are given in brackets. **(**
**A**
**)** Tonic spiking $(\frac{1}{10},\frac{1}{2},\frac{1}{2},\frac{1}{2})$. **(**
**B**
**)** Spike-frequency adaptation $(\frac{1}{10},1,\frac{1}{1000},5)$. **(** **...**

Most cortical excitatory neurons are not ‘fast spiking’, but respond to a step input with a spike train of slowly decreasing frequency, a phenomenon known as ‘spike-frequency adaptation’ (also called ‘regular spiking’). This kind of spiking behavior can be modeled by applying all input to the fast subsystem (*κ* = 1) and choosing *ε* ≪ 1. The adaptation variable then acts as a slow time scale, such that a single spike has a long-lasting effect on the adaptation variable (Fig. 2B). The level of adaptation can be controlled with *γ*.

The excitability of some neurons is temporarily enhanced after they are released from hyperpolarizing current, which can result in the firing of one or more ‘rebound spikes’. Rebound spiking is an important mechanism for central pattern generation for heartbeat and other motor patterns in many neuronal systems (Chik et al. [18]). In the modified Rulkov map, postinhibitory rebound spiking can be modeled by choosing *κ* > 1. In this case, the adaptation variable will become negative while the cell gets hyperpolarized, which can be sufficient to trigger temporary spiking once the inhibitory input is turned off (Fig. 2C). Similarly, excitatory ‘subthreshold’ (*u*_{n} < *θ*) input can elicit temporary spiking if the input is ramped up sufficiently fast (Fig. 2D).

If all input is applied to the slow subsystem (*κ* = 0), there can be a large latency between the input onset and the first spike of the neuron, yielding a delayed response to a pulse input (Fig. 2E). For *κ* < 0, the initial response of the model to changes in input is reversed: excitation initially leads to hyperpolarization of the neuron and inhibition can induce temporary spiking (Fig. 2F). This inhibition-induced spiking is a feature of many thalamo-cortical neurons (Izhikevich [17]).

In the previous section, we illustrated how the parameter *κ* can tune transient spiking responses of the modified Rulkov map to changes in external input. In reality, neurons often receive strong periodic input, e.g. from a synchronous neuronal population nearby. Information transfer between neurons may be optimized by temporal filtering, which is especially important when the same signal transmits distinct messages (Blumhagen et al. [19]). In this section, we study the response of (SNM) to harmonic input

$${u}_{n}=\phi cos(\frac{\omega \pi n}{1000}+\vartheta ),$$

6

with amplitude , phase shift *ϑ* ∈ [0, 2*π*), and where *ω* ∈ [0, 1000] corresponds to the input frequency in Hz assuming that one iteration of (SNM) corresponds to 0.5 ms of time. A Rulkov neuron (SNM) will never spike if

7

In this case, the adaptation reduces to the simple linear equation

8

with explicit solution

$${a}_{n}=-\epsilon (1-\kappa )\sum _{m=1}^{\mathrm{\infty}}{(1-\epsilon )}^{m-1}{u}_{n-m}.$$

9

Inserting (6) into (9) now yields

$$\begin{array}{rl}\kappa {u}_{n}-{a}_{n}& =\kappa \phi cos(\frac{\omega \pi n}{1000}+\vartheta )\\ & \phantom{\rule{1em}{0ex}}+\epsilon (1-\kappa )\phi \sum _{m=1}^{\mathrm{\infty}}{(1-\epsilon )}^{m-1}cos(\frac{\omega \pi (n-m)}{1000}+\vartheta )\\ & =F(\omega )\frac{\phi}{2}{e}^{(\frac{\omega \pi n}{1000}+\vartheta )i}+\stackrel{\u203e}{F(\omega )}\frac{\phi}{2}{e}^{-(\frac{\omega \pi n}{1000}+\vartheta )i}\\ & =|F(\omega )|\phi cos(\frac{\omega \pi n}{1000}+\vartheta +arg(F(\omega ))),\end{array}$$

10

where the overline denotes complex conjugation and the frequency response *F*:[0, 1000] ↦ ℂ is given by

$$F(\omega )=\kappa +\frac{\epsilon (1-\kappa )}{{e}^{\frac{\omega \pi i}{1000}}+\epsilon -1}.$$

11

The absolute value and argument of the frequency response determine the relative magnitude and phase of the output, respectively. It follows that a Rulkov neuron (SNM) receiving periodic input (6) does not spike if

|*F*(*ω*)|*φ* ≤ *θ*.

12

The inverse statement is not true, even if *ω* and in (6) are chosen such that

$$cos(\frac{\omega \pi n}{1000}+\vartheta +arg(F(\omega )))=1\phantom{\rule{1em}{0ex}}\text{for some}n\in \mathbb{N}.$$

13

Since it can take a few iterations of the map to converge to its periodic orbit, a neuron will only spike if its drive is larger than the threshold *θ* for a sufficiently long time. The modulus of the frequency response (11) is given by

$$|F(\omega )|=\sqrt{\frac{{\epsilon}^{2}+2\kappa (\kappa -\epsilon )(1-cos(\frac{\omega \pi}{1000}))}{{\epsilon}^{2}+2(1-\epsilon )(1-cos(\frac{\omega \pi}{1000}))}},$$

14

and it follows that |*F*| is strictly decreasing if and only if *κ* ∈ (−1 + *ε*, 1), and increasing otherwise (Fig. 3). Clearly,

$$F(0)=1,\phantom{\rule{2em}{0ex}}F(1000)=\frac{2\kappa -\epsilon}{2-\epsilon}.$$

15

The input parameter *κ* can therefore be used to model filter properties of the neuron. For −1 + *ε* < *κ* < 1 high frequencies get attenuated and a neuron can act as a low-pass filter in the sense that periodic input within a certain amplitude range only elicits a spiking response if its frequency is low enough (Fig. 4A). Similarly, for *κ* > 1 (and *κ* < −1 + *ε*), high frequencies get amplified and there exists an amplitude range for which the neuron acts as a high-pass filter (Fig. 4B).

Neural field models are based on the assumption that neuronal populations convey all relevant information in their (average) firing rates. If one wants to incorporate certain spiking dynamics, one has to come up with a corresponding rate-reduced formulation first. In this section we present a rate-reduced version of the Rulkov model (SNM) that can be used to extend existing neural field models.

The adaptation variable *a* in the spiking neuron model (SNM) only implicitly depends on the membrane potential *v* via the binary spiking variable *s*. We can therefore decouple the adaption variable from the membrane potential by replacing the binary spiking variable defined in (2) by the instantaneous firing rate (5) of the fast subsystem (FSS), yielding

16

By interpreting (16) as the forward discretization of an ordinary differential equation, we arrive at the continuous time rate-reduced model

$$\frac{1}{\epsilon}\frac{\mathrm{d}a}{\mathrm{d}t}=-a-(1-\kappa )u+\gamma S(\kappa u-a-\theta ).$$

RNM

The rate-reduced neuron model (RNM) preserves the dynamical features of the full model (SNM) and reproduces all previous example spiking patterns (Fig. 5).

Analogously to Sect. 2.3, we now study the response of the rate-reduced model (RNM) to sinusoidal input

$$u(t)=\phi cos(\frac{\omega \pi t}{1000}+\vartheta ).$$

17

Under the assumption that

18

the explicit solution of (RNM) is given by

$$a(t)=-\epsilon (1-\kappa ){\int}_{-\mathrm{\infty}}^{t}{e}^{-\epsilon (t-\tau )}u(\tau )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\tau ,$$

19

cf. (9). Inserting the input (17) into (19) yields

$$\begin{array}{rl}\kappa u(t)-a(t)& =\kappa \phi cos(\frac{\omega \pi t}{1000}+\vartheta )+\epsilon (1-\kappa )\phi {\int}_{-\mathrm{\infty}}^{t}{e}^{-\epsilon (t-\tau )}cos(\frac{\omega \pi \tau}{1000}+\vartheta )\phantom{\rule{0.2em}{0ex}}\mathrm{d}\tau \\ & =G(\omega )\frac{\phi}{2}{e}^{(\frac{\omega \pi t}{1000}+\vartheta )i}+\stackrel{\u203e}{G(\omega )}\frac{\phi}{2}{e}^{-(\frac{\omega \pi t}{1000}+\vartheta )i}\\ & =|G(\omega )|\phi cos(\frac{\omega \pi t}{1000}+\vartheta +arg(G(\omega ))),\end{array}$$

20

where the frequency response *G*:ℝ_{≥0} ↦ ℂ is given by

$$G(\omega )=\kappa +\frac{\epsilon (1-\kappa )}{\epsilon +\frac{\omega \pi i}{1000}}.$$

21

It follows that for the rate-reduced model (RNM) receiving harmonic input (17) we have

22

Because we neglected the transient corresponding to the convergence from fixed point to limit cycle in the rate-reduced model (RNM), the inequality in (22) defines a clear ‘spiking condition’. The modulus of the frequency response (21) is given by

$$|G(\omega )|=\sqrt{\frac{{\epsilon}^{2}+{\kappa}^{2}{(\frac{\pi \omega}{1000})}^{2}}{{\epsilon}^{2}+{(\frac{\pi \omega}{1000})}^{2}}},$$

23

and |*G*| therefore is strictly decreasing if and only if |*κ*| ≤ 1, and increasing otherwise. Revisiting the examples from Sect. 2.3 (Fig. 4), we have

|*G*(1)|*φ* = 0.1696… > *θ* > 0.1255… = |*G*(2)|*φ*,

24

for $(\kappa ,\epsilon ,\theta ,\phi )=(\frac{1}{10},\frac{1}{200},\frac{1}{7},\frac{1}{5})$, and

|*G*(1)|*φ* = 0.1359… < *θ* < 0.1684… = |*G*(2)|*φ*,

25

for $(\kappa ,\epsilon ,\theta ,\phi )=(2,\frac{1}{200},\frac{1}{7},\frac{1}{10})$. Indeed, the rate-reduced model (RNM) reproduces the examples of the full model both qualitatively and quantitatively (Fig. 6). When the rate-reduced model (RNM) is incorporated into existing neural field models, the frequency response of the reduced model can be used to tune the individual temporal filter properties of the different neuronal populations.

Since our neuron model (SNM) is a map, the period *P* of its limit cycle lies in ℕ for all positive suprathreshold drives *ς*. Therefore, the spiking rate function (5) is staircase-like, with points of discontinuity whenever *P* → *P* + 1. Let {*ς*_{1}, *ς*_{2}, …} denote the set of all points of discontinuity of the firing rate function in decreasing order. For *ς* ≥ *ς*_{1} = 1 the ‘reset potential’ *v* = −50 in (FSS) is immediately mapped to a non-negative number, and the neuron is therefore spiking at its maximal frequency of once in three iterations. Similarly, the voltage stays in the left interval for two iterations and the neuron is spiking once in four iterations for ${\varsigma}_{1}>\varsigma \ge {\varsigma}_{2}=\frac{1}{2}(5-\sqrt{17})$. In general, at *ς*_{k}, there is a jump discontinuity of size

$$\underset{\varsigma \to {\varsigma}_{k}^{+}}{lim}S(\varsigma )-\underset{\varsigma \to {\varsigma}_{k}^{-}}{lim}S(\varsigma )=\frac{1}{(k+2)(k+3)},\phantom{\rule{1em}{0ex}}\text{with}S({\varsigma}_{k})=\frac{1}{k+2}.$$

26

The firing rate of the fast subsystem (FSS) can therefore be written as

$$S(\varsigma )=\sum _{k=1}^{\mathrm{\infty}}\frac{H(\varsigma -{\varsigma}_{k})}{(k+2)(k+3)},$$

27

where *H* is the Heaviside step function and

$$\underset{k\to \mathrm{\infty}}{lim}{\varsigma}_{k}=0.$$

28

In large neuronal networks, it is often assumed that the spiking thresholds of the individual neurons are randomly distributed. This ensures heterogeneity and models intrinsic interneuronal differences or random input from outside the network. If we add Gaussian noise to the threshold parameter *θ* in (SNM), it is natural to define an expected firing rate 〈*S*〉:ℝ ↦ ℝ, given by

$$\langle S(\varsigma )\rangle =\frac{1}{\sqrt{2\pi {\sigma}^{2}}}{\int}_{-\mathrm{\infty}}^{\mathrm{\infty}}{e}^{\frac{-{w}^{2}}{2{\sigma}^{2}}}S(\varsigma +w)\phantom{\rule{0.2em}{0ex}}\mathrm{d}w,$$

29

where *σ*^{2} is the variance of the noise. Using (27), we can rewrite (29) as

$$\langle S(\varsigma )\rangle =\frac{1}{6}+\sum _{k=1}^{\mathrm{\infty}}\frac{erf(\frac{\varsigma -{\varsigma}_{k}}{\sqrt{2{\sigma}^{2}}})}{2(k+2)(k+3)},$$

30

where erf denotes the error function. While *S*(*ς*) can readily be computed for any *ς* ∈ ℝ and we derived a concise expression for the expected firing rate, the infinite sum (30) cannot easily be evaluated. For this reason, we approximate 〈*S*(*ς*)〉 by a finite sum of the form

$$\frac{1}{6}+\frac{1}{6N}\sum _{i=1}^{N}erf\left(\frac{\varsigma -{\nu}_{i}}{{\chi}_{i}}\right),$$

31

for some fixed *N* ∈ ℕ and constants *ν*_{i}, *χ*_{i} ∈ ℝ, which are chosen by (numerically) minimizing

$${\parallel \frac{1}{\sqrt{2\pi {\sigma}^{2}}}{\int}_{-\mathrm{\infty}}^{\mathrm{\infty}}{e}^{\frac{-{w}^{2}}{2{\sigma}^{2}}}S(\varsigma +w)\phantom{\rule{0.2em}{0ex}}\mathrm{d}w-\frac{1}{6}-\frac{1}{6N}\sum _{i=1}^{N}erf\left(\frac{\varsigma -{\nu}_{i}}{{\chi}_{i}}\right)\parallel}_{2}.$$

32

For large noise levels *σ*^{2}, the average firing rate (29) has a sigmoidal shape and can be very well approximated with a small value of *N* (Fig. 7).

When large populations of neurons are modeled by networks of individual, interconnected cells, the high dimensionality of state and parameter spaces makes mathematical analysis intractable and numerical simulations costly. Moreover, large network simulations provide little insight into global dynamical properties. A popular modeling approach to circumventing the aforementioned problems is the use of neural field equations. These models aim to describe the dynamics of large neuronal populations, where spikes of individual neurons are replaced by (averaged) spiking rates and space is continuous. Another advantage of neural fields is that they are often well suited to model experimental data. In brain slice preparations, spiking rates can be measured with an extracellular electrode, while intracellular recordings are much more involved. Furthermore, the most common clinical measurement techniques of the brain, electroencephalography (EEG) and functional magnetic resonance imaging (fMRI), both represent the average activity of large groups of neurons and may therefore be better modeled by population equations. The first neural field model can be attributed to Beurle [1], however, the theory really took off with the work of Wilson and Cowan [2, 3], Amari [5, 6], and Nunez [4].

In ‘classical’ neural field models the firing rate of a neuronal population is assumed to be given by its instantaneous input, which is only valid for tonically spiking neurons. With the help of our rate-reduced model (RNM), it is straightforward to augment existing neural field models with more complex spiking behavior. As an example, we will look at the following two-population model on the one-dimensional spatial domain *Ω* = (−1, 1):

$$\begin{array}{rl}(1+\frac{1}{{\alpha}_{1}}\frac{\partial}{\partial t}){u}_{1}(t,x)& ={\int}_{-1}^{1}{J}_{11}(x,{x}^{\prime}){S}_{1}\left({r}_{1}(t,{x}^{\prime})\right)\\ & \phantom{\rule{1em}{0ex}}+{J}_{12}(x,{x}^{\prime}){S}_{2}\left({r}_{2}(t,{x}^{\prime})\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{x}^{\prime},\\ (1+\frac{1}{{\epsilon}_{1}}\frac{\partial}{\partial t}){a}_{1}(t,x)& =-(1-{\kappa}_{1}){u}_{1}(t,x)+{\gamma}_{1}{S}_{1}({r}_{1}(t,x)),\\ (1+\frac{1}{{\alpha}_{2}}\frac{\partial}{\partial t}){u}_{2}(t,x)& ={\int}_{-1}^{1}{J}_{21}(x,{x}^{\prime}){S}_{1}\left({r}_{1}(t,{x}^{\prime})\right)\\ & \phantom{\rule{1em}{0ex}}+{J}_{22}(x,{x}^{\prime}){S}_{2}\left({r}_{2}(t,{x}^{\prime})\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{x}^{\prime},\\ (1+\frac{1}{{\epsilon}_{2}}\frac{\partial}{\partial t}){a}_{2}(t,x)& =-(1-{\kappa}_{2}){u}_{2}(t,x)+{\gamma}_{2}{S}_{2}({r}_{2}(t,x)),\end{array}$$

ANF

where, as before,

33

for *i* ∈ {1, 2}. The differential operators in the left-hand side of the integral equations in (ANF) model exponentially decaying synaptic currents with decay rate *α*_{i}. The connectivity *J*_{ij}(*x*, *x*^{′}) measures the connection strength from neurons of population *j* and position *x*^{′} to neurons of population *i* and position *x*. The connectivity kernels ${J}_{ij}:\stackrel{\u203e}{\Omega}\times \stackrel{\u203e}{\Omega}\mapsto \mathbb{R}$ are assumed to be isotropic and given by

34

where *ρ*_{j} is the density of neurons of type *j*, *η*_{ij} is the maximal connection strength, and *μ*_{ij} is the spatial decay rate of the connectivity. Both firing rate functions *S*_{i}:ℝ ↦ ℝ are chosen to approximate the expected firing rate of Rulkov neurons (SNM) with a noise level of *σ*^{2} = ¼ (Fig. 7),

$${S}_{1}(\varsigma )={S}_{2}(\varsigma )=\frac{1}{6}+\frac{1}{12}erf\left(\frac{\varsigma -0.0335}{0.6890}\right)+\frac{1}{12}erf\left(\frac{\varsigma -0.7099}{0.8213}\right).$$

35

We conclude this section with a simulation of (ANF) for a particular parameter set (Table 1), which illustrates that our augmented neural field can generate interesting spatiotemporal behavior that closely resembles spiking patterns of a network of Rulkov neurons (SNM) with corresponding parameter values (Fig. 8). In the Rulkov network, synaptic input to neuron *i* is given by

$${u}_{n+1}^{(i)}=(1-{\alpha}_{i}){u}_{n}^{(i)}+{\alpha}_{i}\sum _{j=1}^{N}{c}_{ij}{s}_{n}^{(j)},$$

36

where *N* denotes the total number of neurons in the network, and *c*_{ij} is the connection strength from neuron *j* to neuron *i*. To match the parameters in Table 1, we split the total population in two subpopulations of 300 neurons each, which are both equidistantly placed on the interval [−1, 1]. Neurons within the same subpopulation share the same intrinsic parameters, and uncorrelated (in space and time) Gaussian noise is added to the threshold parameters. Finally, the connection strengths in the Rulkov network are given by

37

where *p*_{i} and *x*_{i} are the subpopulation and position of neuron *i*, respectively.

This paper presents a simple rate-reduced neuron model that is based on a variation of the Rulkov map (Rulkov [14], Rulkov et al. [15]), and can be used to incorporate a variety of non-trivial spiking behavior into existing neural field models.

The modified Rulkov map (SNM) is a phenomenological, two-dimensional single neuron model. The isolated dynamics of its fast time scale either generates a stable limit cycle, mimicking spiking activity, or a stable fixed point, corresponding to a neuron at rest (Fig. 1). The slow time scale of the Rulkov map acts as a dynamic spiking threshold and emulates the combined effect of slow recovery processes. The modified Rulkov map can mimic a wide variety of spiking patterns, such as spike-frequency adaptation, postinhibitory rebound, phasic spiking, spike accommodation, spike latency and inhibition-induced spiking (Fig. 2). Furthermore, the model can be used to model neuronal filter properties. Depending on how external input is applied to the model, it can act as either a high-pass or low-pass filter (Figs. 3 and and44).

The rate-reduced model (RNM) is derived heuristically and given by a simple one-dimensional differential equation. On the single cell level, the rate-reduced model closely mimics the spiking dynamics (Fig. 5) and filter properties (Fig. 6) of the full spiking neuron model. While a close approximation of the (expected) firing rate of Rulkov neurons (Fig. 7) is needed to mimic their behavior quantitatively, the types of qualitative dynamics of the rate-reduced model do not depend on the exact choice of firing rate function.

Due to its simplicity, it is straightforward to add the rate-reduced model to existing neural field models. In the resulting augmented equations, parameters can be chosen according to the spiking behavior of a single isolated cell. In our particular example (ANF), the emerging spatiotemporal pattern closely resembles the dynamics of the corresponding spiking neural network (Fig. 8). We believe that this is an elegant way to add more biological realism to existing neural field models, while simultaneously enriching their dynamical structure.

We used a variation of a simple toy model of a spiking neuron (Rulkov [14], Rulkov et al. [15]) to derive a corresponding rate-reduced model. While being purely phenomenological, the model could mimic a wide variety of biologically observed spiking behaviors, yielding a simple way to incorporate complex spiking behavior into existing neural field models. Since all parameters in the resulting augmented neural field equations have a representative in the spiking neuron network (and vice versa), this greatly simplifies the otherwise often problematic translation from results obtained by neural field models back to biophysical properties of spiking networks. An example demonstrated that the augmented neural field equations can produce spatiotemporal patterns that cannot be generated with corresponding ‘classical’ neural fields.

The conclusions of this paper are solely based on mathematical models.

Authors’ contributions

Conceptualization, K.D., Y.K., M.v.P. and S.v.G.; methodology, K.D. and S.v.G.; investigation, K.D.; writing original Draft, K.D.; writing review & Editing, K.D., Y.K., M.v.P. and S.v.G.; visualization, K.D.; supervision, Y.K., M.v.P. and S.v.G. All authors read and approved the final manuscript.

Funding

K.D. was supported by a grant from the Twente Graduate School (TGS).

Our study don’t involve human participants, human data or human tissue.

The authors declare no competing financial interests.

This manuscript does not contain any individual person’s data.

**Publisher’s Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Koen Dijkstra, Email: ln.etnewtu@artskjid.neok.

Yuri A. Kuznetsov, Email: ln.uu.htam@tenzuk.

Michel J. A. M. van Putten, Email: ln.etnewtu@nettupnav.m.a.j.m.

Stephan A. van Gils, Email: ln.etnewtu@slignav.a.s.

1. Beurle RL. Properties of a mass of cells capable of regenerating pulses. Philos Trans R Soc Lond B. 1956;240:55–94. doi: 10.1098/rstb.1956.0012. [Cross Ref]

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

3. Wilson HR, Cowan JD. A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik. 1973;13:55–80. doi: 10.1007/BF00288786. [PubMed] [Cross Ref]

4. Nunez PL. The brain wave equation: A model for the EEG. Math Biosci. 1974;21:279–297. doi: 10.1016/0025-5564(74)90020-0. [Cross Ref]

5. Amari S. Homogeneous nets of neuron-like elements. Biol Cybern. 1975;17:211–220. doi: 10.1007/BF00339367. [PubMed] [Cross Ref]

6. Amari S. Dynamics of pattern formation in lateral-inhibition type neural fields. Biol Cybern. 1977;27:77–87. doi: 10.1007/BF00337259. [PubMed] [Cross Ref]

7. Pinto DJ, Ermentrout GB. Spatially structured activity in synaptically coupled neuronal networks: I. Traveling fronts and pulses. SIAM J Appl Math. 2001;62:206–225. doi: 10.1137/S0036139900346453. [Cross Ref]

8. Pinto DJ, Ermentrout GB. Spatially structured activity in synaptically coupled neuronal networks: II. Lateral inhibition and standing pulses. SIAM J Appl Math. 2001;62:226–243. doi: 10.1137/S0036139900346465. [Cross Ref]

9. Coombes S, Owen MR. Bumps, breathers, and waves in a neural network with spike frequency adaptation. Phys Rev Lett. 2005;94:148102. doi: 10.1103/PhysRevLett.94.148102. [PubMed] [Cross Ref]

10. Kilpatrick ZP, Bressloff PC. Effects of synaptic depression and adaptation on spatiotemporal dynamics of an excitatory neuronal network. Physica D. 2010;239:547–560. doi: 10.1016/j.physd.2009.06.003. [Cross Ref]

11. Kilpatrick ZP, Bressloff PC. Stability of bumps in piecewise smooth neural fields with nonlinear adaptation. Physica D. 2010;239:1048–1060. doi: 10.1016/j.physd.2010.02.016. [Cross Ref]

12. Nicola W, Campbell SA. Bifurcations of large networks of two-dimensional integrate and fire neurons. J Comput Neurosci. 2013;35:87–108. doi: 10.1007/s10827-013-0442-z. [PubMed] [Cross Ref]

13. Visser S, van Gils SA. Lumping Izhikevich neurons. EPJ Nonlinear Biomed Phys. 2014;2:226–243. doi: 10.1140/epjnbp19. [Cross Ref]

14. Rulkov NF. Modeling of spiking-bursting neural behavior using two-dimensional map. Phys Rev E. 2002;65:041922. doi: 10.1103/PhysRevE.65.041922. [PubMed] [Cross Ref]

15. Rulkov NF, Tomofeev I, Bazhenov M. Oscillations in large-scale cortical networks: Map-based model. J Comput Neurosci. 2004;17:203–223. doi: 10.1023/B:JCNS.0000037683.55688.7e. [PubMed] [Cross Ref]

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

17. Izhikevich EM. Which model to use for cortical spiking neurons? IEEE Trans Neural Netw. 2004;15:1063–1070. doi: 10.1109/TNN.2004.832719. [PubMed] [Cross Ref]

Articles from Journal of Mathematical Neuroscience are provided here courtesy of **Springer-Verlag**