Abstract

The intrinsic electrophysiological properties of single neurons can be described by a broad spectrum of models, from realistic Hodgkin-Huxley-type models with numerous detailed mechanisms to the phenomenological models. The adaptive exponential integrate-and-fire (AdEx) model has emerged as a convenient middle-ground model. With a low computational cost but keeping biophysical interpretation of the parameters, it has been extensively used for simulations of large neural networks. However, because of its current-based adaptation, it can generate unrealistic behaviors. We show the limitations of the AdEx model, and to avoid them, we introduce the conductance-based adaptive exponential integrate-and-fire model (CAdEx). We give an analysis of the dynamics of the CAdEx model and show the variety of firing patterns it can produce. We propose the CAdEx model as a richer alternative to perform network simulations with simplified models reproducing neuronal intrinsic properties.

1  Introduction

Computational modeling of large-scale networks requires compromising between biophysical realism and computational cost. This requirement might be satisfied by many two-variable models (Morris & Lecar, 1981; Krinskii & Kokoz, 1973; Fitzhugh, 1961; Izhikevich, 2003; Brette & Gerstner, 2005). In particular, two-variable models that are largely used and studied are the Izhikevich model (Izhikevich, 2003) and the adaptive exponential integrate-and-fire (AdEx) model (Brette & Gerstner, 2005; Naud, Marcille, Clopath, & Gerstner, 2008; Touboul & Brette, 2008). The first variable of these models corresponds to the membrane voltage; the second, usually with slower kinetics, corresponds to neural adaptation and allows achieving more complex dynamics and firing patterns, which can be observed in neural recordings. In both models, the second variable has a form of an additional current flowing into the neuron, which may lead to unrealistic changes of membrane voltage, especially in the case of long and intense neuronal firings such as that observed during seizures (McCormick & Contreras, 2001). Here, we describe this limitation and propose a modification of the second variable dynamics to oppose it. The adaptation variable in our model has the form of a conductance, introducing the conductance-based adaptive experimental integrate-and-fire (CAdEx) model. Previous work has used conductance-based adaptation, but these models used were either more simplified (Treves, 1993) or much more detailed Hodgkin-Huxley-type models with complex channel gating dynamics (Connor & Stevens, 1971). A universal phenomenological model linking f-I curves with adaptation dynamics (Benda & Herz, 2003) has a useful general form, but it did not specify voltage dynamics and post-spike adaptation. The CAdEx model keeps the simplicity of the two-variables models, while extending the repertoire of possible subthreshold dynamics.

2  Conductance-Based Adaptive Exponential Model

The adaptive exponential integrate-and-fire model can reproduce many types of firing patterns: tonic spiking cells, regular spiking adaptive cells, bursting cells, and others (Naud, Marcille, Clopath, & Gerstner, 2008; Clopath, Jolivet, Rauch, Lüscher, & Gerstner, 2007; Touboul & Brette, 2008). However, because the adaptation in the AdEx model is a current, a strong and unrealistic hyperpolarization of the cell can appear after a period of prolonged firing (see Figure 1). Moreover, subthreshold adaptation in this model is linear, which means that in the periods of strong hyperpolarization, the adaptation current does not deactivate and can remain unrealistically strong. For instance, modeling M-current mediated adaptation in this model poses a problem because M-channels are mostly closed when the membrane voltage remains below -60 mV (Brown & Adams, 1980). A similar problem happens for strong depolarization due to lack of saturation of the adaptation current. A corresponding problem exists for inward depolarizing currents such as the Ih current: hyperpolarization-activated cyclic nucleotide-gated (HCN) channels close above -70 mV and saturate below -80 mV (Banks, Pearce, & Smith, 1993).
Figure 1:

Comparison between AdEx and CAdEx models. (a) Both cells were injected with the same current pulse. Their firing rates f and adaptation parameters A were the same during pulse, f=30Hz and A=0.03. However, after pulse, the AdEx cell hyperpolarizes below -150mV while the CAdEx cell hyperpolarizes to proximity of the adaptation reversal potential, EA=-70mV. The mathematical definition of the adaptation parameter, A, is given in section 7. (b) Network simulations. The network comprises 1000 cells with 80% excitatory neurons with the same parameter values as for the corresponding single-cell simulations. More details about these networks are given in appendix A5. Upper plots: Mean membrane voltage of excitatory neurons. Below: Corresponding raster plots for 400 excitatory neurons (below blue line) and 100 inhibitory neurons (above blue line).

Figure 1:

Comparison between AdEx and CAdEx models. (a) Both cells were injected with the same current pulse. Their firing rates f and adaptation parameters A were the same during pulse, f=30Hz and A=0.03. However, after pulse, the AdEx cell hyperpolarizes below -150mV while the CAdEx cell hyperpolarizes to proximity of the adaptation reversal potential, EA=-70mV. The mathematical definition of the adaptation parameter, A, is given in section 7. (b) Network simulations. The network comprises 1000 cells with 80% excitatory neurons with the same parameter values as for the corresponding single-cell simulations. More details about these networks are given in appendix A5. Upper plots: Mean membrane voltage of excitatory neurons. Below: Corresponding raster plots for 400 excitatory neurons (below blue line) and 100 inhibitory neurons (above blue line).

To overcome these problems, we propose a model with a conductance-based adaptation gA and with a sigmoid dependence of subthreshold adaptation on voltage. We named it the conductance-based adaptive exponential model (CAdEx). The equations of the model are as follows:
CdVdt=gL(EL-V)+gLΔTexpV-VTΔT+gA(EA-V)+IsτAdgAdt=g¯A1+expVA-VΔA-gA,
(2.1)
with an after-spike reset mechanism,
ifVVDthenVVRgAgA+δgA,
(2.2)
where C is the membrane capacitance and gL and EL are, respectively, the leak conductance and reversal potential. VT is the spike threshold and ΔT is the slope of the spike initiation, EA is the reversal potential of the adaptation conductance, and Is is an input current. In equation 2.2, τA is the time constant of adaptation, g¯A is the maximal subthreshold adaptation conductance, VA is the subthreshold adaptation activation voltage, and ΔA is the slope of subthreshold adaptation.

A spike is initiated when V approaches VT and the exponential term escalates. After reaching a detection limit VD, the voltage is reset to the reset potential VR, and it remains at this value during a refractory period Δtref. After each spike, gA is incremented by a quantal conductance δgA.

The sigmoidal subthreshold adaptation function is always positive, g¯A0, and it can be a monotonically increasing, ΔA>0, or decreasing, ΔA<0, function of membrane voltage. In the simulations, the value of the detection limit has been set to VD=-40mV and the value of the refractory period to Δtref=5ms.

In appendix A4, we exemplify the parameterizations for IM and Ih currents based on biophysical data (Brown & Adams, 1980; Banks et al., 1993).

The problem of unrealistic hyperpolarization in the AdEx model does not appear in the CAdEx model (see Figure 1). The corresponding problem of unrealistic hyperpolarization can appear also at the network level for AdEx neurons (see Figure 1b). Although both models can be tuned to exhibit very similar raster plots at the network level, the hyperpolarization problem disappears in networks of CAdEx neurons (see Figure 1). For a negative adaptation parameter a, corresponding to the accelerating firing pattern (Brette & Gerstner, 2005), in the AdEx model, an additional problem of infinite hyperpolarization can appear, as shown in appendix A6. Thus, the CAdEx model may be advantageous when the phenomena considered depend on the membrane potential.

All simulations of the models were done in Brian2 neural simulator and Python 3 programming language. The code, which allows running the CAdEx model simulations can be downloaded from https://github.com/neural-decoder/cadex.

Figure 2:

(a) Nullclines of the CAdEx system for three parameterizations (left) for ΔA>0 and at most two fixed points, (center) for ΔA<0 and at most four fixed points; and (right) for ΔA>0 and at most four fixed points. The blue line is a V-nullcline and red line is gA-nullcline. Red filled circles indicate stable fixed points, white filled circles indicate unstable fixed points, and half-filled circles indicate neutral points. Arrows indicate the direction of the vector field. (b) Corresponding S(V) functions. Intersections with lines of constant input current Is indicate the positions of the fixed points.

Figure 2:

(a) Nullclines of the CAdEx system for three parameterizations (left) for ΔA>0 and at most two fixed points, (center) for ΔA<0 and at most four fixed points; and (right) for ΔA>0 and at most four fixed points. The blue line is a V-nullcline and red line is gA-nullcline. Red filled circles indicate stable fixed points, white filled circles indicate unstable fixed points, and half-filled circles indicate neutral points. Arrows indicate the direction of the vector field. (b) Corresponding S(V) functions. Intersections with lines of constant input current Is indicate the positions of the fixed points.

3  Dynamical Analysis of the Model

3.1  Fixed Points and Bifurcations

The CAdEx is a nonlinear dynamical system with a discontinuous post-spike reset mechanism. The V and gA nullclines are as follows (see Figure 2):
gA=gLEL-VV-EA+gLΔTV-EAexpV-VTΔT+IsV-EA,gA=g¯A1+expVA-VΔA.
(3.1)

The shape and the location of the V-nullcline depend on the input current Is (see Figure 11 in appendix A1). The V-nullcline is divided by the vertical asymptote V=EA into the left and right branches (see Figure 2a.)

The intersections of nullclines give the fixed points of the system. From equation 3.1, the solutions of the following equation,
Is=-gL(EL-V)-gLΔTexpV-VTΔT-g¯A1+expVA-VΔA(EA-V),
(3.2)
give the membrane voltages of the fixed points.

The analysis of equation 3.2, which we can write in a simpler form defining a new function S(V) as Is=S(V), gives us the number of fixed points for a given input current (see Figure 2). The function S(V) tends to - for V± and it can have:

  1. A single maximum Vmax, which corresponds to the configuration with at most two possible intersections between nullclines. Let ISN=S(Vmax). For a current input Is<ISN, there are therefore two fixed points V-(Is) and V+(Is); for Is=ISN, there is one fixed point; and for Is>ISN, there is no fixed point. ISN is therefore the current above which the system always spikes spontaneously. The system can start to spike spontaneously for Is<ISN if all fixed points lose stability before reaching ISN in an Andronov-Hopf bifurcation (see below).

  2. Two local maxima Vmax1,2, and one local minimum Vmin, which corresponds to maximally four possible fixed points. For input current Is<S(Vmin), there are two fixed points, one stable and one unstable. For S(Vmin)<Is<minSVmax1,2, there are four fixed points, from which two can be stable. For minSVmax1,2<Is<maxSVmax1,2, there are two fixed points, and for Is>maxSVmax1,2, there is no fixed point and the neuron fires spontaneously.

To analyze the bifurcations of the system, we need to consider the behavior of the Jacobian matrix at the fixed points i as a function of the input current, Ji(Is) (see appendix A2).

The trajectories in the phase space of trace trJi(Is)=tr(Ji(Is)) and determinant detJi(Is)=det(Li(Is)) of the Jacobian give us the type of bifurcation (see Figures 3 and 4).

In case 1, when there are maximally two possible fixed points, the value of trJ at their merging point,
trJ(ISN)=-1τA-g¯A(EA-V±(ISN))CΔAexpVA-V±(ISN)ΔA1+expVA-V±(ISN)ΔA2,
(3.3)
determines the type of bifurcation (see appendix A2). If trJ(ISN)<0, then the bifurcation is of a saddle-node type. If trJ(ISN)=0, the bifurcation is of a Bogdanov-Takens type. If trJ(ISN)>0, the bifurcation is of Andronov-Hopf type. Since the merging point can be located only on the right branch of the V-nullcline, EA<V(ISN); then if ΔA<0, the bifurcation is always of a saddle-node type. If ΔA>0, all types are possible depending on the sign of trJ(ISN). If g¯A=0, the bifurcation is also always of a saddle-node type.
Figure 3:

Bifurcations of the CAdEx system with at most two fixed points. (a) Trajectories of the fixed points on the Poincaré diagram for increasing input current Is. Solid lines denote trajectories of left, initially stable, fixed point with lower membrane voltage; dashed lines show the trajectories of right, unstable fixed point with higher membrane voltage (shown partially). For positive ΔA, blue lines, the stable fixed point loses stability with an Andronov-Hopf bifurcation, g¯A=10and20nS. For negative ΔA, red lines, the stable fixed point merges with the saddle point in a saddle-node bifuraction, g¯A=5and10nS. Without subthreshold adaptation, g¯A=0nS, the black line, the system undergoes a saddle-node bifurcation. (b, c) Response of the cell to a brief current pulse: panel b near the Andronov-Hopf bifurcation and panel c near the saddle-node bifurcation. The parabola, the dashed black line, is a function 14(trJ)2; inside the area limited by the parabola, the system can oscillate around equilibria.

Figure 3:

Bifurcations of the CAdEx system with at most two fixed points. (a) Trajectories of the fixed points on the Poincaré diagram for increasing input current Is. Solid lines denote trajectories of left, initially stable, fixed point with lower membrane voltage; dashed lines show the trajectories of right, unstable fixed point with higher membrane voltage (shown partially). For positive ΔA, blue lines, the stable fixed point loses stability with an Andronov-Hopf bifurcation, g¯A=10and20nS. For negative ΔA, red lines, the stable fixed point merges with the saddle point in a saddle-node bifuraction, g¯A=5and10nS. Without subthreshold adaptation, g¯A=0nS, the black line, the system undergoes a saddle-node bifurcation. (b, c) Response of the cell to a brief current pulse: panel b near the Andronov-Hopf bifurcation and panel c near the saddle-node bifurcation. The parabola, the dashed black line, is a function 14(trJ)2; inside the area limited by the parabola, the system can oscillate around equilibria.

In case 2, when there are maximally four possible fixed points, two fixed points are created in a blue sky bifurcation when the input current reaches a minimum of S(V) function, IBS=S(Vmin) (see Figure 4). A blue sky bifurcation is a saddle-node bifurcation for an increasing input current in which two fixed points appear. If at the creation point trJ(IBS)<0, a pair of stable and unstable fixed points is created. If trJ(IBS)0, a pair of unstable fixed points is created. trJ(IBS) has the form
trJ(IBS)=-1τA-g¯A(EA-Vmin)CΔAexpVA-VminΔA1+expVA-VminΔA2.
(3.4)
If ΔA<0, the blue sky bifurcation point can be located only on the right branch of the V-nullcline and consequently EA<V(IBS) (see the center plot in Figure 2). From equation 3.4, trJ(IBS)<0; hence, a stable fixed point and an unstable fixed point are created. If ΔA>0, the blue sky bifurcation point can be located on case 1, the left branch of the V-nullcline for a strongly hyperpolarized neuron, V(IBS)<EA (see Figure 11 in appendix A1). In this case, trJ(IBS)<0 and a stable fixed point can be created. Case 2 is the right branch of the V-nullcline, V(IBS)>EA (right plot in Figure 4). In this case, both a pair of unstable fixed points or a pair of a stable and an unstable fixed point can be created depending on the sign of trJ(IBS). If a stable fixed point is created, it can subsequently lose its stability in an Andronov-Hopf bifurcation (see the upper right plot in Figure 4).
Figure 4:

Bifurcations of the CAdEx system with at most four fixed points. Trajectories of the fixed points on the Poincaré diagram are shown for increasing input current Is. Trajectories are shown for positive ΔA, blue lines, and for negative ΔA, red lines. Two fixed points are created in a blue sky bifurcation that subsequently merge with a stable and an unstable fixed point in a saddle-node bifurcation. (Upper right) The stable fixed point created in a blue sky bifurcation can lose stability in an Andronov-Hopf bifurcation before merging with the unstable fixed point in a saddle-node bifurcation. The parabola, the dashed black line, is a function 14(trJ)2; inside the area limited by the parabola, the system can oscillate around equilibria.

Figure 4:

Bifurcations of the CAdEx system with at most four fixed points. Trajectories of the fixed points on the Poincaré diagram are shown for increasing input current Is. Trajectories are shown for positive ΔA, blue lines, and for negative ΔA, red lines. Two fixed points are created in a blue sky bifurcation that subsequently merge with a stable and an unstable fixed point in a saddle-node bifurcation. (Upper right) The stable fixed point created in a blue sky bifurcation can lose stability in an Andronov-Hopf bifurcation before merging with the unstable fixed point in a saddle-node bifurcation. The parabola, the dashed black line, is a function 14(trJ)2; inside the area limited by the parabola, the system can oscillate around equilibria.

The trajectories of the equilibria and the corresponding dynamics in the CAdEx model are much different than they are in the AdEx model. In the AdEx model, at most only two fixed points are possible, and the trajectories of equilibria on the Poincaré diagram are linear (see appendix A.7).

3.2  Subthreshold Oscillations

The system can oscillate around equilibrium if the equilibrium is stable, that is, trJ(Is)<0 and detJ(Is)>0, and the eigenvalues of the Jacobian at the fixed point have an imaginary part, trJ2-4detJ<0. In that case, the frequency of oscillations is given by ν=14π4detJ-trJ2 (see Figure 5 and appendix A2).

The sustained subthreshold oscillations due to the emergence of a limit cycle can also be observed. Some examples in the case of multistability are given in the next section.

4  Multistability

In this section, we present a list of examples of multistability. Such behaviors can be observed due to two properties of the model:

  • Subthreshold multistability. Due to the nonlinear form of the subthreshold adaptation, the system may exhibit two or four fixed points. With four fixed points, various configurations are possible: with a positive or negative slope (ΔA) of the adaptation nullcline. The stability of these fixed points, depending on the input current, is shown in Figure 4.

  • Superthreshold stability. The model has a discontinuity that occurs after a spike and corresponds to the reset of the membrane voltage and the increment of adaptation conductance. It can lead to steady-state firing behavior, which can be treated as an attractor of the model. The reset discontinuity has, in many aspects, the same effect as a third variable. It allows chaotic behaviors (as described in the next section) normally impossible in a continuous two dimensional system.

Here, we give three examples of multistability: (1) with four fixed points and a negative slope of adaptation, (2) with four fixed points and a positive slope of adaptation, and (3) with two fixed points and a positive slope of adaptation:

  1. With a negative and small enough value of ΔA, the nullclines can have four crossings. In this case, the system has two stable fixed points with different values of the membrane potential. In this situation, three stable steady states are observed (see Figure 6a): two possible resting states of membrane potential and a self-sustained spiking state. Small perturbations permit a switch between these steady states, as shown in Figure 6b. By applying synaptic noise through conductance-based inhibitory and excitatory synapses (see appendix A3), transitions between these stable states can be observed (see Figure 6c).

  2. With a positive and small enough value of ΔA, the system can also have four fixed points. In this case, the emergence of a stable limit cycle is observed. It leads to three possible stable behaviors: resting state, subthreshold oscillations, and regular spiking, as shown in Figure 7a. As in the previous condition, transitions between these states can be observed under synaptic bombardment (see Figure 7b).

  3. In a situation with two fixed points, the emergence of an unstable limit cycle is observed, leading to bistability near the threshold between spiking and nonspiking damped oscillations, as shown in Figure 8.

Figure 5:

(a) Subthreshold oscillations of the membrane voltage around the stable equilibrium for different constant input currents. (b) Frequency of subthreshold membrane voltage oscillations in the CAdEx model as a function of a constant input current. The crosses signify numerical results and the dashed line a theoretical prediction, ν=14π4detJ-trJ2. Here, a pulse of current of amplitude 10 pA and duration of 10 ms was applied when the membrane potential was at a stable equilibrium.

Figure 5:

(a) Subthreshold oscillations of the membrane voltage around the stable equilibrium for different constant input currents. (b) Frequency of subthreshold membrane voltage oscillations in the CAdEx model as a function of a constant input current. The crosses signify numerical results and the dashed line a theoretical prediction, ν=14π4detJ-trJ2. Here, a pulse of current of amplitude 10 pA and duration of 10 ms was applied when the membrane potential was at a stable equilibrium.

Figure 6:

Multistability observed with a negative slope of gA. (a) The dynamics for a single neuron with two stable fixed points and a stable limit cycle due to the reset. (b) Transitions between possible steady states due to excitatory or inhibitory short-time perturbations. (c) Transition between the up and down state of a single neuron receiving input through conductance-based synapses.

Figure 6:

Multistability observed with a negative slope of gA. (a) The dynamics for a single neuron with two stable fixed points and a stable limit cycle due to the reset. (b) Transitions between possible steady states due to excitatory or inhibitory short-time perturbations. (c) Transition between the up and down state of a single neuron receiving input through conductance-based synapses.

Figure 7:

Multistability observed with a positive slope of gA. (a) The dynamics for a single neuron with a stable fixed point: a stable limit cycle due to the subthreshold dynamics and a limit cycle due to the reset. (b) Transitions between the up and down states of single neurons receiving input through conductance-based synapses.

Figure 7:

Multistability observed with a positive slope of gA. (a) The dynamics for a single neuron with a stable fixed point: a stable limit cycle due to the subthreshold dynamics and a limit cycle due to the reset. (b) Transitions between the up and down states of single neurons receiving input through conductance-based synapses.

Figure 8:

Multistability observed with two fixed points. (a) Damped subthreshold oscillations are observed near the spiking behavior. (b) Enlarged concerned region (blue square). Two simulations initiated: inside (red) and outside (dark gray) the unstable limit cycle.

Figure 8:

Multistability observed with two fixed points. (a) Damped subthreshold oscillations are observed near the spiking behavior. (b) Enlarged concerned region (blue square). Two simulations initiated: inside (red) and outside (dark gray) the unstable limit cycle.

5  Firing Patterns

The CAdEx model is able to reproduce a large repertoire of intrinsic electrophysiological firing patterns. In this section, we detail the dynamics of firing patterns in CAdEx, reproducing those found in AdEx. In order to facilitate the use of this model to reproduce specific behaviors, we describe which parameters influence specific firing patterns. The initial values of the two variables in the simulations were V(0)=-60mV for the membrane potential and gA(0)=0nS for the adaptation, except when ΔA<0 when gA is initialized at gA(-60mV). For the sets of parameters used in this section, the system does not have fixed points. Therefore, independent of the starting point location, the system converges to the same steady-state firing. However behaviors in the transient regime may change, which is crucial for spike frequency adaptation. The parameters for each pattern of Figure 9 are given in Table 1.

Table 1:
Parameters of the CAdEx Model for the Firing Patterns Shown in Figure 9.
Cm [pF]EA [mV]EL [mV]Is [pA]VA [mV]ΔA [mV]
Adaptive spiking 200 −70 −60 200 −50 
Tonic spiking 200 −70 −70 192 −45 
Bursting 200 −60 −58 150 −45 
Delayed bursting 200 −70 −60 100 −45 
Accelerated spiking 200 −70 −60 130 −60 −5 
Chaotic-like spiking1111 200 −70 −58 90 −40 
 VR [mV]1111111 VT [mV] δgA [nS] g¯A [nS] gL [nS] τA [ms] 
Adaptive spiking1111 −55 −50 10 10 200 
Tonic spiking −56 −50 10 40 
Bursting −46 −50 10 10 200 
Delayed bursting −46 −50 12 100 
Accelerated spiking −58 −48 10 300 
Chaotic-like spiking −47 −50 10 10 25 
Cm [pF]EA [mV]EL [mV]Is [pA]VA [mV]ΔA [mV]
Adaptive spiking 200 −70 −60 200 −50 
Tonic spiking 200 −70 −70 192 −45 
Bursting 200 −60 −58 150 −45 
Delayed bursting 200 −70 −60 100 −45 
Accelerated spiking 200 −70 −60 130 −60 −5 
Chaotic-like spiking1111 200 −70 −58 90 −40 
 VR [mV]1111111 VT [mV] δgA [nS] g¯A [nS] gL [nS] τA [ms] 
Adaptive spiking1111 −55 −50 10 10 200 
Tonic spiking −56 −50 10 40 
Bursting −46 −50 10 10 200 
Delayed bursting −46 −50 12 100 
Accelerated spiking −58 −48 10 300 
Chaotic-like spiking −47 −50 10 10 25 

Adaptive spiking. The adaptation can cause a decrease of spiking frequency (see Figure 9a). Correspondingly, the value of gA increases until it reaches a steady state. The spike frequency adaptation is mainly affected by the spike triggering parameters (VT and ΔT); by the membrane properties, C and gL; as well as by the adaptation parameters, g¯A and strongly by δgA.

Tonic spiking. In the absence of adaptation, that is, g¯A=0, the model exhibits tonic spiking. However, tonic spiking can also be observed with nonzero adaptation, as in Figure 9b, if the steady-state firing can be reached after the first spike. The steady-state spiking frequency is influenced by the membrane capacitance, C, the leak conductance, gL, and the slope of spike initiation ΔT, and also by the adaptation time constant τA. The after-hyperpolarization following spiking depends on both: reset value VR and adaptation time constant τA.

Bursting. To obtain bursting behaviors (see Figure 9c), the reset value, VR, has to be higher than the voltage of the minimum of the V-nullcline. The system spikes until it crosses the V-nullcline, ending the burst. Then, above the V-nullcline, the system is driven to lower values of adaptation and membrane potential. After a second crossing of the V-nullcline, the system spikes again, starting a new bursting cycle. Bursts can be characterized by intra- and interburst activities. The interburst activity is affected by adaptation; g¯A and τA strongly determine the interburst time interval and after-burst hyperpolarization.

Delayed spiking. All firing patterns can occur with a time delay. Figure 9d gives an example for bursting and Figure 9b an example for low-frequency spiking. The distance between V and gA nullclines is determined by Is. By changing this distance and the time constant τA, it is possible to obtain a region of slow flow (small dVdt) and consequently to increase the delay and the interevent interval. The reset VR and δgA also affects the time required for the system to go around the V nullcline and then to the next event.

Accelerated spiking. A negative slope of gA (ΔA<0) associated with δgA=0, leads to firing rate acceleration, as shown in Figure 9e.

Chaotic-like spiking. Because the CAdEx model is not continuous due to the after-spike reset of voltage and incrementation of adaptation, chaotic-like spiking behavior can be observed, as shown in Figure 9f. These phenomena occur when the reset is close to the right branch of the V nullcline. To verify that this irregularity is not due to numerical error, we used various integration methods and various time steps in the Brian2 simulator (Stimberg, Brette, & Goodman, 2019).

Figure 9:

Firing patterns of the CAdEx model: (a) adaptive spiking, (b) delayed low-frequency tonic spiking, (c) bursting, (d) delayed bursting, (e) accelerated spiking, (f) chaotic-like spiking (see Table 1 for parameters).

Figure 9:

Firing patterns of the CAdEx model: (a) adaptive spiking, (b) delayed low-frequency tonic spiking, (c) bursting, (d) delayed bursting, (e) accelerated spiking, (f) chaotic-like spiking (see Table 1 for parameters).

6  Adaptation and Irregularity

A convenient way to measure spike frequency adaptation is given by the adaptation index
A=1N-1i=1N-1ISIi+1-ISIiISIi+1+ISIi,
(6.1)
where N is the number of subsequent interspike intervals, ISIi. Note that A(-1,1) and is positive for decelerating spike trains and negative for accelerating spike trains.
To measure the irregularity of spiking, we used the coefficient of variation of ISI,
CVISI=σISIISI,
(6.2)
where ISI and σISI are the mean and the standard deviation of ISIs, respectively. According to this definition, CVISI[0,), where CVISI=0 for regular tonic spiking.

Both subthreshold g¯A and postspike δgA adaptation parameters affect the adaptation index (see Figure 10a). This allows the model to reproduce the wide range of A index values observed in neurons (Allen Brain Institute, 2015).

The irregular spiking (like chaotic-like spiking and bursting) is especially pronounced in the transition zone between slow and fast regular spiking regions (see Figure 10b). On the phase diagram, slow, regular spiking corresponds to a postspike reset occurring on the left side of the right branch of the V-nullcline and, consequently, leading to longer interspike intervals, while fast tonic spiking corresponds to a reset occurring on the right side, leading to fast subsequent spiking. In the transition zone, alternations between resets on the left and right sides of the V-nullcline can lead to highly irregular spiking (see the chaotic-like spiking pattern Figure 9f) and bursting (see Figures 9c and 9d).

Figure 10:

Spike frequency adaptation and variation of interspike intervals in the CAdEx model. (a) Adaptation index A and corresponding firing rate as a function of the maximal subthreshold adaptation g¯A and postspike adaptation δgA. (b) Coefficient of variation of interspike intervals CVISI and corresponding firing rate as a function of reset potential VR and postspike adaptation δgA.

Figure 10:

Spike frequency adaptation and variation of interspike intervals in the CAdEx model. (a) Adaptation index A and corresponding firing rate as a function of the maximal subthreshold adaptation g¯A and postspike adaptation δgA. (b) Coefficient of variation of interspike intervals CVISI and corresponding firing rate as a function of reset potential VR and postspike adaptation δgA.

This irregularity of spiking originates from the postspike reset mechanism, which also occurs in the AdEx model (Touboul & Brette, 2008). The correspondence between parameters of both models is intuitive: δgA corresponds to b in the AdEx model, while the reset potential VR has the same meaning in both models.

7  Discussion

In this article, we have proposed a new integrate-and-fire model with two variables, which can produce a large repertoire of electrophysiological patterns while still allowing for clear mathematical insights and large-scale simulations. This CAdEx model is completely specified with 12 biophysical parameters, and reproduces qualitatively similar pattern as the AdEx model (Naud et al., 2008), because the gA nullcline may be considered as locally linear and approximates that of the AdEx model. While the dynamics of the CAdEx model is comparable to the AdEx model for moderate input and firing, the CAdEx model does not suffer from an unnaturally strong hyperpolarization after prolonged periods of strong firing. This can be very advantageous for modeling of highly synchronized rhythms and firings, like slow-wave oscillations or epileptic seizures. The different behaviors of the averaged membrane voltages between the models were demonstrated here at the network level. Moreover, the sigmoidal subthreshold adaptation function allows one to model the dynamics of voltage-dependent ion channels in more detail while retaining the overall computational simplicity. The sigmoidal form of the adaptation function enriches the dynamics, allowing a wider repertoire of multistabilities. An important difference between CAdEx and AdEx models is that the CAdEx model includes a new bifurcation structure with an Andronov-Hopf bifurcation in the four-fixed-point configuration (see Figure 4, upper right). In this regime, there are two simultaneously stable fixed points—one can be in an “integrator” mode (and so going through saddle-node bifurcation when losing stability) with the other one in a “resonator” mode—and the strongly resonating fixed point can then undergo an Andronov-Hopf bifurcation (Izhikevich, 2007). Thus, the CAdEx model can display a bistability between integrator and resonator modes for the same set of parameters. It would be interesting to see if such a predicted bistability can be found in other models or experimentally.

The CAdEx model also solves a problem of divergence of the AdEx model to minus infinity for a negative slope of subthreshold adaptation. This divergence problem can be a source of instabilities in network simulations. It is automatically fixed by the CAdEx model because the voltage is bound to the reversal value of the adaptation current.

Our model also has some limitations. First, the adaptation has the form of a noninactivating current (such as IM potassium current), which limits the description of a class of inactivating ionic channels. It also includes only one type of subthreshold adaptation. In comparison to the AdEx model, the computational cost of our model may be slightly higher due to the form of the adaptation variable and, more specifically, the introduction of an exponential function. Also, as more parameters are introduced, more thorough dynamical studies and explorations of the parameter space are needed.

Appendix

A.1  Shape of V-Nullclines While Varying Input Current

The V-nullcline changes its shape when the input current Is changes. When the input current is strong and hyperpolarizing, the V-nullcline is inverted near an asymptote V=EA (see Figure 11).

A.2  Bifurcation Analysis

A.2.1  Minima and Maxima of the S(V) Function

The derivative of the function S(V) is given by
S'(V)=gL-gLexpV-VTΔT+g¯A1+expVA-VΔA-1ΔAg¯AexpVA-VΔA1+expVA-VΔA2(EA-V).
(A.1)
The solutions of the transcendental equation S'(V)=0 give these locations:
  1. A single global maximum of S(V) function, S(Vmax), for two possible intersections between the V and gA nullclines.

  2. Two maxima and one minimum of S(V) function, for four possible intersections between the V and gA nullclines.

Figure 11:

(a) The V-nullclines while varying input current Is. For strong and hyporpolarizing input, the current V-nullcline is inverted near an asymptote V=EA=-70mV. (b) The four fixed points of the system after inversion of the V-nullcline by a strong, hyperpolarizing current, Is=-200mV.

Figure 11:

(a) The V-nullclines while varying input current Is. For strong and hyporpolarizing input, the current V-nullcline is inverted near an asymptote V=EA=-70mV. (b) The four fixed points of the system after inversion of the V-nullcline by a strong, hyperpolarizing current, Is=-200mV.

A.2.2  Local Linearized Dynamics around Equilibria

The Jacobian of the CAdEx system around equilibrium i, located at (V(i),gA(i)), for an input current Is, has the form
Ji(Is)=-gLC+gLCexpV(i)(Is)-VTΔT-gA(i)(Is)CEA-V(i)(Is)Cg¯AτAΔAexpVA-V(i)(Is)ΔA1+expVA-V(i)(Is)ΔA2-1τA.
(A.2)
The trace trJi(Is) and the determinant detJi(Is) of the Jacobian are
trJi(Is)=-gLC+gLCexpV(i)(Is)-VTΔT-gA(i)(Is)C-1τA,
(A.3)
detJi(Is)=gLτAC-gLτACexpV(i)(Is)-VTΔT+gA(i)(Is)τAC,-g¯A(EA-V(i)(Is))τACΔAexpVA-V(i)(Is)ΔA1+expVA-V(i)(Is)ΔA2.
(A.4)
From the above equations, we obtain:
detJi=-1τAtrJi-1τA2-g¯A(EA-V(i)(Is))τACΔAexpVA-V(i)(Is)ΔA1+expVA-V(i)(Is)ΔA2.
(A.5)
Consider the case with at most two possible fixed points {-,+}. For Is-, V of the left (-) fixed point tends to -. Since gA0, from equation A.2limIs-trJ-(Is)<0. Similarly, from equation A.3, limIs-detJ-(Is)>0. It means that the trajectory of the left fixed point on the Poincaré diagram {trJ-(Is),detJ-(Is)} always starts in the stable quadrant (trJ<0,detJ>0).
The trajectory of the left fixed point ends when it merges with an unstable fixed point on the axis detJ=0. It happens for an input current ISN. If trJ(ISN)<0, the bifurcation is of the saddle-node type. If trJ(ISN)=0, the bifurcation is of the Bogdanov-Takens type. If at the merging point, trJ(ISN)>0, the trajectory crossed the axis trJ=0 previously; therefore, the fixed point lost its stability in an Andronov-Hopf bifurcation. From equation A.4,
trJ(ISN)=-1τA-g¯A(EA-V-(ISN))CΔAexpVA-V-(ISN)ΔA1+expVA-V-(ISN)ΔA2.
(A.6)
The eigenvalues of the Jacobian matrix at fixed points are then given by
λ12=12trJ±(trJ)2-4detJ.
(A.7)
If the eigenvalues are complex, that is, (trJ)2-4detJ<0, then the system oscillates around equilibrium. The imaginary part of an eigenvalue is equal to the angular frequency of the oscillation, Im(λ)=ω=2πν. Consequently, the frequency of oscillations is given by
ν=14π4detJ-(trJ)2.
(A.8)

A.3  Conductance-Based Synapses

To describe the behavior of the system receiving synaptic input, as in section 4 concerning multistability, we used a conductance-based model of synaptic inputs. The synaptic input current to our model is given by,
Isyn=gE(EE-V)+gI(EI-V),
(A.9)
where EE=0mV is the reversal potential of excitatory synapses and EI=-80mV is the reversal potential of inhibitory synapses. gE and gI are, respectively, the excitatory and inhibitory conductances, which increase by quantity QE=4nS and QI=1.5nS for each incoming spike. The increment of conductance is followed by an exponential decrease according to
dgE/Idt=-gE/Iτsyn,
(A.10)
where τsyn=5ms. The spike trains are generated by a Poissonian process with the firing rate modulated by an Ornstein-Ulhenbeck stochastic process (Fourcaud & Brunel, 2002).

A.4  IM and Ih Currents

The currents IM and Ih are proposed as an example of currents that can be simulated using the CAdEx model. For the hyperpolarizing current Ih, the parameters have been identified in Banks et al. (1993), and for IM, the parameters have been identified in Brown and Adams (1980). The values of parameters for these currents are listed in Table 2.

Table 2:
Parameters of the CAdEx to Model Neurons with IM and Ih Currents.
Cm [pF]EA [mV]EL [mV]Is [pA]VA [mV]ΔA [mV]
Neuron with IM 200 −90 −60 350 −35 
Neuron with Ih1111 200 −43 −60 100 −75.7 −5.7 
 VR [mV]1111111 VT [mV] δgA [nS] g¯A [nS] gL [nS] τA [ms] 
Neuron with IM1111 −58 −45 16 10 550 
Neuron with Ih −70 −50 1.5 43 12 800 
Cm [pF]EA [mV]EL [mV]Is [pA]VA [mV]ΔA [mV]
Neuron with IM 200 −90 −60 350 −35 
Neuron with Ih1111 200 −43 −60 100 −75.7 −5.7 
 VR [mV]1111111 VT [mV] δgA [nS] g¯A [nS] gL [nS] τA [ms] 
Neuron with IM1111 −58 −45 16 10 550 
Neuron with Ih −70 −50 1.5 43 12 800 

A.5  Network Simulations

The parameters of the AdEx and CAdEx networks for which results are presented in Figure 1 are shown in Tables 3, 4, and 5. The input to both networks was given by current noise injected to all neurons and had the form of an Ornstein-Uhlenbeck process, with the same parameters: a variance, σ=2.4pA, and a time constant, τ=100ms. The neurons were connected by conductance-based synapses (see appendix A3).

Table 3:
AdEx Network Parameters.
AdEx
Excitatory CellsInhibitory Cells
Population 800 200 
Capacitance C 150 pF 
Leak conductance gL 10 nS 
Leak reversal potential EL −63 mV −65 mV 
Spike threshold VT -50 mV 
Spike initiation slope ΔT 2 mV 0.5 mV 
Adaptation time constant τw 500 mV n.a. 
Post-spike adaptation increment b 107 pA n.a. 
Subthreshold adaptation a 0 nS n.a 
Quantal synaptic conductance QE/I 1.2 nS 5 nS 
Synaptic reversal potential EE/I 0 mV −75 mV 
AdEx
Excitatory CellsInhibitory Cells
Population 800 200 
Capacitance C 150 pF 
Leak conductance gL 10 nS 
Leak reversal potential EL −63 mV −65 mV 
Spike threshold VT -50 mV 
Spike initiation slope ΔT 2 mV 0.5 mV 
Adaptation time constant τw 500 mV n.a. 
Post-spike adaptation increment b 107 pA n.a. 
Subthreshold adaptation a 0 nS n.a 
Quantal synaptic conductance QE/I 1.2 nS 5 nS 
Synaptic reversal potential EE/I 0 mV −75 mV 
Table 4:
CAdEx Network Parameters.
CAdEx
Excitatory CellsInhibitory Cells
Population 800 200 
Capacitance C 150 pF 
Leak conductance gL 10 nS 
Leak reversal potential EL −63 mV −65 mV 
Spike threshold VT -50 mV 
Spike initiation slope ΔT 2 mV 0.5 mV 
Adaptation time constant τA 500 mV n.a. 
Postspike adaptation increment δgA 5 nS n.a. 
Maximal subthreshold adaptation g¯A 0 nS n.a 
Adaptation reversal potential EA −70 mV n.a 
Quantal synaptic conductance QE/I 1.2 nS 5 nS 
Synaptic reversal potential EE/I 0 mV −75 mV 
CAdEx
Excitatory CellsInhibitory Cells
Population 800 200 
Capacitance C 150 pF 
Leak conductance gL 10 nS 
Leak reversal potential EL −63 mV −65 mV 
Spike threshold VT -50 mV 
Spike initiation slope ΔT 2 mV 0.5 mV 
Adaptation time constant τA 500 mV n.a. 
Postspike adaptation increment δgA 5 nS n.a. 
Maximal subthreshold adaptation g¯A 0 nS n.a 
Adaptation reversal potential EA −70 mV n.a 
Quantal synaptic conductance QE/I 1.2 nS 5 nS 
Synaptic reversal potential EE/I 0 mV −75 mV 
Table 5:
Network Connectivity.
Excitatory CellsInhibitory Cells
Excitatory cells 12% 10% 
Inhibitory cells 10% 12% 
Excitatory CellsInhibitory Cells
Excitatory cells 12% 10% 
Inhibitory cells 10% 12% 

A.6  Unbounded Hyperpolarization of the AdEx Model

In the AdEx model, the membrane voltage can be unbounded from below, that is, it tends to minus infinity. This situation can occur if a is negative and gL<|a| (see Figure 12). In this situation, the slope of the oblique asymptote of the V nullcline, which equals gL, is less steep than the slope of the adaptation nullcline. In the shaded area between nullclines in Figure 12, the voltage tends toward minus infinity. This divergence of the AdEx model to minus infinity is automatically fixed by the CAdEx model because the voltage is bound to the reversal value of the adaptation current.
Figure 12:

Unbounded hyperpolarization of the AdEx model. (a) Network simulation. The mean membrane voltage of excitatory AdEx neurons tends to minus infinity. The parameters of the network are the same as in Figure 1b (see Table 3). The only difference is the value of a parameter, here a=-15nS. (b) Corresponding phase plot. Blue line: voltage nullcline; black line: adaptation nullcline; red line: trajectory of the system.

Figure 12:

Unbounded hyperpolarization of the AdEx model. (a) Network simulation. The mean membrane voltage of excitatory AdEx neurons tends to minus infinity. The parameters of the network are the same as in Figure 1b (see Table 3). The only difference is the value of a parameter, here a=-15nS. (b) Corresponding phase plot. Blue line: voltage nullcline; black line: adaptation nullcline; red line: trajectory of the system.

A.7  Mathematical Overview of the AdEx Model

The equations of the adaptive exponential integrate-and-fire (AdEx) model are as follows:
CdVdt=gL(EL-V)+gLΔTexpV-VTΔT-w+Isτwdwdt=a(V-EL)-w,
(A.11)
with the after-spike reset mechanism:
ifVVDthenVVRww+b.
(A.12)
The nullclines of the AdEx model are
w=gLEL-V+gLΔTexpV-VTΔT+Is,w=a(V-EL).
(A.13)
The Jacobian of the AdEx system around equilibrium i, located at (V(i),w(i)) for an input current Is, has the form
Ji(Is)=-gLC+gLCexpV(i)(Is)-VTΔT-1Caτw-1τw.
(A.14)
The trace trJi(Is) and the determinant detJi(Is) of the Jacobian are
trJi(Is)=-gLC+gLCexpV(i)(Is)-VTΔT-1τw,detJi(Is)=gLCτw(1-expV(i)(Is)-VTΔT)+aCτw.
(A.15)
From the above equations, we get
detJi(Is)=-1τwtrJi(Is)-1τw2+aCτw,
(A.16)
which tells us that the relationship between detJ and trJ is linear in the AdEx model, that is, the trajectories of the equilibria on the Poincaré diagram are straight lines.

Acknowledgments

The work was supported by CNRS, the European Community (Human Brain Project, H2020-785907), and the ICODE excellence network.

References

Allen Brain Institute
. (
2015
).
Allen cell types database.
https://celltypes.brain-map.org/data
Banks
,
M. I.
,
Pearce
,
R. A.
, &
Smith
,
P. H.
(
1993
).
Hyperpolarization-activated cation current (IH) in neurons of the medial nucleus of the trapezoid body: Voltage-clamp analysis and enhancement by norepinephrine and camp suggest a modulatory mechanism in the auditory brain stem
.
Journal of Neurophysiology
,
70
(
4
),
1420
1432
.
Benda
,
J.
, &
Herz
,
A. V. M.
(
2003
).
A universal model for spike-frequency adaptation
.
Neural Computation
,
15
(
11
),
2523
2564
. doi:10.1162/089976603322385063
Brette
,
R.
, &
Gerstner
,
W.
(
2005
).
Adaptive exponential integrate-and-fire model as an effective description of neuronal activity
.
Journal of Neurophysiology
,
94
(
5
),
3637
3642
.
Brown
,
D.
, &
Adams
,
P.
(
1980
).
Muscarinic suppression of a novel voltage-sensitive k+ current in a vertebrate neurone
.
Nature
,
283
(
5748
), 673.
Clopath
,
C.
,
Jolivet
,
R.
,
Rauch
,
A.
,
Lüscher
,
H.-R.
, &
Gerstner
,
W.
(
2007
).
Predicting neuronal activity with simple models of the threshold type: Adaptive exponential integrate-and-fire model with two compartments
.
Neurocomputing
,
70
(
10–12
),
1668
1673
.
Connor
,
J. A.
, &
Stevens
,
C. F.
(
1971
).
Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma
.
J. Physiol. (Lond.)
,
213
(
1
),
31
53
.
Fitzhugh
,
R.
(
1961
).
Impulses and physiological states in theoretical models of nerve membrane
.
Biophys. J.
, (
6
),
445
466
.
Fourcaud
,
N.
, &
Brunel
,
N.
(
2002
).
Dynamics of the firing probability of noisy integrate-and-fire neurons
.
Neural Computation
,
14
(
9
),
2057
2110
.
Izhikevich
,
E. M.
(
2003
).
Simple model of spiking neurons
.
IEEE Transactions on Neural Networks
,
14
(
6
),
1569
1572
.
Izhikevich
,
E.
(
2007
).
Dynamical systems in neuroscience: The geometry of excitability and bursting
.
Cambridge, MA
:
MIT Press
.
Krinskii
,
V. I.
, &
Kokoz
,
I. u. M.
(
1973
).
Analysis of equations of excitable membranes. I. Reduction of the Hodgkin-Huxley equations to a second-order system
.
Biofizika
,
18
(
3
),
506
511
.
McCormick
,
D. A.
, &
Contreras
,
D.
(
2001
).
On the cellular and network bases of epileptic seizures
.
Annu. Rev. Physiol.
,
63
,
815
846
.
Morris
,
C.
, &
Lecar
,
H.
(
1981
).
Voltage oscillations in the barnacle giant muscle fiber
.
Biophys. J.
,
35
(
1
),
193
213
.
Naud
,
R.
,
Marcille
,
N.
,
Clopath
,
C.
, &
Gerstner
,
W.
(
2008
).
Firing patterns in the adaptive exponential integrate-and-fire model
.
Biological Cybernetics
,
99
(
4–5
), 335.
Stimberg
,
M.
,
Brette
,
R.
, &
Goodman
,
D. F.
(
2019
).
Brian 2, an intuitive and efficient neural simulator
.
eLife
,
8
, e47314. doi:10.7554/eLife.47314
Touboul
,
J.
, &
Brette
,
R.
(
2008
).
Dynamics and bifurcations of the adaptive exponential integrate-and-fire model
.
Biol. Cybern.
,
99
(
4–5
),
319
334
.
Treves
,
A.
(
1993
).
Mean-field analysis of neuronal spike dynamics
.
Network: Computation in Neural Systems
,
4
(
3
),
259
284
. doi:10.1088/0954-898X_4_3_002

Author notes

T.G. and D.D. contributed equally.