PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of jmathneurSpringerOpen.comSubmit OnlineRegisterThis journalThis article
 
J Math Neurosci. 2017; 7: 13.
Published online 2017 December 11. doi:  10.1186/s13408-017-0055-3
PMCID: PMC5725415

A Rate-Reduced Neuron Model for Complex Spiking Behavior

Abstract

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.

Introduction

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.

Single Spiking Neuron Model

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

{vn+1=f(vn,vn1,κunanθ),an+1=anε(an+(1κ)unγsn),
SNM

where the piecewise continuous function f:ℝ3 → ℝ is given by

f(x1,x2,x3)={2500+150x150x1+50x3if x1<0,50+50x3if 0x1<50+50x3x2<0,50otherwise.
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

sn={1if the neuron spiked at iteration n,0otherwise.
2

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

vn ≥ 0  ∧  (vn ≥ 50 + 50(κun − an − θ)  ∨  vn−1 ≥ 0).
3

The dependence of vn+1 on vn−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 unφ, the neuron spikes persistently if and only if φ > θ. After a change of variable an → an + (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 κun − an, 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.

Fast Dynamics

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 κun − an − θς is constant. In this case, (SNM) reduces to the fast subsystem

vn+1={2500+150vn50vn+50ςif vn<0,50+50ςif 0vn<50+50ς,50otherwise.
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

vs=25(ς2ς28ς),vu=25(ς2+ς28ς),
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 vu 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(ς)={0for ς0,1P(ς)for ς>0,
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]).

Fig. 1
Illustration of the fast subsystem (FSS) of (SNM). ( A ) For ς=110 there exist a stable (green) and unstable (orange) fixed point. ( B ) For ς=110 the system will settle into a stable periodic orbit (dashed green ...

Spiking Patterns

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.

Tonic Spiking/Fast Spiking

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>ε>110) 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 κ.

Fig. 2
Different types of spiking patterns generated by the single neuron model (SNM). Corresponding parameter values (θκεγ) are given in brackets. ( A ) Tonic spiking (110,12,12,12). ( B ) Spike-frequency adaptation (110,1,11000,5). ( ...

Spike-Frequency Adaptation/Regular Spiking

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

Rebound Spiking and Accommodation

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’ (un < θ) input can elicit temporary spiking if the input is ramped up sufficiently fast (Fig. 2D).

Spike Latency and Inhibition-Induced Spiking

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]).

Neuronal Filtering

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

un=φcos(ωπn1000+ϑ),
6

with amplitude [var phi], 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

κun − an ≤ θn.
7

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

an+1 = (1 − ε)an − ε(1 − κ)un
8

with explicit solution

an=ε(1κ)m=1(1ε)m1unm.
9

Inserting (6) into (9) now yields

κunan=κφcos(ωπn1000+ϑ)+ε(1κ)φm=1(1ε)m1cos(ωπ(nm)1000+ϑ)=F(ω)φ2e(ωπn1000+ϑ)i+F(ω)φ2e(ωπn1000+ϑ)i=|F(ω)|φcos(ωπn1000+ϑ+arg(F(ω))),
10

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

F(ω)=κ+ε(1κ)eωπi1000+ε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 [theta] in (6) are chosen such that

cos(ωπn1000+ϑ+arg(F(ω)))=1for some nN.
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(ω)|=ε2+2κ(κε)(1cos(ωπ1000))ε2+2(1ε)(1cos(ωπ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,F(1000)=2κε2ε.
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).

Fig. 3
Illustration of the frequency response (11) for different values of ε. ( A ) For κ=110 high frequencies get attenuated. ( B ) For κ = 2 high frequencies get amplified. Note the similarity, which is caused by the fact that ...
Fig. 4
Responses of (SNM) to periodic input, illustrating neuronal filter properties. ( A ) For κ=110 the neuron acts as a low-pass filter. Input with an amplitude of φ = ⅕ elicits a spiking response for ω = 1, whereas the neuron is quiescent ...

The Rate-Reduced Neuron Model

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

an+1an − ε(an + (1 − κ)un − γS(κun − an − θ)).
16

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

1εdadt=a(1κ)u+γS(κuaθ).
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).

Fig. 5
Different types of spiking behavior generated by the rate-reduced model (RNM). Top traces show the firing rate with r(t) = κu(t)−a(t)−θ. Corresponding parameter values (θκεγ) are given in brackets. For small values of ε (i.e. a large time ...

Frequency Response of the Reduced Model

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

u(t)=φcos(ωπt1000+ϑ).
17

Under the assumption that

κu(t)−a(t) ≤ θt
18

the explicit solution of (RNM) is given by

a(t)=ε(1κ)teε(tτ)u(τ)dτ,
19

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

κu(t)a(t)=κφcos(ωπt1000+ϑ)+ε(1κ)φteε(tτ)cos(ωπτ1000+ϑ)dτ=G(ω)φ2e(ωπt1000+ϑ)i+G(ω)φ2e(ωπt1000+ϑ)i=|G(ω)|φcos(ωπt1000+ϑ+arg(G(ω))),
20

where the frequency response G:ℝ≥0 ↦ ℂ is given by

G(ω)=κ+ε(1κ)ε+ωπi1000.
21

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

S(κu(t)−a(t)−θ) = 0 ∀t if and only if |G(ω)|φ ≤ θ.
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(ω)|=ε2+κ2(πω1000)2ε2+(πω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 (κ,ε,θ,φ)=(110,1200,17,15), and

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

for (κ,ε,θ,φ)=(2,1200,17,110). 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.

Fig. 6
Responses of the rate-reduced model (RNM) to periodic input. Top traces show the firing rate with r(t) = κu(t)−a(t)−θ. ( A ) For κ=110 the model acts as a low-pass filter. Input with an amplitude of φ = ⅕ yields a response in the firing ...

The Firing Rate Function

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 ς1>ςς2=12(517). In general, at ςk, there is a jump discontinuity of size

limςςk+S(ς)limςςkS(ς)=1(k+2)(k+3),with S(ςk)=1k+2.
26

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

S(ς)=k=1H(ςςk)(k+2)(k+3),
27

where H is the Heaviside step function and

limkς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

S(ς)=12πσ2ew22σ2S(ς+w)dw,
29

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

S(ς)=16+k=1erf(ςςk2σ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

16+16Ni=1Nerf(ςνiχi),
31

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

12πσ2ew22σ2S(ς+w)dw1616Ni=1Nerf(ςνiχi)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).

Fig. 7
Expected firing rate for a noise level of σ2 = ¼. Shown are a numerical integration of (29) (blue) and its approximation (31) for N = 2 and (ν1ν2χ1χ2) = (0.0335, 0.7099, 0.6890, 0.8213) (orange)

Augmenting Neural Fields

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):

(1+1α1t)u1(t,x)=11J11(x,x)S1(r1(t,x))+J12(x,x)S2(r2(t,x))dx,(1+1ε1t)a1(t,x)=(1κ1)u1(t,x)+γ1S1(r1(t,x)),(1+1α2t)u2(t,x)=11J21(x,x)S1(r1(t,x))+J22(x,x)S2(r2(t,x))dx,(1+1ε2t)a2(t,x)=(1κ2)u2(t,x)+γ2S2(r2(t,x)),
ANF

where, as before,

ri(xt) = κiui(tx)−ai(tx)−θi
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 Jij(xx) measures the connection strength from neurons of population j and position x to neurons of population i and position x. The connectivity kernels Jij:Ω×ΩR are assumed to be isotropic and given by

Jij(xx) = ρjηijeμij|xx|
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 Si:ℝ ↦ ℝ are chosen to approximate the expected firing rate of Rulkov neurons (SNM) with a noise level of σ2 = ¼ (Fig. 7),

S1(ς)=S2(ς)=16+112erf(ς0.03350.6890)+112erf(ς0.70990.8213).
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

un+1(i)=(1αi)un(i)+αij=1Ncijsn(j),
36

where N denotes the total number of neurons in the network, and cij 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

cijηpipjeμpipj|xixj|
37

where pi and xi are the subpopulation and position of neuron i, respectively.

Fig. 8
Spatio-temporal spiking patterns. ( A ) Simulation of the augmented neural field (ANF) with parameter values given in Table Table1.1. Shown is the firing rate S1(κ1u1(tx)−a1(tx)−θ1) of the first population. ( B ) Simulation ...
Table 1
Parameter overview for the neural field (ANF)

Discussion

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.

Conclusions

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.

Availability of data and materials

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

Authors’ contributions

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

Funding

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

Notes

Ethics approval and consent to participate

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

Competing interests

The authors declare no competing financial interests.

Consent for publication

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

Footnotes

Publisher’s Note

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

Contributor Information

Koen Dijkstra, ln.etnewtu@artskjid.neok.

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

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

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

References

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]
18. Chik DTW, Coombes S, Wang ZD. Clustering through postinhibitory rebound in synaptically coupled neurons. Phys Rev E. 2004;70:011908. doi: 10.1103/PhysRevE.70.011908. [PubMed] [Cross Ref]
19. Blumhagen F, Zhu P, Shum J, Schärer Y-PZ, Yaksi E, Deisseroth K, Friedrich RW. Neuronal filtering of multiplexed odour representations. Nature. 2011;479:493–498. doi: 10.1038/nature10633. [PubMed] [Cross Ref]

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