Dynamical systems analysis of spike-adding mechanisms in transient bursts

Transient bursting behaviour of excitable cells, such as neurons, is a common feature observed experimentally, but theoretically, it is not well understood. We analyse a five-dimensional simplified model of after-depolarisation that exhibits transient bursting behaviour when perturbed with a short current injection. Using one-parameter continuation of the perturbed orbit segment formulated as a well-posed boundary value problem, we show that the spike-adding mechanism is a canard-like transition that has a different character from known mechanisms for periodic burst solutions. The biophysical basis of the model gives a natural time-scale separation, which allows us to explain the spike-adding mechanism using geometric singular perturbation theory, but it does not involve actual bifurcations as for periodic bursts. We show that unstable sheets of the critical manifold, formed by saddle equilibria of the system that only exist in a singular limit, are responsible for the spike-adding transition; the transition is organised by the slow flow on the critical manifold near folds of this manifold. Our analysis shows that the orbit segment during the spike-adding transition includes a fast transition between two unstable sheets of the slow manifold that are of saddle type. We also discuss a different parameter regime where the presence of additional saddle equilibria of the full system alters the spike-adding mechanism.

How a single spike or a burst of spikes is generated and regulated for neuron cells is one of the most fundamental questions in neuroscience [1]. Spike generation is closely related to neuronal excitability, which is the ability of the cell’s membrane potential to undergo a large excursion, called an action potential or a spike, when subjected to a sufficiently strong stimulus [1-3]. An excitable cell either responds in full to such a stimulus or not at all, which allows for a reliable transmission of information. Therefore, the different mechanisms for excitability and bursting have been widely studied; we refer to Izhikevich [4] for a comprehensive overview of mechanisms for neurons. Excitable behaviour has also been reported to occur in many other types of cells [1-3] as well as in physical systems, such as lasers [5,6], electronic circuits [7-10] and chemical reactions [11]. Neuronal excitability can be very sensitive to even relatively small changes in, for example, the biophysical properties of a neuron [12-14] or its morphology [15]. This parameter sensitivity indicates that dynamical systems theory is particularly suited for explaining the rich dynamics found in excitable systems.

The classification of different bursting mechanisms was pioneered by Rinzel [16], who used a system decomposition into slow and fast subsystems. He showed that the burst can be divided into active (spiking) and silent phases, which follow different types of attractors of the fast subsystem. Hence, a classification of the bursting oscillators is provided by the structure of the bifurcation diagram of the fast subsystem. Rinzel’s classification was extended in Izhikevich’s studies [1,4]. As an alternative approach, Golubitsky et al. [17] used singularity theory to classify the bifurcation diagrams of the fast subsystem and, hence, the different bursting mechanisms; see also [18]. These studies are primarily for systems with one slow variable; see Smolen et al. [19] for an extension to two slow variables and Ermentrout and Terman [3] for a summary of these ideas along with new results.

A classification of bursting mechanisms, however, does not answer questions about the number of spikes in a particular burst of the same type nor does it explain possible transitions between spiking and bursting. While the latter has received considerable attention over the years, the former has hardly been addressed. Terman [20] analysed transitions between bursting and tonic (continuous) spiking in a pancreatic β-cell model. He recognised the importance of connecting classical slow-fast analysis with full system bifurcation analysis and identified bifurcations of the periodic bursting solutions that organise the transitions between different parameter regimes. He further studied chaotic spiking that can arise in between such transitions [21]. Recently, Benes et al. [22] and Kramer et al. [23] reported that a new type of torus canard can play a role in the transition from spiking to bursting in a model of Purkinje cells. Other recent studies particularly focus on spike-adding in a periodic bursting oscillator; for example, see Govaerts and Dhooge [24], Guckenheimer and Kuehn [25], Tsaneva-Atanasova et al. [26] and Linaro et al. [27]. These studies show that the spike-adding mechanism is formed by a pair of saddle-node bifurcations of periodic orbits of the full system; bursts with different numbers of spikes are, in fact, different periodic attractors of the full system that may coexist only if the number of spikes differs by one [20,26]. The recent study by Teka et al. [28] explains how one may predict the precise number of spikes in a burst; here, the use of two slow variables is essential. Methods to regulate the number of spikes are reported by Ghigliazza and Holmes [29], where a minimal Hodgkin-Huxley model of bursting is proposed to analyse spike-adding and transitions between bursting and tonic spiking in a more general context.

As discussed by Izhikevich [4], the mechanisms for generating spikes do not depend on whether the neuron exhibits spiking only as a transient phenomenon when subjected to a strong enough stimulus or whether it is spiking or bursting continuously. This is indeed the case if transient bursting is organised by the applied stimulus. Applying a stimulus has the effect of changing the right-hand side of the underlying system of ODEs. The bifurcation diagram of the corresponding fast subsystem typically no longer has a stable equilibrium that corresponds to the resting potential, and a spike or burst arises from new attractors that exist only when the stimulus is on. As soon as the stimulus is switched off, the system relaxes back to the resting potential. Therefore, the mechanism is due to a change in the structure of the bifurcation diagram, which depends on the strength of the applied stimulus. Studies of this nature, where the type of burst is studied in dependence on the strength or duration of the stimulus, have been done, for example, by Tran et al. [30], Kim et al. [31] and Stern et al. [32].

In this paper, we investigate spike-adding in a transient burst in the model of hippocampal pyramidal neurons from Nowacki et al. [13]. In contrast to the above studies, the spike-adding occurs after the applied stimulus has been switched off. Hence, the bursting behaviour is governed by the underlying bifurcation structure of the original system, for which a stable equilibrium exists. The strength of the applied stimulus must be such that a first action potential is generated, but the stimulus is not responsible for generating any additional spikes. Our investigation can be compared to excitability in laser systems [5,6] where the response after an applied stimulus is explained by the existence of a nearby homoclinic bifurcation with respect to a parameter. This mechanism is due to the presence of a saddle equilibrium that coexists with the stable equilibrium (i.e., the resting potential). The interesting aspect about our model is that there exists only a stable equilibrium, and we have been unable to identify any nearby saddle-type invariant object in the parameter region of interest that could organise homoclinic bifurcations. We explain the transient behaviour following the ideas of Geometric Singular Perturbation Theory (GSPT) [33-36].

Another important difference with existing studies is that we do not rely on simulations to study the effects of an applied stimulus in numerical experiments. Instead, we use a continuation-based approach and show the process of spike-adding transitions during a transient burst in unprecedented detail, which is not possible using brute-force simulation. Our numerical method is based on the continuation of orbit segments as solutions of a two-point boundary value problem; this approach has already been applied to the bifurcation analysis of periodic orbits, including homoclinic or heteroclinic bifurcations [37], and more recently for the computation of invariant manifolds [38] and so-called slow manifolds in systems with multiple time scales [39,40]. We divide the system into two separate orbit segments, with and without current injection, which are coupled only by the boundary conditions. This allows us to continue the orbit segments in a chosen parameter and analyse the precise nature of their continuous deformation even over exponentially small parameter variations. We complement this continuation with the computation of the two-dimensional critical manifold of the fast subsystem, which comprises all the equilibria of the fast subsystem parametrised by two slow variables. The critical manifold is folded and consists of a number of stable and unstable sheets. In our setting, the unstable sheets are of saddle type, and we will refer to them as saddle-unstable sheets. The fast subsystem also has families of periodic orbits that emanate from Hopf bifurcations on the critical manifold, which give rise to the spiking behaviour during a burst. Our analysis shows that, during the spike-adding transition, orbit segments trace saddle-unstable slow manifolds that lie very close to corresponding saddle-unstable sheets of the critical manifold; the distance between these two manifolds is of the same order as the ratio between the contraction/expansion rates towards and on the manifold, which is organised by the difference between the slow and fast time scales [33,41]. This canard-like behaviour is very similar to behaviour during a spike-adding transition of a periodic burst [25], but it does not involve bifurcations, and coexistence of bursts with different numbers of spikes is not possible here. Thus, our analysis indicates that it is the presence of canard-like behaviour that organises the spike adding.

As for periodic bursting, the spike-adding transition in our model occurs over an exponentially small parameter interval [42-44]. Within this exponentially small parameter interval, we find an even smaller parameter interval during which the canard-like orbit segment includes a fast transition from a saddle-unstable to another saddle-unstable slow manifold. This phenomenon is similar to the so-called fold-initiated canards that have been observed for periodic orbits [45]. To understand this behaviour, we study the associated slow flow on the critical manifold and identify the effect of folds and folded singularities on the behaviour of the orbit segment. Our findings complement the study for planar systems by Guckenheimer et al. [45] and confirm that such fold-initiated canard-like behaviour occurs robustly during a spike-adding transition due to continuity of the vector field with respect to the parameter.

Our study concerns the analysis of a spike-adding mechanism where the full system has a unique stable equilibrium that does not undergo any bifurcations. We find that, for different values of model parameters, the system can have additional unstable equilibria that alter the nature of the spike-adding mechanism. More precisely, the appearance of two saddle equilibria on the critical manifold suppresses the fold-initiated transition between saddle-unstable sheets and changes the behaviour of the orbit segment. The presence of the saddle equilibria give rise to a mechanism that is reminiscent of the excitability studies by Wieczorek et al. [5] and Krauskopf et al. [6]; a new spike is added via a heteroclinic connection of the perturbed state to a saddle equilibrium of the full system. However, this saddle equilibrium lies outside a neighbourhood of the resting potential, and we observe a canard-like transition before and after the heteroclinic connection that is similar to the canard-like transition that occurs when no additional equilibria are present. We compute the critical manifold for this situation and study the associated slow flow to explain this phenomenon; we find that the saddle equilibrium point lies in the middle of a saddle-unstable sheet of the critical manifold, and it is an attractor with respect to the slow flow on the critical manifold.

This paper is organised as follows: In the next section we present our model of the study. In Section 3, we numerically identify the mechanism of spike-adding via continuation of the orbit segment. Next, in Section 4, we calculate the critical manifold of the fast subsystem along with the slow flow to explain the transition to bursting. Here, we also investigate the transition between two unstable slow manifolds of the saddle type and the changes in the spike-adding mechanism when additional equilibria of the full system are present. We end with a discussion in Section 5.

We apply our analysis of transient bursts to pyramidal neuron cells from the CA1 and CA3 regions of the hippocampus. A detailed model of such neurons in Hodgkin-Huxley formalism was presented by Nowacki et al. [13], but for the purpose of this paper, we study a reduced version of this model. The simplified model consists of four ionic currents, namely, fast and slow inward currents, denoted as and , respectively, and fast and slow outward currents, denoted as and , respectively. Inward currents are responsible for the depolarisation or increase of the membrane potential, whereas outward currents hyperpolarise or decrease the membrane potential and return the cell back to its resting state (a stable equilibrium) [2,46]. The fast inward current represents the fastest class of spike-generating - and -currents. The rates of change of these currents are usually similar to that of the membrane potential. Therefore, we assume that the gating of is instantaneous [13,46-48]. The slow inward current mainly corresponds to the transient T-type -current [13,49-51] and represents the low-voltage activated currents responsible for shaping the subthreshold behaviour of the model. The fast outward current represents high-voltage activated fast -currents that we base on the delayed rectifier -current [13,46,47]. Finally, represents muscarinic-sensitive -current [50,52,53], which has an activation rate of the same order as that of .

We only consider the following variables as dynamic: the membrane potential V, the gating variables , and that govern activation of the respective currents and the gating variable for inactivation of , which we denote by . Hence, our reduced model is five dimensional and has the form

Here, is the non-dimensionalised state vector, and is an applied current that stimulates (perturbs) the cell model when it is non-zero. We specifically indicate further parameter dependencies with the parameter vector for some integer . In this paper, we primarily focus on how the system depends on the maximal conductances of and ; these parameters are likely to vary between neurons due to different sizes and numbers of channels in different cells (even among the same types of neurons). The right-hand side of (1) has the specific form that is well known from Hodgkin-Huxley formalism: the dynamics of the membrane potential is organised by the equations for the ionic currents, modelled as

where is the membrane capacitance. Here, with are maximal conductances of the currents, and and are Nernst potentials of the inward and outward currents, respectively. Note that only depends on V, that is, as defined in (4) below. The dynamics of the gating variables is modelled by

the corresponding activation and inactivation steady-state functions of the respective currents, as well as , are given in Boltzmann form as:

Unless specified otherwise, the default values that we use for the parameters of this simplified model are summarised in Table 1.

Figure 1 illustrates three classes of the responses of the simplified model obtained by changing the maximal conductance . These correspond to cell responses that are typically observed experimentally. During the simulations, the model is perturbed from its stable equilibrium by a short current injection whose duration guarantees that the rapidly rising membrane potential will reach and cross its local maximum, creating a fully developed spike; see [12,13,47,48] for more details. Two of the three typical responses shown in Figure 1 exhibit a positive deflection of the membrane potential characterised by a ‘hump’ in the time trace of the membrane potential at the end of the burst; this is called after-depolarisation (ADP), which can exist provided that [13]. Only the first response (lower curve) is a spike without ADP. Note that the last trace, which corresponds to , the highest value of in the example, has sufficiently strong to enable the membrane potential to cross the excitability threshold during the ADP so that additional spikes are fired.

System (1) defined by Equations (2)-(4) evolves on multiple time scales because (as an approximation of the time scale for V) and with have different orders of magnitude. As indicated in Table 1 and are slow variables that vary on a time scale that is (roughly) 10 times slower than and and 100 times slower than V. In particular, this means that our model is capable of firing an arbitrarily large number of spikes during the ADP. More precisely, an increase in , as in Figure 1 and throughout this paper, has the net effect that the slow variable becomes even slower so that more spikes can be fired during the time it takes for to relax back to its equilibrium value. In this paper, we are not interested in the exact nature of this process, but we mention here that a large number of spikes will also be accompanied by a noticeable reduction in oscillation amplitudes. As slows down, the dynamics will resemble more and more the behaviour organised by slow passage through a Hopf bifurcation [26,54]. We focus on the process of spike adding and take advantage of the difference in time scales in Section 3, where it suffices to consider the time-scale separation between the three fast variables V and , and the two slow variables and .

The gating variables express the fractions of channels in a given state and naturally range over the interval . The natural range of the membrane potential V is bounded by the two Nernst potentials [2,46], i.e., , where and . It is beneficial for the numerical analysis if all variables vary over a similar range. Therefore, the computations are done using the scaled membrane potential , where . For our numerical investigations, we used the continuation package AUTO[55,56] for solving the boundary value problems. All numerical simulations were done with XPP [57] using the front-end package XPPy [58] in Python [59], and visualisations were done in Python using Matplotlib [60] and Mayavi [61].

Spike adding happens after a current injection, that is, in the regime where . Hence, any numerical investigation of the transient behaviour must take into account a discontinuous jump from to on the right-hand side of Equation (2). We view the orbit as a concatenation of two orbit segments that are the solution of two boundary value problems and define appropriate boundary conditions to account for the discontinuity in .

More precisely, we consider two successive orbit segments denoted as and , during which current is injected () and during which it is not (), respectively; the concatenation of the two orbit segments and gives the orbit segment that characterises the solution of interest. An illustration of this idea is given in Figure 2, where is the segment coloured red, and is the segment coloured blue. Both and are solution segments of (1) but for different values of and for different integration times and , respectively. Using the set-up that is standard in AUTO[55,56], we formulate a boundary value problem using scaled equations such that the total integration time for both segments is 1. That is, and are solutions of

In order to obtain a unique solution pair , we must impose boundary conditions. The boundary conditions for (5) are determined by the fact that the current injection perturbs system (1) from its resting potential for a fixed duration , as indicated by horizontal black and vertical red dashed lines. Hence, (5) is effectively an initial value problem with , where is the stable equilibrium of (1) with ; we solve for implicitly in AUTO[55,56], and the boundary condition becomes

Equations (5) and (7) uniquely define the orbit segment as a function of λ for fixed . The orbit segment continues on from , but now, . Hence, is again effectively the solution of an initial value problem with initial condition

Throughout this paper, we use for a total duration , which is long enough to drive the system past its threshold for the constants as in Table 1. We fix so that the total integration time of the orbit segment is , which is long enough for to be (approximately) at the resting potential. System (5)-(8) is now well posed and uniquely defines a λ-dependent solution family.

As illustrated by the example in Figure 1, we expect that increasing leads to a spike-adding transition; a new spike is added on top of ADP when it reaches a critical threshold of the membrane potential V. Hence, we set in (5)-(8) and solve it by continuation in AUTO[55,56], starting from . Figure 3a shows the resulting solution branch using the standard -norm of AUTO[55,56] as a measure. We observe that the solution norm exhibits a series of fairly constant ‘plateaus’ that are separated by sharp downward peaks. This behaviour seems similar to that of spike-adding phenomena of periodic bursting solutions, which is organised by pairwise saddle-node bifurcations of periodics [19,26]. However, our numerical set-up imposes a fixed initial condition rather than a periodicity constraint. Hence, the uniqueness of the solutions of (5)-(8) prevents the possibility of coexisting orbit segments, that is, the branch in Figure 3a cannot have folds with respect to . The orbit segments of selected solutions along the branch are shown in Figure 3b,c,d,e,f; note that we present the time series up to for clarity of the presentation.

Figure 3b shows our starting solution, i.e., a single spike followed by ADP. Along the first plateau of the solution branch up to the first downward peak, all orbit segments are qualitatively like Figure 3b; in particular, the ADP is a small hump. As we follow the solutions into the downward peak, the hump of the ADP for the orbit segments stretches out as shown in Figure 3c, which lies at the bottom of the downward peak. Interestingly, as we follow the solution back up along the downward peak, the orbit segment generates a double step in the ADP, as shown in Figure 3d; we selected the orbit segment with the longest double step (with respect to time). As we continue to trace the solution up along the downward peak, the small spike at the end of the orbit segment grows into a fully developed spike, while the stretched double step retracts; the orbit segment shown in Figure 3e is representative of such a solution, and all orbit segments along the second plateau in Figure 3a are qualitatively like Figure 3e. Figure 3f represents a solution along the next plateau, which exhibits three spikes that are created via the same process as explained above for the two-spike burst. In fact, the same spike-adding process takes place for all spike-adding transitions via the downward peaks in Figure 3a. We emphasise here that the stretched single- and double-step ADPs, as shown in Figure 3c,d, only exist along the downward peaks in Figure 3a, that is, in the exponentially small parameter interval during which a spike-adding transition occurs. Hence, such solutions are unlikely to be observed in actual experiments, and they are also very difficult to find in numerical experiments that use initial-value integration methods. The fact that a spike-adding transition happens over an exponentially small parameter interval, during which the solution measure changes rapidly, suggests that it is organised by the difference in time scales present in system (1). Therefore, in order to obtain a better understanding, we proceed by using GSPT [34-36,41,62,63].

As mentioned in Section 2, the full five-dimensional system (1) contains a three-dimensional fast subsystem with variables V, and . Since and are much slower, the idea of GSPT is to assume that and do not change at all and treat them as parameters. More precisely, we consider the singular limit of system (1) and analyse the dynamics of the layer equation

Furthermore, because we are interested in spike-adding phenomena after the brief current injection, we set .

The important objects of study in the singular limit are equilibria and periodic orbits. Since and are parameters, these invariant objects occur in two-parameter families. The -dependent families of equilibria are known as the critical manifold, which we denote by S. The equilibria on S can be stable or unstable, determined with respect to the three-dimensional fast subsystem, and are typically separated by curves of fold or Hopf bifurcations. Similarly, we can expect the existence of -dependent families of periodic orbits that emanate from a curve of Hopf bifurcations on the critical manifold; these periodic orbits can again be stable or unstable with respect to the three-dimensional fast subsystem. Typically, the stable periodic orbits of the fast subsystem organise the spiking phase of the bursting oscillators [16,19,26].

The critical manifold S, when considered in the full five-dimensional phase space of system (1), is a two-dimensional surface, or collection of surfaces, and the associated families of periodic orbits form a three-dimensional manifold, or collection of manifolds, that we denote by P. Together, these collections of manifolds organise the behaviour of solutions of (1). If and vary slowly enough, then GSPT guarantees that a solution of (1) (with ) will trace attracting sheets of S or P that correspond to the -dependent families of attractors of the fast subsystem [41]. For example, the transient spikes of system (1) trace the manifold that corresponds to the family of attracting periodic orbits of the fast subsystem, while and are slowly varying [16,64]. More precisely, solutions of (1) lie on so-called slow manifolds that are perturbations of the different sheets of S and P from the singular limit [41]. Solutions of (1) are characterised by fast transitions between, followed by exponential contraction onto the slow manifolds. The essential difference in behaviour during a spike-adding transition is the fact that the solution of (1) contains a segment that traces a slow manifold associated with a sheet of S that is unstable (of saddle type) rather than attracting; see also [25,34,36,63]. While technically the solutions of (1) trace slow manifolds, we will abuse notation and write ‘sheet of S’ where we mean ‘slow manifold corresponding to the sheet of S.’

The geometry of S and P depends on the values of the other parameters in the system, such as the conductance . In order to illustrate the spike generation, we consider the fast subsystem at the fixed value , which is approximately at the first downward peak in Figure 3a where the solution changes from a one-spike to a two-spike transient burst. Figure 4 shows the critical manifold S for this value of from two different viewpoints; in both views, the embedding into the five-dimensional phase space of (1) is projected onto the -coordinates. The surface was obtained as follows: for ten fixed values of uniformly distributed in the interval , we computed the -dependent curves of equilibria via standard equilibrium continuation with AUTO[55,56], where we allowed to extend outside its physiological range of ; the surface S was obtained via concatenation of this collection of ten -slices, and it is shown in Figure 4 with restricted to the interval for the sake of presentation.

The critical manifold S in Figure 4 forms a single manifold, containing four fold curves, and can be divided into six different sheets depending on the stability type of the equilibria; the stable sheets are coloured black and the unstable ones red. The bottom (black) sheet is labelled , and it contains the resting potential as a stable equilibrium on S that is an actual equilibrium of the full five-dimensional system (1) with . The sheet is connected via a curve of fold bifurcation points to the sheet labelled in Figure 4; this fold curve lies outside and is not shown in Figure 4. The sheet is a two-parameter family of equilibria with two stable and one unstable eigenvalues. Hence, has a four-dimensional stable and a three-dimensional unstable manifold. The next two sheets labelled and have the same stability types as and , respectively; is connected to via the fold curve and and are separated by the fold . Note that the sheet is nearly vertical (with respect to V), as shown in Figure 4; this is not an artefact of the chosen projection. The sheet is connected to via the fold curve , and this sheet consists of equilibria with one stable and two unstable eigenvalues, that is, has a three-dimensional stable and a four-dimensional unstable manifold. The sheet ends at the curve H of Hopf bifurcations, after which it becomes stable again and is labelled . An overview of the different sheets and their stability properties is provided in Table 2.

The maxima and minima of the families of periodic orbits originating from H are shown in Figure 4c,d, using the same two viewpoints as in panels (a) and (b), respectively. The Hopf bifurcation is subcritical along the entire curve so that the emanating family of periodic orbits is unstable (of saddle type), with four-dimensional stable and unstable manifolds; we coloured this family magenta and labelled it . The family of periodic orbits becomes stable via a fold of periodic orbits, after which it is coloured blue and labelled , and ends in a homoclinic bifurcation involving equilibria on the sheet . We refer again to Table 2 for an overview of the different families and their stability properties.

Figure 5 illustrates how orbit segments selected from the first downward peak in Figure 3a trace the different sheets of the critical manifold S for ; these orbit segments are all for virtually the same values that differ only in the seventh decimal point, during which the manifolds S and P hardly change at all. This extreme sensitivity of is a serious challenge for numerical computations, which we overcome by using continuation in AUTO[55,56] of the orbit segments defined by (5)-(8). For clarity of the presentation, we show only the segment , that is, after the current injection, depicted as a blue gradient (cyan to magenta) to visualise its evolution in time. In each panel of Figure 5, the orbit segment starts from , which is located at the top-left in each panel, above the sheet ; since hardly changes, is virtually the same point in all of these panels. The orbit segments traverse the critical manifold S before reaching the stable equilibrium of the full system (1), which lies on the bottom stable sheet . Recall that the full phase space is five dimensional, and Figure 5 may show intersections that are due to projection onto the -coordinates; Figure 5 gives the best possible projection and viewpoint to illustrate the location of the one-dimensional orbit segments relative to S. Some anomalous intersections remain at isolated points, e.g., the intersection with the sheet of S, but the observation that traces sheets of S is real, also in the full five-dimensional space.

The spike-adding process occurs along the downward peak in Figure 3a, during which the parameter remains almost fixed, but the orbit segments of system (5)-(8) change dramatically. We observe the formation of a stretched ADP, which initially gets increasingly longer and shortens again as we follow the orbit segments along the downward peak in Figure 3a. This transition is initiated by the fact that, at the special value , the injected current perturbs the orbit segment such that lies almost on the four-dimensional stable manifold of the saddle-unstable sheet ; more precisely, has a corresponding slow manifold with a corresponding finite-time stable manifold that comes close to. The closer lies to this four-dimensional manifold, the closer the corresponding orbit segment comes to the slow manifold associated with and the longer it will trace this slow manifold. We approximate the slow manifold by the critical manifold , and Figure 5a shows the orbit segment from Figure 3c, which traces almost up to the fold before it drops down to and converges to the resting potential.

Figure 5b,c,d,e,f illustrates orbit segments for the second upward part of the downward peak in Figure 3a. Interestingly, a double-step ADP is created via a transition from to as shown. The orbit that was previously shown in Figure 3d traces both saddle-unstable sheets to of S all the way up to the fold . This transition from to is, in fact, a robust part of the spike-adding process that is a continuous (parameter-dependent) variation of the case where the folds and are absent, that is, where connects directly to ; see Section 4.2 for more details. After reaching the top fold , the membrane potential V initially increases instead of immediately decreasing down to the stable sheet , and a small spike is created. As we continue to follow the solution up along the downward peak, the spike part of the orbit segment grows and moves back towards the attracting periodic orbit family , as illustrated in Figure 5d,e. Finally, as shown in Figure 5f, the orbit segment traces for only a very short time before the second spike occurs; this orbit segment is selected almost at the end of the downward peak, after which orbit segments stop tracing altogether, and the transition from a one- to two-spike transient burst ends. We refer to Govaerts and Dhooge [24], Guckenheimer and Kuehn [25], and Osinga et al. [65] for similar transitions of periodic bursting oscillators.

We remark here that the manner of eventual convergence to the resting potential depends on the nature of the lift-off from the slow manifolds that correspond to the sheet or . Recall that both sheets and of the critical manifold S have three-dimensional stable manifolds; see Table 2. This means that the associated slow manifolds have a one-dimensional repelling fast component [41], and orbit segments that trace these saddle-unstable slow manifolds can leave it only along a single fast direction. We can see this in Figure 5 as a lift-off from ‘down’ in V, shown in Figure 5a, or a lift-off from ‘up’ in V, shown in Figure 5f; this uniquely defined change in direction along the one-dimensional repelling fast component is real and not just due to the projection onto -space. The same holds for the sheet , for which Figure 5b,e is a good example that also shows the required lift-off up from in order to reach . In what follows, the notions up and down are with respect to this uniquely-defined change in direction.

The behaviour of the orbit segment of system (1) in relation to the critical manifold S of the fast subsystem that corresponds to the first downward peak in Figure 3a is representative for what happens along the other downward peaks in Figure 3a. Each time reaches a special value such that the orbit segment comes close enough to the four-dimensional stable manifold of , it gets trapped onto (or, more precisely, the corresponding saddle-unstable slow manifold) for increasingly longer times, and the next spike-adding transition begins. For the parameters of Table 1, we found that this process always includes a transition between two saddle-unstable sheets, which organises the double-step ADP solutions. As mentioned before, the two sheets and together are perturbations of a single sheet connected to that can be obtained continuously via a small parameter variation (using a suitable parameter from Table 1), such that the fold curves and disappear in a curve of cusp bifurcations. Therefore, in the next section, we first explain the jump at the end of the canard-like behaviour, that is, the behaviour near the fold that separates the two saddle-unstable sheets and . We then discuss the transition between and in Section 4.2. Section 4.3 illustrates how the spike-adding mechanism can change when additional equilibria are present.

4.1 Slow flow on the critical manifold near

Let us first focus our attention on the behaviour near the fold , that is, the transition from Figure 5c to Figure 5d. The behaviour near folds can be explained by analysis of the slow flow on the critical manifold S[63]. The slow flow on S is defined by the differential algebraic system

Here, we always have . Recall that the gating variables of (1) are only coupled through the membrane potential V. In fact, it is easy to solve equations and explicitly, which gives us the solutions for the fast gating variables and . We substitute these solutions into to obtain

that is, the slow flow on the two-dimensional critical manifold S is defined by two ordinary differential equations for and and a single algebraic constraint . Unfortunately, S is folded with respect to V so that and do not uniquely define V; however, the algebraic constraint does uniquely define or from given pair or , respectively; compare also Figure 4a,b. Hence, it is advantageous to express the slow flow in terms of only one of the slow variables, or , together with the fast variable V.

We choose to work with V and . If we formally differentiate the algebraic constraint, we obtain

(10)

where is uniquely determined from . We refer to Desroches et al. [63] for more details on this step. Note that (10) becomes singular when , that is, precisely where S has folds with respect to V. We can desingularise the flow by scaling time with the factor . This rescaling reverses the direction of the time whenever , and we obtain the desingularised slow flow in the form

(11)

The actual slow flow on S is now defined by the desingularised slow flow (11) where we must take into account the time reversal in the regimes where . Figure 6 illustrates this for a neighbourhood of the fold on S that separates the sheets and ; we have on and on . The phase portraits in Figure 6 are projected onto the -plane. Figure 6a shows how a trajectory of (11) near (grey line) is attracted to a focus equilibrium of the desingularised slow flow, marked with a black dot on . Figure 6b shows the corresponding projection of the slow flow (9) on S; note the change in direction of the flow for the region where . The fold in Figure 6b is now divided into two parts, a repelling segment on the left side of the focus equilibrium (light-green line) and an attracting segment on the right side of the focus equilibrium (dark-green line). In fact, the focus equilibrium is no longer a focus; it has become a folded singularity or, more precisely, a folded focus. We refer to Wechselberger [36] and Desroches et al. [63] for more details. Figure 7 shows the sheets and of the critical manifold near the folded singularity, with the orbit segment from Figure 5c depicted by a blue-gradient curve as before; panels (a) and (b) provide two opposite viewpoints, and panels (c) and (d), two corresponding close-up views. The slow flow is visualised as a vector field on S, where hotter colours depict vectors with a higher magnitude (the length of the vectors is constant for clarity of presentation). The fold in Figure 7 is coloured the same dark and light green as in Figure 6b.

Figure 7a,b shows a projection onto -coordinates of how the orbit segment follows the slow flow on as it approaches . The first part of the orbit segment lies slightly above (relative to this projection) , and after a jump, another part of this orbit segment lies slightly below . In a neighbourhood of the folded focus on the fold curve , the slow flow has the form of large semi-cycles that cause the orbit segment to trace laterally and, at the same time, push it toward . Since the flow on the top sheet also points towards , as shown in Figure 7c,d, the orbit segment cannot pass and reaches a so-called jump point; compare also with Figure 6b. At the jump point, the fast directions of the flow take over, which causes the formation of a small spike as the orbit segment leaves S; see also Figure 5c. Let us emphasise here that the behaviour of the orbit segments near does not involve interactions with the slow flow on ; the small spike and subsequent drop down to do not intersect the surfaces and . As mentioned earlier in this section, the actual spike formation develops as soon as an orbit segment has reached . After reaching , the orbit segment will lie slightly above (relative to this projection) and experience a lift-off up from (or later only from as shown in Figure 5f). The spike-formation takes place on the fast time scale, and any perceived intersections with and are due to the projection onto -coordinates.

4.2 Slow flow of the critical manifold near the folds and

The formation of a new well-developed spike occurs over an exponentially small parameter interval for which the effect of the injected current is precisely such that the orbit segment comes close to the four-dimensional stable manifold of . The behaviour of the orbit segment near the top fold corresponds to the onset of such a new spike, but the process of reaching , as illustrated in Figure 5a,b, as well as the further development of the spike, as illustrated in Figure 5d,e,f, involves the creation of a double-step ADP; this behaviour is organised by a (fast) jump from to another saddle-unstable sheet . Such a jump, which is actually a jump between the two corresponding saddle-unstable slow manifolds, is a phenomenon that occurs robustly as part of the spike-adding mechanism and has previously been observed for periodic orbits in planar systems; it was reported as a new type of canard called fold-initiated canards in a study by Guckenheimer et al. [45], and a slightly different version termed ducks with relaxation is discussed by Arnol’d [66], Ch.4, Sec.5.4]. In fact, the behaviour we observe in our model is essentially planar and very similar to the example discussed by Guckenheimer et al. [45]. Indeed, there are no folded singularities on the folds and , which means that each fold point has the same effect on the dynamics, and the slow flow is essentially one dimensional. Furthermore, the repelling fast component of the slow manifolds associated with and is one dimensional as well.

The presence of folds and results in the formation of a double-step ADP during the spike-adding process; see Figure 5a,b. The double-step creation is a direct consequence of the fact that the direction of the slow flow on S is transverse to both fold curves and ; this is illustrated in Figure 8. Figure 8a shows the same orbit segment as in Figure 5a, and Figure 8b,c shows two subsequent orbit segments that both occur before the case shown in Figure 5b. As before, the orbit segments are depicted as blue-gradient curves, and the colour-coded vectors indicate the slow flow on S. We observe that the orbit segment in Figure 8a exhibits a lift-off down from before a fast jump down to returns the system to its resting potential, while the two orbit segments in Figure 8b,c exhibit a lift-off up from . These orbit segments are all part of the same continuous one-parameter family of orbit segments that trace ; each orbit segment corresponds to a unique value of even though we always have and the variation is exponentially small, occurring only in the seventh decimal place.

Let us consider this continuous one-parameter family of orbit segments as identified by the moment of lift-off (first down and then up) from . At the start of the spike-adding process, orbit segments trace only the saddle-unstable sheet before a lift-off down to returns the system to its resting potential; the example in Figure 8a shows an orbit segment that almost reaches . As increases continuously (but only exponentially small), the orbit segments come increasingly closer to until one actually reaches ; these orbit segments grow increasingly longer stretched ADPs.

Using the analysis via the desingularised slow flow (11) as derived in Section 4.1, we can decide what happens when an orbit segment reaches . We find that the desingularised slow flow (11) does not have any equilibria in the neighbourhood of the two folds and , which means that there are no folded singularities on either or ; the fold curve is uniformly attracting, which we indicated by a dark-green colour, and is uniformly repelling, indicated by a light-green colour. Hence, upon reaching , the orbit segment simply jumps down toward the resting potential, and subsequent orbit segments exhibit a lift-off up from . Since the sheet on the other side of is attracting, the fast directions will push these orbit segments toward , provided that the lift-off up from occurs not too far away from . Note that the slow flow on points back to , so orbit segments that reach will flow to and again drop down to ; an example is shown in Figure 8b.

As we continue to increase ever so slightly, orbit segments will converge to closer and closer near , until the lift-off up from happens earlier, far enough from , such that they may reach . In order to determine what happens when an orbit segment reaches , we must remember that system (1) and, hence, the slow flow on S depend continuously on . This means that the family of orbit segments is also continuous, and the orbit segment that reaches is a continuous variation of an orbit segment that experiences a lift-off up from later so that it still converges to and flows back to as well as an orbit segment that experiences a lift-off up from earlier so that it misses altogether and forms a well-developed spike as illustrated in Figure 5f. Therefore, the family of orbit segments must include a subfamily of orbit segments that experience a lift-off from . Just as for , the lift-off will first be down as a continuous variation from the orbit segments that flow along and then up to transform into an orbit segment that misses altogether. This subfamily exists in a (doubly) exponentially small parameter regime of fold-initiated canard behaviour [45] and consists of orbit segments that exhibit a double-step ADP.

Another way to understand this phenomenon is in terms of singularity theory, which suggests that and are part of the same surface that unfolds a cusp singularity. Recall that the critical manifold S is a collection of families of equilibria of the fast subsystem of (1); the folds and are curves of saddle-node bifurcations that exist robustly in the -parameter plane organised by the two slow variables. Since and are typically not parallel, they will meet; note that and may need to be extended into an unphysical parameter regime of the -plane. Singularity theory tells us that the two fold curves and typically meet at a cusp point and end there, which means that the two sheets and merge into one in this region of the -plane. The existence of a double-step ADP then merely depends on the location relative to the cusp point of the interaction between the orbit segments and the critical manifold S. We remark that a small change in one or more of the parameters given in Table 1 may move the cusp point into the physical regime or such that and no longer exist for physiologically realistic values of and ; in the latter case, the spike-adding transition will not feature a double-step ADP.

4.3 Spike-adding when additional equilibria are present

It turns out that the spike-adding mechanism organised by canard-like behaviour during the downward peaks of Figure 3a always features a double-step ADP stage involving a jump between and . Hence, each downward peak in Figure 3a corresponds to a qualitatively similar transition as discussed for the first one at . If we increase from the fixed value that was used in Figure 3 to the new value , then the nature of spike adding changes due to the presence of additional unstable equilibria of system (1). If we again continue the two-point boundary value problem (5)-(8) as before with but set to its new value, we get a bifurcation diagram similar to the one for shown in Figure 3. In fact, the spike-adding mechanism for the first four additional spikes involves a double-step ADP stage as we have seen in the previous section. However, for , that is, just before the transition from five to six spikes, a saddle-node bifurcation occurs on the saddle-unstable sheet . This creation of two new (unstable) equilibria prevents a double-step ADP stage; the spike-adding mechanism only involves orbit segments exhibiting a stretched ADP with a singe step, and there is no longer a jump between saddle-unstable slow manifolds.

Let us focus on the transition from a burst with five to one with six spikes, which takes place at . For this value of , there exists three equilibria, but only one is stable so that there is no bistability. The stable equilibrium is the resting potential on . The other two equilibria are saddles, one with one and one with two unstable eigenvalues denoted and , respectively; these additional saddle equilibria are located on . We calculate the critical manifold S of the fast subsystem for ; it is shown in Figure 9 projected onto ()-space. Figure 9 illustrates that the critical manifold does not change qualitatively for higher values of and ; compare with Figure 4. Two orbit segments, one selected from the falling slope and one from the rising slope of the downward peak at , are superimposed onto S; see Figure 9a,b with enlargements in Figure 9c,d, respectively. As before, only the part of the orbit segments that starts after the current injection is shown, so only the downward part of the first of the five spikes is visible. The enlargements in Figure 9c,d also show the slow flow on S in a neighbourhood of the two equilibria and better visualise the interaction of the two orbit segments with and . The equilibria and are both saddles, but with respect to the slow flow on S, the equilibrium is stable (black dot) and is a saddle (red dot).

The value for this case with is again special because at the end of the oscillations, when the orbit segment reaches the family of homoclinic orbits where ends, it lies extremely close to the four-dimensional stable manifold of so that it drops down and traces the saddle-unstable sheet of S. The difference with the spike-adding mechanism illustrated in Figure 5 is that the behaviour of the orbit segment on is affected by the presence of the equilibria and . With respect to the two-dimensional slow flow on , the equilibrium is an attractor, and all orbit segments on converge to provided that they lie in its basin of attraction, which is bounded by the one-dimensional stable manifold of the saddle . In terms of the full five-dimensional flow, is obviously unstable, and orbit segments that come close enough to will behave as dictated by the slow flow for only a finite amount of time; this means that convergence to will eventually be followed by a fast repulsion away from . The orbit segment in Figure 9a,c enters a close enough neighbourhood of in the region of the basin of attraction of ; hence, during the time that it is following the slow flow, it converges to , but we can clearly see in Figure 9c that the fast directions take over before it reaches . Since this orbit segment was selected from the falling slope of the downward peak of the spike-adding mechanism, the orbit segment jumps straight down toward , where it converges to the resting potential. Orbit segments on this slope that lie closer to the minimum of the downward peak would come closer to but still jump down toward when the fast directions take over. On the other hand, orbit segments from the rising slope of the downward peak eventually experience a lift-off ‘up’ from so that a large action potential occurs before converging back to the resting potential; this change of direction corresponds to the onset of a new spike, which is more dramatic and abrupt than the gradual increase in V followed by a small-amplitude spike as illustrated in Figure 5.

Continuity of the vector field (1) implies that there exists an orbit segment that actually converges to the saddle and never relaxes back to the resting potential. This happens when is exactly at the value where the perturbed trajectory lies on the four-dimensional stable manifold of the equilibrium . In contrast to the four-dimensional stable manifold of , the manifold is invariant under the flow of (1), and this heteroclinic connection is a well-defined bifurcation for system (1).

We remark here that the presence of additional equilibria, such as and in the example discussed, only affects the spike-adding mechanism if the orbit segments that trace enter the basin of attraction of . If such orbit segments trace on the other side of the stable manifold of , then a double-step ADP stage would occur. We know from our further model analysis (not shown) that the unstable equilibria persists for higher values of as well as , and in all cases that we investigated, these additional equilibria on affect the spike generation in the way described above.

In this paper, we performed a detailed analysis of the mechanisms of spike generation and spike-adding in a transient burst. Based on a reduction of our previous model [13], we identify these mechanisms using numerical continuation of orbit segments that are solutions to a well-posed boundary value problem. In our analysis, we utilised the separation of time scales in system (1). We calculated the two-dimensional critical manifold S of the fast subsystem, which organises the behaviour of the system. The spike-generation process is characterised by the fact that orbit segments trace saddle-unstable slow manifolds that correspond to saddle-unstable sheets of S. More precisely, there are two saddle-unstable sheets, and , with four-dimensional stable and three-dimensional unstable manifolds; this means that the lift-off from the associated slow manifolds is characterised by a uniquely defined direction. The changes in sign of this direction mark the different phases of the spike-adding transition.

By considering the slow flow on S, we were able to explain the onset of a spike as well as the double-step stretched ADP that leads up to it. For the value of that we considered, the onset of a spike is organised by the top fold of S. This fold contains a folded-focus singularity, but it is not accessible and the spikes are due to (regular) jump points. The folds and that are involved in the double-step ADP do not contain any folded singularities, and they are uniformly attracting and repelling, respectively. Therefore, the first step in the stretched ADP ends at a regular jump point. The second step occurs due to a type of fold-initiated canard because the slow flow points away from ; the fold-initiated canard-like behaviour forms a robust part of the spike-adding mechanism that has also been observed for periodic orbits [45]. The actual spike-adding transition occurs in an exponentially small parameter regime that is very difficult, if not impossible, to find with brute-force integration routines; therefore, it is also highly unlikely to observe anything like a stretched ADP or double-step ADP in experiments.

We found that the nature of the spike-adding mechanism may change if increases slightly. For higher values of , an increase in causes the appearance of two equilibria, and , on that form a saddle-node pair with respect to the slow flow. As it turns out, the presence of these equilibria prevents the double-step ADP. Instead, orbit segments that come close to during the spike-generation process flow towards the attracting equilibrium before a lift-off in the fast direction. This means that the onset of a spike is now organised by rather than the fold , and the stretched ADP involves only a single step. A spike generated by is dramatically different from one generated by . While both spike generations happen in an exponentially small parameter interval, the increase in amplitude of a spike generated by is gradual and should be viewed as a variation of the orbit segment that depends continuously on . On the other hand, a spike generated by is not a continuous variation, and a large-amplitude spike appears abruptly as is increased (on an exponentially small scale). The two families, one with and one without the additional spike, are separated by a heteroclinic connection (via a current injection) from the resting potential to the saddle ; our numerical method for continuation of the family only gets past this discontinuity because we do not impose relaxation back to the resting potential but keep fixed instead. Similar spike-adding mechanisms can be organised by the presence of a saddle periodic orbit or other saddle-type invariant object, but we did not observe this in our model.

In theory, it should be possible to have a double-step ADP as part of the spike-adding transition even when additional equilibria are present. The occurrence of a double-step ADP in this case only depends on whether tracing of commences in the basin of attraction of (restricted to ) or not, which is determined by the stable manifold of the saddle equilibrium . In our numerical explorations, the orbit segments always commence tracing in the basin of attraction of . Hence, we may conclude that and do not have a profound influence on the relative location where orbit segments begin to trace during the spike-adding process. However, other parameters of the system may alter this relative location and provoke a double-step ADP even in the presence of additional equilibria. While this observation indicates a challenge for a precise definition of spike-onset in our context, the two different mathematical notions seem to have the same biological effect.

We believe that the canard-like transition involving saddle-unstable sheets of the critical manifold lies at the heart of any spike-adding mechanism when no invariant saddle-type objects are present. The different phases during the transition, however, could be organised by features other than regular jump points and fold-initiated canards. For example, it should be expected that other folded singularities may appear due to variations in the slow flow on S. An investigation of all possibilities remains an interesting and challenging project for future work.

The authors declare that they have no competing interests.

JN, HO, and KT performed the analysis and wrote or rewrote drafts of the manuscript. JN did all numerical computations and visualisations with assistance from HO and KT. All authors read and approved the final manuscript.

The authors thank Pablo Aguirre, Mathieu Desroches, Bernd Krauskopf, John Guckenheimer and Peter De Maesschalck for their helpful discussions. JN was supported by the grant EP/E032249/1 from the Engineering and Physical Sciences Research Council (EPSRC). HMO was supported by an EPSRC Advanced Research Fellowship grant, and KT-A was supported by the EPSRC grant EP/I018638/1.

1. Izhikevich EM: Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. MIT Press, Cambridge; 2006.

2. Keener JP, Sneyd J: Mathematical Physiology: Cellular Physiology. Springer, New York; 2008.

3. Ermentrout GB, Terman DH: Mathematical Foundations of Neuroscience. Springer, New York; 2010.

4. Izhikevich EM: Neural excitability, spiking and bursting.

   Int J Bifurc Chaos Appl Sci Eng 2000, 10(6):1171-1266. Publisher Full Text

5. Wieczorek S, Krauskopf B, Lenstra D: Multipulse excitability in a semiconductor laser with optical injection.

   Phys Rev Lett 2002, 88(6):1-4.

6. Krauskopf B, Schneider K, Sieber J, Wieczorek S, Wolfrum M: Excitability and self-pulsations near homoclinic bifurcations in semiconductor laser systems.

   Opt Commun 2003, 215(4-6):367-379. Publisher Full Text

7. Nagumo J, Arimoto S: An active pulse transmission line simulating nerve axon.

   Proc IRE 1962, 50(10):2061-2070.

8. Douglas R, Mahowald M, Mead C: Neuromorphic analogue VLSI.

   Annu Rev Neurosci 1995, 18:255-281. PubMed Abstract | Publisher Full Text

9. Guckenheimer JM, Hoffman K, Weckesser W: The forced van der Pol equation I: the slow flow and its bifurcations.

   SIAM J Appl Dyn Syst 2003, 2(1):1-35. Publisher Full Text

10. Indiveri G, Linares-Barranco B, Hamilton TJ, Van Schaik A, Etienne-Cummings R, Delbruck T, Liu SC, Dudek P, Häfliger P, Renaud S, Schemmel J, Cauwenberghs G, Arthur J, Hynna K, Folowosele F, Saighi S, Serrano-Gotarredona T, Wijekoon J, Wang Y, Boahen K: Neuromorphic silicon neuron circuits.

    Front Neurosci 2011, 5:73. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

11. Brøns M, Bar-Eli K: Canard explosion and excitation in a model of the Belousov-Zhabotinskii reaction.

    J Phys Chem 1991, 95(22):8706-8713. Publisher Full Text

12. Brown JT, Randall AD: Activity-dependent depression of the spike after-depolarization generates long-lasting intrinsic plasticity in hippocampal CA3 pyramidal neurons.

    J Physiol 2009, 587(6):1265-1281. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

13. Nowacki J, Osinga HM, Brown JT, Randall AD, Tsaneva-Atanasova K: A unified model of CA1/3 pyramidal cells: an investigation into excitability.

    Prog Biophys Mol Biol 2011, 105(1-2):34-48. PubMed Abstract | Publisher Full Text

14. Brown JT, Chin J, Leiser SC, Pangalos MN, Randall AD: Altered intrinsic neuronal excitability and reduced Na(+) currents in a mouse model of Alzheimer’s disease.

    Neurobiol Aging 2011, 32:2109.e1-2109.e14. Publisher Full Text

15. Van Elburg RAJ, Van Ooyen A: Impact of dendritic size and dendritic topology on burst firing in pyramidal cells.

    PLoS Comput Biol 2010., 6(5):

16. Rinzel J: A formal classification of bursting mechanisms in excitable systems. In Mathematical Topics in Population Biology, Morphogenesis, and Neurosciences. Edited by Teramoto E, Yamaguti M. Springer, Berlin; 1987:267-281.

17. Golubitsky M, Josić K, Kaper T: An unfolding theory approach to bursting in fast-slow systems. In Global Analysis of Dynamical Systems: Festschrift Dedicated to Floris Takens for His 60th Birthday. Edited by Broer H, Krauskopf B, Vegter G. Institute of Physics Publishing, Bristol; 2001:277-308.

18. Osinga HM, Sherman A, Tsaneva-Atanasova K: Cross-currents between biology and mathematics: The codimension of pseudo-plateau bursting.

    Discrete Contin Dyn Syst, Ser A 2012, 32(8):2853-2877.

19. Smolen P, Terman DH, Rinzel J: Properties of a bursting model with two slow inhibitory variables.

    SIAM J Appl Math 1993, 53(3):861-892. Publisher Full Text

20. Terman DH: The transition from bursting to continuous spiking in excitable membrane models.

    J Nonlinear Sci 1992, 2(2):135-182. Publisher Full Text

21. Terman DH: Chaotic spikes arising from a model of bursting in excitable membranes.

    SIAM J Appl Math 1991, 51(5):1418-1450. Publisher Full Text

22. Benes GN, Barry AM, Kaper TJ, Kramer Ma, Burke J: An elementary model of torus canards.

    Chaos 2011., 21(2):

23. Kramer M, Traub R, Kopell N: New dynamics in cerebellar purkinje cells: torus canards.

    Phys Rev Lett 2008, 101(6):68103.

24. Govaerts W, Dhooge A: Bifurcation, bursting and spike generation in a neural model.

    Int J Bifurc Chaos Appl Sci Eng 2002, 12(8):1731-1741. Publisher Full Text

25. Guckenheimer J, Kuehn C: Computing slow manifolds of saddle type.

    SIAM J Appl Dyn Syst 2009, 8(3):854-879. Publisher Full Text

26. Tsaneva-Atanasova K, Osinga HM, Rieß T, Sherman A: Full system bifurcation analysis of endocrine bursting models.

    J Theor Biol 2010, 264(4):1133-1146. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

27. Linaro D, Champneys A, Desroches M, Storace M: Codimension-two homoclinic bifurcations underlying spike adding in the Hindmarsh-Rose burster. SIAM J Appl Dyn Syst, in press. arXiv; 2011. [http://arxiv.org/abs/1109.5689]Linaro D, Champneys A, Desroches M, Storace M: Codimension-two homoclinic bifurcations underlying spike adding in the Hindmarsh-Rose burster. SIAM J Appl Dyn Syst, in press. arXiv; 2011. [http://arxiv.org/abs/1109.5689]

28. Teka W, Tabak J, Vo T, Wechselberger M, Bertram R: The dynamics underlying pseudo-plateau bursting in a pituitary cell model.

    J Math Neurosci 2011., 1:

29. Ghigliazza RM, Holmes PJ: Minimal models of bursting neurons: how multiple currents, conductances, and timescales affect bifurcation diagrams.

    SIAM J Appl Dyn Syst 2004, 3(4):636-670. Publisher Full Text

30. Tran D, Sato D, Yochelis A, Weiss J, Garfinkel A, Qu Z: Bifurcation and chaos in a model of cardiac early afterdepolarizations.

    Phys Rev Lett 2009, 102(25):1-4.

31. Kim MY, Aguilar M, Hodge A, Vigmond E, Shrier A, Glass L: Stochastic and spatial influences on drug-induced bifurcations in cardiac tissue culture.

    Phys Rev Lett 2009, 103(5):1-4.

32. Stern JV, Osinga HM, LeBeau A, Sherman A: Resetting behavior in a model of bursting in secretory pituitary cells: distinguishing plateaus from pseudo-plateaus.

    Bull Math Biol 2008, 70(1):68-88. PubMed Abstract | Publisher Full Text

33. Jones CKRT: Geometric singular perturbation theory. In Dynamical Systems. Springer, Heidelberg; 1995:44-118.

34. Dumortier F: Techniques in the theory of local bifurcations: blow-up, normal forms, nilpotent bifurcations, singular perturbations. In Bifurcations and Periodic Orbits of Vector Fields (Montreal, PQ, 1992). Kluwer Academic, Dordrecht; 1993:19-74.

35. Szmolyan P, Wechselberger M: Canards in .

    J Differ Equ 2001, 177(2):419-453. Publisher Full Text

36. Wechselberger M: Existence and bifurcation of canards in in the case of a folded node.

    SIAM J Appl Dyn Syst 2005, 4(1):101-139. Publisher Full Text

37. Champneys AR, Kuznetsov YA, Sandstede B: A numerical toolbox for homoclinic bifurcation analysis.

    Int J Bifurc Chaos Appl Sci Eng 1996, 6(5):867-887. Publisher Full Text

38. Krauskopf B, Osinga HM: Computing invariant manifolds via the continuation of orbit segments. In Numerical Continuation Methods for Dynamical Systems: Path Following and Boundary Value Problems. Edited by Krauskopf B, Osinga HM, Galán-Vioque J. Springer, Dordrecht; 2007:117-154.

39. Desroches M, Krauskopf B, Osinga HM: Mixed-mode oscillations and slow manifolds in the self-coupled FitzHugh Nagumo system.

    Chaos 2008., 18(1):

40. Desroches M, Krauskopf B, Osinga HM: The geometry of slow manifolds near a folded node.

    SIAM J Appl Dyn Syst 2008, 7(4):1131-1162. Publisher Full Text

41. Fenichel N: Geometric singular perturbation theory for ordinary differential equations.

    J Differ Equ 1979, 31:53-98. Publisher Full Text

42. Benoît É, Callot JL, Diener F, Diener M: Chasse au canard.

    Collect Math 1981, 31-32(1-3):37-119.

43. Dumortier F, Roussarie R: Canard cycles and center manifolds.

    Mem Am Math Soc 1996, 121(577):1-101.

    [With an appendix by Cheng Zhi Li.]

    

44. Lee E, Terman D: Uniqueness and stability of periodic bursting solutions.

    J Differ Equ 1999, 158:48-78. Publisher Full Text

45. Guckenheimer JM, Hoffman K, Weckesser W: Numerical computation of canards.

    Int J Bifurc Chaos Appl Sci Eng 2000, 10(12):2669-2687. Publisher Full Text

46. Hodgkin AL, Huxley AF: A quantitive description of membrane current and its application to conduction and excitation in nerve.

    J Physiol 1952, 105(117):500-544.

47. Golomb D, Yue C, Yaari Y: Contribution of persistent current and M-type current to somatic bursting in CA1 pyramidal cells: combined experimental and modeling study.

    J Neurophysiol 2006, 96(4):1912-1926. PubMed Abstract | Publisher Full Text

48. Yue C, Remy S, Su H, Beck H, Yaari Y: Proximal persistent channels drive spike afterdepolarizations and associated bursting in adult CA1 pyramidal cells.

    J Neurosci 2005, 25(42):9704. PubMed Abstract | Publisher Full Text

49. Jaffe DB, Ross WN, Lisman JE, Lasser-Ross N, Miyakawa H, Johnston D: A model for dendritic accumulation in hippocampal pyramidal neurons based on fluorescence imaging measurements.

    J Neurophysiol 1994, 71(3):1065-1077. PubMed Abstract | Publisher Full Text

50. Yaari Y, Yue C, Su H: Recruitment of apical dendritic T-type channels by backpropagating spikes underlies de novo intrinsic bursting in hippocampal epileptogenesis.

    J Physiol 2007, 580(2):435-450. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

51. Blackmer T, Kuo SP, Bender KJ, Apostolides PF, Trussell LO: Dendritic calcium channels and their activation by synaptic signals in auditory coincidence detector neurons.

    J Neurophysiol 2009, 102(2):1218-1226. PubMed Abstract | Publisher Full Text | PubMed Central Full Text

52. Yue C, Yaari Y: KCNQ/M channels control spike afterdepolarization and burst generation in hippocampal neurons.

    J Neurosci 2004, 24(19):4614-4624. PubMed Abstract | Publisher Full Text

53. Yue C, Yaari Y: Axo-somatic and apical dendritic Kv7/M channels differentially regulate the intrinsic excitability of adult rat CA1 pyramidal cells.

    J Neurophysiol 2006, 95(6):3480-3495. PubMed Abstract | Publisher Full Text

54. Baer SM, Erneux T, Rinzel J: The slow passage through a Hopf bifurcation: delay, memory effects, and resonance.

    SIAM J Appl Math 1989, 49:55-71. Publisher Full Text

55. Doedel EJ: AUTO: a program for the automatic bifurcation analysis of autonomous systems.

    Congr Numer 1981, 30:265-284.

56. Doedel EJ, Oldeman BE: AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations; 2007. [cmvl.cs.concordia.ca/auto/]Doedel EJ, Oldeman BE: AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations; 2007. [cmvl.cs.concordia.ca/auto/]

57. Ermentrout GB: Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. SIAM, Philadelphia; 2002.

58. Nowacki J: XPPy; 2011. [http://seis.bris.ac.uk/~enxjn/xppy]Nowacki J: XPPy; 2011. [http://seis.bris.ac.uk/~enxjn/xppy]

59. Oliphant T: Python for scientific computing.

    Comput Sci Eng 2007, 9(3):10-20.

60. Hunter J: Matplotlib: a 2D graphics environment.

    Comput Sci Eng 2007, 9(3):90-95.

61. Varoquaux G, Ramachandran P: Mayavi: making 3D data visualization reusable. In Proceedings of the 7th Python in Science Conference. SciPy, Pasadena; 2008:51-56.

62. Hek G: Geometric singular perturbation theory in biological practice.

    J Math Biol 2010, 60(3):347-386. PubMed Abstract | Publisher Full Text

63. Desroches M, Guckenheimer JM, Krauskopf B, Kuehn C, Osinga HM, Wechselberger M: Mixed-mode oscillations with multiple time scales.

    SIAM Rev 2012, 54(2):211-288. Publisher Full Text

64. Rinzel J, Ermentrout GB: Analysis of neural excitability and oscillations. In Methods in Neuronal Modelling. Edited by Koch C, Sagev I. MIT Press, Cambridge; 1998:251-292.

65. Osinga HM, Tsaneva-Atanasova KT: Dynamics of plateau bursting depending on the location of its equilibrium.

    J Neuroendocrinol 2010, 22(12):1301-1314. PubMed Abstract | Publisher Full Text

66. Arnol’d VI: Dynamical systems V. In Encyclopaedia of Mathematical Sciences. Springer, Berlin; 1994.