A Rate-Reduced Neuron Model for Complex Spiking Behavior

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 B K. Dijkstra koen.dijkstra@utwente.nl

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 v n+1 = f (v n , v n−1 , κu n − a n − θ), a n+1 = a n − ε(a n + (1 − κ)u n − γ s n ), where the piecewise continuous function f : R 3 → R is given by otherwise. (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 = 1 if the neuron spiked at iteration n, 0 otherwise.
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 v n ≥ 0 ∧ v n ≥ 50 + 50(κu n − a n − θ) ∨ v n−1 ≥ 0 .
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.

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 κu n − otherwise.

(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 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 : R → Q for the fast subsystem (FSS), given by where P : R >0 → N 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]).

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

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 'spikefrequency 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' (u n < θ) 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 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 κu n − a n ≤ θ ∀n.
In this case, the adaptation reduces to the simple linear equation with explicit solution Inserting (6) into (9) now yields κu n − a n = κϕ cos ωπn 1000 where the overline denotes complex conjugation and the frequency response F : [0, 1000] → C is given by 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 The inverse statement is not true, even if ω and ϑ in (6) are chosen such that 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 and it follows that |F | is strictly decreasing if and only if κ ∈ (−1 + ε, 1), and increasing otherwise (Fig. 3). Clearly, 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 lowpass 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).

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 a n+1 = a n − ε a n + (1 − κ)u n − γ S(κu n − a n − θ) .
By interpreting (16) as the forward discretization of an ordinary differential equation, we arrive at the continuous time rate-reduced model The rate-reduced neuron model (RNM) preserves the dynamical features of the full model (SNM) and reproduces all previous example spiking patterns (Fig. 5).

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 Under the assumption that the explicit solution of (RNM) is given by cf. (9). Inserting the input (17) into (19) yields where the frequency response G : R ≥0 → C is given by It follows that for the rate-reduced model (RNM) receiving harmonic input (17) we have 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  (16) 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 for (κ, ε, θ, ϕ) = ( 1 10 , 1 200 , 1 7 , 1 5 ), and for (κ, ε, θ, ϕ) = (2, 1 200 , 1 7 , 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. Input with an amplitude of ϕ = 1 5 yields a response in the firing rate for ω = 1, whereas the firing rate remains zero for ω = 2. In the former case, the integral of the spiking rate during one period is approximately 4.55, while there are 5 spikes per period in the full model (Fig. 4A). (B) For κ = 2, the reduced model acts as a high-pass filter. Input with amplitude ϕ = 1 10 elicits a firing rate response for ω = 2, whereas a lower input frequency of ω = 1 does not. In the former case, the integral of the spiking rate during one period is approximately 3.14, while there are 3 spikes per period in the full model (Fig. 4B)

The Firing Rate Function
Since our neuron model (SNM) is a map, the period P of its limit cycle lies in N for all positive suprathreshold drives ς . Therefore, the spiking rate function (5) is staircaselike, 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 nonnegative 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 = 1 2 (5 − √ 17). In general, at ς k , there is a jump discontinuity of size The firing rate of the fast subsystem (FSS) can therefore be written as where H is the Heaviside step function and lim k→∞ ς k = 0.
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 : R → R, given by where σ 2 is the variance of the noise. Using (27), we can rewrite (29) as where erf denotes the error function. While S(ς) can readily be computed for any ς ∈ R 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 Fig. 7 Expected firing rate for a noise level of σ 2 = 1 4 . 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) for some fixed N ∈ N and constants ν i , χ i ∈ R, which are chosen by (numerically) minimizing 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).

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): where, as before, 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 : Ω × Ω → R are assumed to be isotropic and given by 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 : R → R are chosen to approximate the expected firing rate of Rulkov neurons (SNM) with a noise level of σ 2 = 1 4 (Fig. 7), 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  Table 1. Shown is the firing rate  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 where p i and x i are the subpopulation and position of neuron i, respectively.

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 spikefrequency 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 4).
The rate-reduced model (RNM) is derived heuristically and given by a simple onedimensional 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.
Funding K.D. was supported by a grant from the Twente Graduate School (TGS).

Availability of data and materials
The conclusions of this paper are solely based on mathematical models.