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 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 where the piecewise continuous function is given by 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 A Rulkov neuron spikes at iteration n if its membrane potential is reset to in the next iteration. It follows from (1) that the spiking condition in (2) is satisfied if and only if The dependence of on 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 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 , the neuron spikes persistently if and only if . After a change of variable 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 , κ 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 , changes in external input are initially magnified for . Asymptotically, this is then counterbalanced by additional adaption. Finally, for , 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 is a slow-fast system, and we can explore the fast spiking dynamics of the model by assuming the suprathreshold drive is constant. In this case, (SNM) reduces to the fast subsystem The map (FSS) undergoes a saddle-node bifurcation at (Fig. 1). For 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 (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 always lies to the right of the ‘reset potential’ , a stable fixed point and a periodic orbit can never coexist. This guarantees that we can define a firing rate function for the fast subsystem (FSS), given by where 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 , the ‘middle’ interval , and the ‘right’ interval . 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 there exist a stable (green) and unstable (orange) fixed point. ( B ) For the system will settle into a stable periodic orbit (dashed green line) with period

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’ 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 . ( B ) Spike-frequency adaptation . ( C ) Rebound spiking . ( D ) Accommodation . ( E ) Spike latency . ( F ) Inhibition-induced spiking

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 () and choosing . 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 . 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’ () 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 (), 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 , 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 , and where 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 In this case, the adaptation reduces to the simple linear equation with explicit solution Inserting (6) into (9) now yields where the overline denotes complex conjugation and the frequency response 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 is strictly decreasing if and only if , and increasing otherwise (Fig. 3). Clearly, The input parameter κ can therefore be used to model filter properties of the neuron. For 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 (and ), 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 high frequencies get attenuated. ( B ) For high frequencies get amplified. Note the similarity, which is caused by the fact that is an odd function of

Fig. 4

Responses of (SNM) to periodic input, illustrating neuronal filter properties. ( A ) For the neuron acts as a low-pass filter. Input with an amplitude of elicits a spiking response for , whereas the neuron is quiescent for . ( B ) For , the neuron acts as a high-pass filter. Input with amplitude elicits a spiking response for , whereas a lower input frequency of does not. In both examples,

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

Fig. 5

Different types of spiking behavior generated by the rate-reduced model (RNM). Top traces show the firing rate with . Corresponding parameter values are given in brackets. For small values of ε (i.e. a large time scale separation), there is excellent agreement with the corresponding examples of the full model (Fig. 2), which is quantified by comparing the integral of the spiking rate in the reduced model to the number of spikes in the full model. ( A ) Tonic spiking ; . ( B ) Spike-frequency adaptation ; . ( C ) Rebound spiking ; . ( D ) Accommodation ; . ( E ) Spike latency ; . ( F ) Inhibition-induced spiking ;

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 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 and therefore is strictly decreasing if and only if , and increasing otherwise. Revisiting the examples from Sect. 2.3 (Fig. 4), we have for , and for . 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 . ( A ) For the model acts as a low-pass filter. Input with an amplitude of yields a response in the firing rate for , whereas the firing rate remains zero for . 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 , the reduced model acts as a high-pass filter. Input with amplitude elicits a firing rate response for , whereas a lower input frequency of 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). In both examples,

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 . Let denote the set of all points of discontinuity of the firing rate function in decreasing order. For the ‘reset potential’ 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 . In general, at , 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

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 , given by where is the variance of the noise. Using (27), we can rewrite (29) as where erf denotes the error function. While 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 by a finite sum of the form for some fixed and constants , which are chosen by (numerically) minimizing For large noise levels , 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 . Shown are a numerical integration of (29) (blue) and its approximation (31) for and (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 : where, as before, for . The differential operators in the left-hand side of the integral equations in (ANF) model exponentially decaying synaptic currents with decay rate . The connectivity measures the connection strength from neurons of population j and position to neurons of population i and position x. The connectivity kernels are assumed to be isotropic and given by where is the density of neurons of type j, is the maximal connection strength, and is the spatial decay rate of the connectivity. Both firing rate functions are chosen to approximate the expected firing rate of Rulkov neurons (SNM) with a noise level of (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 where N denotes the total number of neurons in the network, and 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 . 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 and are the subpopulation and position of neuron i, respectively.

Table 1

Parameter overview for the neural field ()

Parameter	α_i	θ_i	κ_i	ε_i	γ_i	ρ_i	η_i1	η_i2	μ_i1	μ_i2
Population 1	120	12	2	11000	5	150	23	−13	4	1
Population 2	110	45	110	1100	2	150	1115	−1130	174	1110

Fig. 8

Spatio-temporal spiking patterns. ( A ) Simulation of the augmented neural field (ANF) with parameter values given in Table 1. Shown is the firing rate of the first population. ( B ) Simulation of a corresponding network of 300 excitatory and 300 inhibitory Rulkov neurons, all-to-all coupled via simple exponential synapses. Both populations are equidistantly placed on the interval . Uncorrelated (in space and time) Gaussian noise with variance is added to the threshold parameter of each neuron. Shown is the spiking activity of the excitatory population. Each spike is denoted by a black dot

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

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.