Lumping Izhikevich neurons

Abstract

We present the construction of a planar vector field that yields the firing rate of a bursting Izhikevich neuron can be read out, while leaving the sub-threshold behavior intact. This planar vector field is used to derive lumped formulations of two complex heterogeneous networks of bursting Izhikevich neurons. In both cases, the lumped model is compared with the spiking network. There is excellent agreement in terms of duration and number of action potentials within the bursts, but there is a slight mismatch of the burst frequency. The lumped model accurately accounts for both intrinsic bursting and post inhibitory rebound potentials in the neuron model, features which are absent in prevalent neural mass models.

Background

The Izhikevich neuron is a computationally simple, yet biologically realistic neuron model [1, 2], which belongs to the class of nonlinear integrate-and-fire neurons. Bifurcation analysis of a general class of nonlinear integrate-and-fire neurons is given in [3, 4], whereas a more detailed analysis of the quadratic form is given in [5].

When coupled in a complex network, knowledge of the individual neurons does not suffice to predict the network behavior. As large networks are not suitable for bifurcation analysis, one seeks to formulate low dimensional reductions that capture the space-time-averaged behavior of the network. An example is the Wilson-Cowan reduction [6, 7] where the mean population activities are governed by a low dimensional system of ordinary differential equations (ODEs) or integro-differential equations respectively. Just as many other lumped formulations of neural activity, these models pivot on the relation between the individual cells’ firing rates and aggregate network activity [8–10].

Indeed, the most traditional form of a lumped model for a single population of neurons has the form

$$\mathit{\text{Qx}}=f(x,I)$$

(1)

with x representing either the mean membrane potential or the postsynaptic conductances/currents.

Q is a linear differential operator which characterizes the synaptic dynamics. Here, f depicts the firing rate as function of the population’s own activity x, due to connections within the population, and some external input I. This basic framework is easily extended to include multiple populations or other features, such as spike frequency adaptation and spatial structures (e.g. neural fields). Yet, it should be noted that models of this form only describe the dynamic evolution of populations averaged quantities.

Population density methods, on the other hand, describe the evolution of the distribution of neurons over the state space, which depends on the quantities that are used in the underlying neuron model [11–14]. In [14] application of the moment closure method is investigated, whose steady-state solution, in its simplest form (closure at the first moment), is the so-called mean-field firing-rate model. Next, in [15] this method is applied to arrive at a planar system of switching ODEs for all-to-all coupled networks of Izhikevich neurons with different topologies. This system is subsequently analyzed by means of numerical bifurcation methods [16] and direct simulation. Its variables are the mean adaptation and the synaptic input (here the adjective ‘mean’ can be omitted due to the all-to-all coupling). One assumption is that the mean adaptation follows the mean firing rate, rather than the individual firing rates. Furthermore a separation of time scales is assumed: the adaptation time scale is assumed to be much larger than the membrane time-constant.

The analogy with [15] is that we also consider networks of Izhikevich neurons. The novelty is the construction of the Izhikevich Planar Vector Field (IPVF) from which the firing rate function can be read out. The resulting reduced vector field is three dimensional in contrast to the two-dimensional vector field obtained by [15]; we keep the mean membrane potential as a state variables. This state variable follows naturally from the construction from the associated planar vector field. We claim that it is necessary to retain the membrane potential, particularly to deal with post-inhibitory rebound spikes/ bursts.

This paper is structured as follows: Izhikevich’ spiking neuron model is studied with parameters corresponding to an intrinsically bursting neuron. Using phase plane analysis, the discontinuous resets are approximated with a smooth branch, which relates to the firing rate of the cell. Most importantly, the proposed approximation leaves all of subthreshold dynamics intact. The resulting model is shown to match quantitatively with the original spiking model, with respect to both spike count and spike timing. Upon making several assumptions relating to network connectivity and synchrony, a random network of Izhikevich neurons can be collapsed into a system of ODEs. Finally, we give two examples of random networks, not necessarily all-to-all connected, of heterogeneous Izhikevich neurons where we compare the true network behavior with the properties of the corresponding reduced vector field.

Methods

Izhikevich neuron in the phase plane

The Izhikevich neuron is given by the following equations [1]

$$\left\{\begin{array}{l}\dot{V}=0.04V^2+5V+140-u+I_{\text{ext}}\left(t\right)\\ \dot{u}=a(\mathit{\text{bV}}-u)\hfill \end{array}\right.,V\left(t^{-}\right)=30\Rightarrow \left\{\begin{array}{c}V\left(t^{+}\right)=c\\ u\left(t^{+}\right)=u\left(t^{-}\right)+d\end{array}\right.$$

(2)

Indeed, with four parameters the model is able to reproduce many of the key features observed in real neural tissue, e.g. intrinsic spiking/bursting, spike-frequency adaptation, and rebound spikes/bursts. Due to a reset process, present in all other integrate-and-fire models, the model does not require a fast recovery variable relating to a resonating current. Hence, the variable u is often seen as a slow variable within the neuron, e.g. Calcium concentration.

Since the Izhikevich model has two variables, phase plane analysis is a natural starting point for characterizing the model. Figure 1 depicts the phase portrait of an intrinsically bursting neuron: Starting on the V-nullcline at point A, the dynamics are exclusively governed by the u component which are, in this case, directed downwards. As u decreases slowly, the fast dynamics in the V direction push the system towards the V-nullcline, causing the orbit to closely follow the quadratic V-nullcline. After passing point B, V increases fast until the threshold V=30 is reached and V is reset onto the green line with increased u. As long as the system’s reset point is right of the unstable branch of the V-nullcline, another spike is generated. Note that this condition is equivalent to the resets points lying below point C: the intersection point of the V-nullcline and the reset line V=c (in green). Upon being reset above C, the dynamics in the V-direction are pointing towards the left and the system evolves in that direction until the limit cycle is completed when hitting A.

Figure 1

V − u phase plane of the Izhikevich neuron. (Blue) nullclines of V and u, (green) reset value c, and (red) periodic solution of the system. The red dots correspond with resets. See text for description of the points A, B, and C.

The main objective is to replace the fast spiking dynamics by a slowly changing variable which represents the neurons firing rate, whilst leaving the subthreshold dynamics intact. In particular:

• The right branch of the V-nullcline corresponds with the firing threshold.

• The left branch of the V-nullcline characterizes the fast dynamics in the subthreshold regime.

• The u-nullcline left of the firing threshold determines the slow dynamics in the subthreshold regime.

• Point B corresponds to the entrance of the bursting regime.

• Point C corresponds to the exit of the bursting regime.

• Point A corresponds to the undershoot after burst termination.

An interpretation of these criteria, expressed as a set of ODEs, is as follows.

Reduction

First, the model’s ability to produce spikes, i.e. V quickly reaches 30mV, is taken away by adding another branch to the V-nullcline, see Figure 2. Orbits close to this additional branch correspond with a neuron being in burst-state and, since the condition (non-) bursting is binary, the precise shape and position of this artificial branch are of little importance.

Figure 2

V − u phase plane of the reduced model. (Blue) nullclines of V and u, (green) reset value c, and (red) periodic solution of the system. See text for description of the points A, B, C, and C^′.

In the bursting regime, the time between consecutive spikes is approximated as follows. Assuming that J = − u + I_ext (t) varies sufficiently slow to be considered as a constant parameter for V, the reduced membrane potential v satisfies:

$$\dot{v}=0.04v^2+5v+140+J,v\left(0\right)=c$$

(3)

where the fact is used that v is reset to c after each spike. Albeit nonlinear, this equation has an analytic solution v (t;J,c) which is used to find the time-to-spike T _s for which v (T _s ;J,c) = 30:

$$T_s(J;c)=T(J;30)-T(J;c),$$

with function T given by:

$$T(J;\nu )=\frac{10}{\sqrt{4J-65}}\text{arctan}\left(\frac{2\nu +125}{5\sqrt{4J-65}}\right)$$

In order to have T _s (J;c) positive, it is required that (3) has a positive right hand side at the initial value:

$$0.04c^2+5c+140+J>0.$$

(4)

Inside the bursting regime the instantaneous firing rate is now readily determined as the reciprocal of the inter spike times.

For reasons explained later, the firing rate inside the bursting regime is capped at a minimal value:

$$f_b(J;c)=\left\{\begin{array}{l}\text{max}\left(\frac{1}{T_s(J;c)},f_{\text{min}}\right)0.04c^2+5c+140+J>0,\\ f_{\text{min}}\text{otherwise}\hfill \end{array}\right.$$

In order to distinguish between bursting or non-bursting regime, the actual firing rate (where J is substituted) is set to

$$f(V,u,I_{\text{ext}};c)=H_{\omega }(V-c)f_b(-u+I_{\text{ext}};c)$$

(5)

with H _ω a smooth approximation to the Heaviside step function with steepness-parameter ω. The step function aligned at the reset value now distinguishes between the bursting and non-bursting regimes of the neuron. In the remainder of this text, the following is used:

$$H_{\omega }\left(V\right)=\frac{1}{1+\text{exp}(-\omega V)}$$

with ω = 20.

This firing rate is not only the main output variable of the neuron, for it also finds its way back into in the dynamical system. Namely, every time a spike is generated, the recovery variable is increased by d, which yields an time-averaged increase rate of fd in the reduced model. Hence, the recovery variable in the reduced model takes the form:

$${\dot{u}}_r=a({\mathit{\text{bV}}}_r-u_r)+d·f(V_r,u_r,I_{\text{ext}};c).$$

As a final step, the V-nullcline has to be adjusted such that another branch appears in the phase portrait, c.f. Figure 2. The most natural choice would be to have this branch follow the time-averaged value of V during a spike, i.e.

$$\bar{v}=\frac{1}{T_s(J;c)}{\int }_0^{T_s(J;c)}v\left(t\right)\mathit{\text{dt}}$$

Although a closed expression can be found for this integral, this approach is not pursued here. The primary reason is that the precise position and shape of the additional branch are merely insignificant. Instead, a much simpler expression is chosen which does not obscure the resulting equations as much as the expression for $\bar{v}$ would.

Despite the fact that the shape of the new branch is of little importance, the local maximum at C^′ is not. Since at this point the stability of the slow manifold, i.e. V-nullcline, changes due to a fold bifurcation, it corresponds with the termination of the burst. For that reason it is important that this maximum is in line with the original termination point C, such that the undershoot at A is preserved accurately. Taking into account the requirement that subthreshold behavior should not be modified too much, the new branch is conveniently introduced by an exponential. In this case, the full reduced system takes the form:

$$\begin{array}{ll}{\dot{V}}_r& =g(V_r,u_r,I_{\text{ext}}):=0.04V_r^2+5V_r+140-u_r+I_{\text{ext}}-\xi \text{exp}(\eta (V_r-c\left)\right)\\ {\dot{u}}_r& =h(V_r,u_r,I_{\text{ext}}):=a({\mathit{\text{bV}}}_r-u_r)+d·f(V_r,u_r,I_{\text{ext}};c)\end{array}$$

(6)

with η and ξ parameters characterizing the additional branch. Appropriate values for these parameters can be determined as follows: Have u=n _v (V) depict the V-nullcline of the original Izhikevich model, then the level of point C in Figure 2 is given by n _v (c). For u = n _vr (V _r ) = n _v (V _r ) − ξ exp (η (V _r − c)) the V _r -nullcline of the reduced model, the following equations should hold to have point C^′ at the same level as C :

$$\begin{array}{ll}{n}_{\mathit{\text{vr}}}(c+\kappa )& =n_v\left(c\right)\\ n_{\mathit{\text{vr}}}^{\prime }(c+\kappa )& =0\end{array}$$

with κ representing the horizontal separation between C and C^′, and n vr ′ denoting the derivative of n _vr . Solving the equations results in:

$$\eta =\frac{n_v^{\prime }(c+\kappa )}{n_v(c+\kappa )-n_v\left(c\right)}\xi =\frac{n_v(c+\kappa )-n_v\left(c\right)}{e^{\eta \kappa }}.$$

where $n_v^{\prime }$ is the derivative of n _v . It is undesirable to choose low values of κ since the additional branch will be very steep and the stiffness of the dynamical system increases unnecessarily. Large values, on the contrary, will result in a branch which is too gradual, such that both the timing and the subthreshold behavior are affected. Values of κ in the interval [ 0.5,2] give merely identical results; here is chosen for κ=0.8.

Remark

Since the appearance f_min in the reduction is unexpected, it is worthwhile to discuss both the necessity and interpretation for this parameter. In the reduced vector field, an equilibrium on the right branch of the V _r -nullcline corresponds to a tonically firing neuron, while absence of such an equilibrium corresponds to an intrinsically bursting neuron. In the case of the bursting neuron, passage from the right branch of the V _r -nullcline to the left is guaruanteed by raising the u _r -nullcline such that possible equilibria are removed. This raise is achieved by the introduction of f_min.

Physiologically, f_min represents the minimum the firing rate for a neuron at which it can still be considered ‘bursting’. While this depends on the type of cell at hand, it is generally expected that large cells will have lower firing rates (also within bursts) than small cells.

Single cell comparison

In order to assess the quality of the proposed reduction (6), simulations are performed for different parameters of the model. Not only the membrane potentials and recovery variables of both models, also the output firing rate of the reduced model is matched with the discrete spikes of the original model. For that, define:

$$\begin{array}{ll}N\left(t\right)& :=\text{number of spikes in}[0,t]\end{array}$$

(7a)

$$\begin{array}{ll}{N}_r\left(t\right)& :={\int }_0^{t}f\left(V_r\right(\tau ),u_r(\tau ),I_{\text{ext}}(\tau \left)\right)d\tau \end{array}$$

(7b)

for the Izhikevich and the reduced model respectively.

Figure 3 shows a typical comparison of a tonically bursting neuron. The general impression is that the reduction is very successful in mimicking the original Izhikevich model. Although in the latter half of the simulation the inter-burst frequency of the reduced model is higher, which results in a phase difference between both models, the average spiking behavior, on the other hand, is accurately preserved by the reduction.

Figure 3

Tonic bursting neuron, which is suppressed with a negative current. Top, middle, and bottom diagrams correspond with V, u and N (c.f. (7)) respectively. Parameters as follows: a = 0.02, b = 0.2, c = − 55, d = 2, κ = 0.8, and f_min = 0.1. Furthermore, a background current I₀ is given of magnitude 6, which is set to −6 between 200–500 ms, resulting in a suppression of the bursting pattern.

A comparison for the generation of post-inhibotory rebound potentials of a single cell in given in Appendix A: Single cell comparison with post-inhibitory rebound potentials.

A population of Izhikevich neurons

With the reduction of a single neuron into place, the next natural step is to consider the global activity patterns of a network of neurons. A single population of merely identical cells is considered first, which can be generalized to multiple populations afterwards.

For now, consider a network of N neurons, which are randomly connected. It is assumed that each spike arriving at a synaptic terminal at time $\widehat{t}$ results in an identical response of the synapse, say $q(t-\widehat{t})$, yielding the postsynaptic current

$$I_{\text{post}}\left(t\right)=\bar{g}q(t-\widehat{t})\left(V_{\text{post}}\right(t)-E_{\text{syn}})$$

with $\bar{g}$ representing the maximal conductance, E_syn the associated reversal potential and V_post the membrane potential of the postsynaptic neuron. Furthermore, let q be the impulse response of a linear operator Q such that

$$\mathit{\text{Qq}}=\delta$$

where δ is the Dirc delta mass centered at 0 and equality holds in the sense of distributions. For any neuron i, let $T_i\left(t\right):=\{{\widehat{t}}_j,j=1,\dots ,n|0<{\widehat{t}}_1<\dots <{\widehat{t}}_n\le t\}$ be the (possibly empty) set of spike times of neuron i up to and including time t. Define the synaptic activation s _i of postsynaptic channels as the convolution of q with the spike times T _i , then the following two formulations are equivalent

$$s_i\left(t\right):=\sum _{\widehat{t}\in T_i\left(t\right)}q(t-\widehat{t})\iff {\mathit{\text{Qs}}}_i\left(t\right)=\sum _{\widehat{t}\in T_i\left(t\right)}\delta (t-\widehat{t})$$

Next, have w _ij represent the synaptic weight of the connection from neuron j onto neuron i and note that w _ij = 0 implies that neuron j does not connect with i. The postsynaptic current arriving at cell i is then given by

$$I_{\text{syn},i}\left(t\right)=\bar{g}\sum _{j=1}^N{w}_{\mathit{\text{ij}}}{s}_j\left(t\right)\left(V_i\right(t)-E_{\text{syn}})$$

(8)

with V _i the membrane potential of neuron i.

In order to deduce the average network activity, a number of assumptions has to be made. Since all of these assumptions are common practice in neural mass modeling (c.f. [17]), we only state them briefly in the following part. The reader, however, is encouraged to challenge these assumptions and reflect on their implications.

• Temporal averaging. First the individual spike trains are averaged by representing them with the firing rate. That is, every neuron is described by the reduced system (6) as described in the previous section. The corresponding synaptic variable satisfies:

  $${\mathit{\text{Qs}}}_i=f(V_i,u_i,I_{\text{ext}}-I_{\text{syn},i})$$

• Randomly connected network. Let W be a random variable representing the synaptic weight of connections, such that all synaptic weights w _ij are assumed to be independently identically distributed as W. Note that if P (W = 0) > 0, the network is not necessarily all-to-all connected. Define

  $$\bar{w}:=\mathbb{EW},S\left(t\right):=\sum _{i=1}^N{s}_i\left(t\right)$$

which gives rise to the expected postsynaptic current at cell i

$$I_{\text{syn},i}\left(t\right)=\bar{g}\bar{w}S\left(t\right)\left(V_i\right(t)-E_{\text{syn}})$$

(9)

• Synchrony. Neurons receiving a similar input S are assumed to be reside closely together in phase space. Hence, the averages

  $$V\left(t\right):=\frac{1}{N}\sum _{i=1}^N{V}_i\left(t\right),u\left(t\right):=\frac{1}{N}\sum _{i=1}^N{u}_i\left(t\right)$$

  (10)

• suffice to approximate, c.f. (6)

  $$\begin{array}{ll}\dot{V}& =\frac{1}{N}\sum _{i=1}^{N}g(V_i,u_i,I_{\text{ext}}-I_{\text{syn},i})\approx g(V,u,I_{\text{ext}}-I_{\text{syn}})\end{array}$$

  (11a)
  $$\begin{array}{ll}\dot{u}& =\frac{1}{N}\sum _{i=1}^{N}h(V_i,u_i,I_{\text{ext}}-I_{\text{syn},i})\approx h(V,u,I_{\text{ext}}-I_{\text{syn}})\end{array}$$

  (11b)

with

$$I_{\text{syn}}=I_{\text{syn}}\left(t\right)=\bar{g}\bar{w}S\left(t\right)\left(V\right(t)-E_{\text{syn}})$$

Finally, due to both synchrony and linearity of Q, one finds

$$\mathit{\text{QS}}=\mathit{\text{Nf}}(V,u,I_{\text{ext}}-I_{\text{syn}})$$

(11c)

The reduced model (11) now represents the lumped network activity.

Extensions to multiple populations are easily incorporated in a similar fashion. An inhibitory population, for instance, can be included by introducing another differential operator, say Q₂, which corresponds with a slower postsynaptic current, combined with an appropriate reversal potential.

Results

Comparison

Now that the general derivation of a neural mass model based on Izhikevich neurons is discussed, the correspondence with the network of spiking neurons they represent is studied next. Therefore, simulations of both the neural mass and the detailed model are performed and the relation between them is investigated.

First, only the dynamics of a single population of merely identical excitatory neurons are considered, since we are primarily interested in the quantitative and qualitative comparison of both formulations. A model with multiple populations is studied afterwards.

Setup

To start, exponential synapses with a decay rate α are chosen for all connections in all network, hence

$$Q:=\frac{d}{\mathit{\text{dt}}}+\alpha$$

Furthermore, in the spiking neuron network we set p_con as the probability for one neuron to make a connection to another. The synaptic weights of these connections will all be equal and set to w₀. The expected weight for a connection, is now given by:

$$\bar{w}=w_0{p}_{\text{con}}$$

(12)

The random generation of connections between cells will lead to slight inhomogeneities in the network. Additional inhomogeneities are introduced by randomizing the parameters of the individual cells. This type of perturbation will ensure that the synchronization patterns we observe are due to the network interactions, rather than the intrinsic properties of the neurons. Within the reduced model, the expected value of each of these parameters is used.

Simulations

Simulations of the spiking neuron network are performed with Norns — Neural Network Studio^a: a dedicated C++ simulation tool with an intuitive Matlab interface. Networks are generated with their connections as described above. The spike trains obtained from Norns are post-processed within Matlab, which we discuss in detail below.

The formulated neural mass model can, depending on the parameters, become a stiff system of ODEs. In that case, Matlab’s ode23s is used to for numerical integration, while ode45 is used otherwise.

Comparison

To assess the similarity between the models, it is desirable to compare the same physical quantities within both models. In this light, the synaptic activation S appears to be the best choice: it can be determined from individual spike trains and it is the key variable in the lumped model.

For the lumped model, S is readily available from the simulations. In the network of spiking neurons, we reconstruct the synaptic activation of all neurons by convolving their spike trains with the shape of a postsynaptic response.

One excitatory population

Simulations of a population of excitatory intrinsically firing neurons are shown in Figure 4. The top diagram shows the behavior of a single isolated neuron in the the network and its corresponding firing rate reduction, c.f. Figure 3. The raster plot of a 1000 randomly connected such neurons is shown in the middle and the total activity of both the network and the lumped model are shown in the bottom diagram. It can be seen that, although all connections are excitatory, the frequency of bursting of the network is lower than that of the single cell.

Figure 4

Simulations for an isolated cell and a 1000 cell network of excitatory intrinsically bursting neurons. Top: membrane potential of isolated Izhikevich neuron with average parameters showing the original spiking neuron (blue) and the approximation by the reduced model (red). A random network of 1000 such spiking neurons is created, of which spike trains of 50 neurons are depicted in the raster plot (middle). Bottom diagram depicts the convoluted spike trains of all neurons in the network (blue), c.f. Section ‘Comparison’. The synaptic variable S of the corresponding lumped model of this network is also shown (red). Neuron parameters: a ∼ N (0.02,0.001), b ∼ N (0.2,0.01), c ∼ N (− 55,1), d ∼ N (2,0.1), κ = 0.8, f_min = 0.2, and I_ext ∼ N (6,0.2). Network parameters: α = 10^{− 1}, $\bar{g}=0.002$, E_syn = 0, p_con = 0.03, and w₀ = 1.

This is explained as follows: as all connections in the network are excitatory, the V-nullcline is effectively lifted during the network bursts due to the external input from synapses. When the network burst terminates, the V-nullcline drops to its initial position. This creates the larger range of u which must be traversed for both the spiking and quiescent phase of the bursting, resulting in more spikes and longer interburst intervals respectively.

These network bursts are, apart from a small difference in interburst-frequency, accurately reproduced by the network reduction (11). Although the interburst-frequency is off, the last network burst (starting at t=900 ms) reveals that the amplitude, shape, and duration are each mimicked in great detail by the neural mass model.

Network with rebound bursts

One of the key features of the reduction is that the subthreshold dynamics of the Izhikevich neuron are left intact, and, as a consequence, it is able to generate post-inhibitory rebound potentials, c.f. Appendix A: Single cell comparison with post-inhibitory rebound potentials. To illustrate this distinct feature of the model, simulations are performed with a network whose rhythms relies heavily on rebound activity.

A network consisting of three populations of neurons, one inhibitory (I) of 100 neurons and two excitatory (E₁ and E₂) of 300 neurons each, is set up with connections as in Figure 5A. All parameters of the network are given in Appendix B: Details of three population network, but here we give a global outline of its properties. Due to inhomogeneities, neurons in E₁ fire either sporadically or are quiescent otherwise and, furthermore, these neurons are capable of firing rebound bursts in response to an inhibitory barrage by neurons in population I. Population E₂ contains intrinsically bursting neurons, which, similar to the example in Section ‘One excitatory population’, produce regular network bursts. The slow parameter a, however, is decreaed such that the burst-frequency is lower. These bursts will activate the inhibitory neurons in population I.

Figure 5

Simulations for a network with post-inhibitory rebound potentials. (A) depicts the network structure of the three populations: two excitatory populations, E₁ and E₂, each consisting of 300 cells and one inhibitory population I of 100 cells. (B) shows raster plots for 20 neurons within the populations I, E₁ and E₂ respectively. Note that some cells in E₁ are quiescent. The activation of synapses S for each population is shown in (C), where (blue) corresponds to the spiking neuron network and (red) represents the reduced model. See text for a description of network activity. E_inh = − 80 and other parameters in Appendix B: Details of three population network.

A simulation of this particular network is shown in Figure 5B by means of raster plots. The pattern arises because a population bursts in E₂ gives rise to activation of I, which, in turn, causes a synchronized post-inhibitory burst in E₁. This burst in E₁ is large enough to initiate the next population burst in E₂ before intrinsic mechanisms would. This completes the cycle. Figure 5C shows the summed synaptic variables of each population in blue.

The population reduction from this model follows readily from the single population discussed before; for clarity though, the model and its parameters are stated in Appendix B: Details of three population network. The time evolution of the resulting reduced model is shown in red in Figure 5C. Apart from a change in amplitude, the frequency of the reductio60n matches closely with the original spiking neuron network. The fact that the period of this solution matches better than in the single population example, see previous section, could be explained by the following: In the single population model, c.f. Figure 4, the network bursts are initiated when the slow variable u reaches a certain threshold value and are, therefore, dependent on the value of the slow variable at burst termination. Since the network of Izhikevich neurons increases the slow variable u with discrete steps d, a discrepancy is to be expected with respect to the reduced model, which increases u in a time-continuous manner. Although the slow variable plays a role in the rebound bursts in the network at hand, the timing of the (rebound-)bursts depends more so on the synaptic time constants — which are identical in both models.

Next, we note that, apart from a difference in amplitudes, the shapes of the bursts in both populations I and E₂ are accurately mimicked. The discrepancy is largest in population E₁, with respect to both amplitude and burst duration. It appears that it is this gain in E₁ in amplitude which gives rise to an increased excitatory input onto E₂, resulting in a burst with increased magnitude. Subsequently, the stronger burst in E₂ explains the amplitude gain of the inhibitory population I.

The natural question to ask now is why the post-inhibitory rebound bursts in E₁ are poorly reproduced by the reduced model. The fact that some neurons are quiescent during the bursts, as follows from the raster plot in Figure 5, suggests that inhomogeneities play a role in this; either due to variations in single cell parameters, or to differences in number of network connections. Indeed, further investigation has revealed that, at the level of a single cell, the generation of rebound bursts is particularly sensitive to variations in parameter b. Simulations (not shown) of the spiking neuron network in which for all cells b is set to a constant — rather than drawn from a normal distribution, c.f. Appendix B: Details of three population network — show that the bursts of neurons in E₁ are more alike than in Figure 5. As a consequence, the result lies closer to the response of the reduced model (in which all inhomogeneities of the original network are averaged); not only for population E₁, also for I and E₂.

Finally, we illustrate the models capability to distinguish between different types of inhibition. Figure 6 shows simulation results for an identical network structure as before, with the exception that the reversal potential of the inhibitory synapses is set to E_inh=−65. For this value is close to the resting membrane potentials of the neurons in E₁, the inhibition becomes effectively of the shunting type. Because shunting inhibition does not alter the membrane potential of a cell, it is unable to activate the slow currents responsible for the post-inhibitory rebound potentials and, hence, the neuron is unable to fire such action potentials. This matches with the results in Figure 6B: the raster plot of E₁ shows not only less activity, but it appears more irregular too. This is better seen by the synaptic variable in Figure 6C, which is 3–4 times smaller than in 5C. Here it is also visible that this background activity of E₁ is absent during inhibitory activity of population I, hence the shunting inhibition successfully suppresses this activity.

Figure 6

Similar to Figure 5, except that the post-inhibitory rebound potentials are taken away by raising the reversal potential of the inhibitory synapses to E _inh =− 65 . (A) depicts the network structure, with the reduced inhibition denoted by the dashed line. Raster plots of the three populations are in (B) and the net synaptic activation of the populations in (C). Here, the original network is plotted in (blue) and the lumped model in (red).

Without the rebound potentials of E₁, neurons in E₂ receive less excitatory input and the frequency of network bursts is merely determined by the intrinsic dynamics of population E₂. that the neurons in E₂ are not subject to bursts, it remains unclear why the difference in amplitude between the network and corresponding reduction is still in place. A partial explanation is that the excitatory background activity in E₁, which is missed by the population reduction, lowers the threshold of neurons, such that network bursts are generated faster than without input. With a lower inter-burst time, the slow recovery variable u has not decayed as much as it would have if the bursts were further apart. As a result of this, the number of spikes within the burst is expected to be lower. Although this effect plays a role, it is expected that other, yet unidentified mechanisms, are also involved.

Discussion

The work presented here shows how the global activity of networks of intrinsically bursting Izhikevich neurons can be described with a lumped model. This lumped model relies, as the majority of reduced models, on the firing rate of individual neurons. Assisted by phase plane analysis of a single neuron an additional term is proposed, which abolishes the need for discontinuous resets. Effectively, the new term divides the phase space in two regions, corresponding with bursting and non-bursting. Within the bursting regime, the neuron’s firing can, once a separation of timescales is assumed, be determined analytically, while in the non-bursting regime, the neuron’s original subthreshold behavior is kept intact. We note that the resulting model is Lipschitz smooth and its phase portrait shares key features with the Fitzhugh-Nagumo [18, 19] and Morris-Lecar [20] reductions. In our case, however, ‘active periods’ correspond to an entire burst of action potentials rather than a single pulse. Using simulations, the proposed single cell reduction is shown to match well with the original Izhikevich model. Although a small mismatch is present between the models’ bursting frequency, both duration and number of action potentials within the burst are accurately represented.

With the single cell model into place, the scene is set for lumping networks of bursting neurons. The cells’ planar dynamics are augmented with another state variable, which corresponds to the activation of synapses within the population. By employment of assumptions common in neural mass modeling, a reduced network model is formulated, which consists of three state variables representing the population averages of membrane potential, slow current, and synaptic activation. Similar as for the single cell model, simulations are performed to assess the quality of the newly proposed model. For one population of excitatory intrinsically bursting neurons, the reduced model matches closely with the network dynamics. Indeed, the network bursts, which last longer than the bursts of isolated cells, are mimicked by the lumped model in terms of amplitude, duration, and shape. While the model’s burst frequency is lower than the corresponding spiking neuron network, both models’ bursting frequencies in the network are notably lower than the isolated cells.

Here it is important to point out that a very similar result for a single population has been obtained in [15]. The fact that their formulation is expressed in only two state variables, suggests that the third state variable in our formulation is superfluous. We believe, however, that this third state variable, which corresponds to the mean membrane potential, is essential for keeping the neurons’ subthreshold behavior intact. To demonstrate this particular trait of our reduction, a network is simulated whose rhythm depends critically on the subthreshold behavior (c.f. Figure 5). More precisely, neurons belonging to one of the three populations are capable of generating post-inhibitory rebound bursts, which affect the network’s behavior at a global scope. Apart from a discrepancy in amplitude, simulations show that the reduced model is able to accurately describe the network activity. By raising the inhibitory reversal potential, such that the inhibition becomes shunting, we show that this particular rhythm is, in both the original spiking neuron network and the reduced model, indeed reliant on the presence of post-inhibitory rebound bursts (c.f. Figure 6).

The single cell reduction proposed in this article considers the bursting Izhikevich neuron merely as an example. We indeed assumed the reset value c to be higher than the local minimum of the V nullcline (i.e. c > − 62.5) and d>0. Although this particular formulation might not be applicable when the parameters at hand do not correspond with the bursting type, we believe that the procedure here can be generalized for all relevant choices of parameters. Going even further, it appears to us that, besides the Izhikevich and other integrate-and-fire models, also higher dimensional models can be reduced in a similar manner; that is, leaving the subthreshold dynamics intact by only adjusting the dynamics in the suprathreshold regime. In this regime, one should seek to eliminate the (very) fast dynamics and determine their contribution on the slow variables. Yet, these fast variables can play a prominent role in the subthreshold dynamics (e.g. to generate post-inhibitory rebound potentials). Elimination of all fast variables could therefore have an unanticipated (and hence undesirable) effect on the subthreshold behavior. So caution is required when constructing firing rate reductions.

Endnote

^a Software freely available via ModelDB entry 154739: http://senselab.med.yale.edu/modeldb/showmodel.asp?model=154739.

Appendix A: Single cell comparison with post-inhibitory rebound potentials

Figure 7 shows a concrete example of post-inhibitory rebound potentials in both the spiking model (black) and rate model (gray). An negative current of −10 is injected between 200–500 ms, which lowers the membrane potential. After the inhibiting current is taken away, the spiking neuron fires three action potentials. Not only the count of action potentials is accurately reproduced by the rate model, also the time course of all other state variables is mimicked precisely.

Figure 7

Negative current induces post-inhibitory rebound potentials in both the spiking neuron and the reduction. Top, middle, and bottom diagrams correspond with V, u and N (c.f. (7)) respectively. Parameters as follows: a=0.03, b = 0.25, c = − 52, d = 2, κ = 0.8, and f_min = 0.1. Furthermore, I_ext = − 10 between 200 – 500 ms, after which both models fire (the equivalent of) three action potentials.

Appendix B: Details of three population network

This appendix provides the details of the three population network described in Section ‘Network with rebound bursts’. The spiking neuron is characterized first, followed by the reduced model.

B.1 Spiking neuron network

Three population of neurons are generated, one inhibitory I and two excitatory E₁ and E₂. Parameters for neurons within each population are i.i.d. according to the normal distributions given in Table 1.

Table 1 Population parameters

Connections between neurons are determined randomly, depending on which population they belong to. The population-specific probability p of forming a connection from a source neuron to a destination neuron is denoted in Table 2; as is the related maximal conductance $\bar{g}$. All connections formed have synaptic weight w₀=1, hence weight distribution W is given by

$$W=\left\{\begin{array}{l}1\text{w.p.}p\hfill \\ 0\text{w.p.}1-p\end{array}\right.$$

Table 2 Connectivity parameters

Synaptic dynamics of all cells, both excitatory and inhibitory, are described by the first order process

$$\frac{\mathit{\text{ds}}}{\mathit{\text{dt}}}+\mathrm{\alpha s}=\delta (t-T)$$

with T the time of a presynaptic spike. Values of α, as well as the corresponding reversal potential of the resulting postsynaptic current, are in Table 3. The total synaptic current arriving at a particular cell is now determined as in (8).

Table 3 Synaptic parameters

B.2 Reduction

Following the reduction procedure outlined in Section ‘A population of Izhikevich neurons’, each population of neurons is reduced to a system of three differential equations. For k ∈ {i,1,2} representing the populations I, E₁, E₂ respectively, the state variables are denoted by V _k , u _k , and S _k . Parameters for the right hand side functions f _k , g _k , and h _k , c.f. (11), are for each population set to their expected values given in Table 1. Interactions between populations are exclusively via the synaptic current I_syn, which we state in full detail for clarity. Let I_syn,k represent the postsynaptic current arriving at cells of population k = {i,1,2}, then one finds:

$$\begin{array}{ll}{I}_{\text{syn},i}& ={\bar{g}}_{i2}{p}_{i2}{S}_2(V_i-E_{\text{exc}})\\ I_{\text{syn},1}& ={\bar{g}}_{1i}{p}_{1i}{S}_i(V_1-E_{\text{inh}})\\ I_{\text{syn},2}& =({\bar{g}}_{21}{p}_{21}{S}_1+{\bar{g}}_{22}{p}_{22}{S}_2)(V_2-E_{\text{exc}})\end{array}$$

where p _kl and g _kl for k,l ∈ {i,1,2} are as in Table 2 and we used (12).