Theta neuron subject to delayed feedback: a prototypical model for self-sustained pulsing

We consider a single theta neuron with delayed self-feedback in the form of a Dirac delta function in time. Because the dynamics of a theta neuron on its own can be solved explicitly -- it is either excitable or shows self-pulsations -- we are able to derive algebraic expressions for existence and stability of the periodic solutions that arise in the presence of feedback. These periodic solutions are characterized by one or more equally spaced pulses per delay interval, and there is an increasing amount of multistability with increasing delay time. We present a complete description of where these self-sustained oscillations can be found in parameter space; in particular, we derive explicit expressions for the loci of their saddle-node bifurcations. We conclude that the theta neuron with delayed self-feedback emerges as a prototypical model: it provides an analytical basis for understanding pulsating dynamics observed in other excitable systems subject to delayed self-coupling.


Introduction
The theta neuron model, or theta neuron for short, is a mathematical model designed to capture essential features of spiking or bursting neurons [1]. It takes the form of a differential equation for a single angular variable θ(t) ∈ (−π, π], representing a phase point moving over an attracting periodic orbit; by convention, the theta neuron produces an output spike or pulse when the angle θ(t) moves through π.
We study here a single theta neuron that receives self-feedback of strength κ after a time delay τ , where we take the feedback to act instantaneously as modeled by a Dirac delta function in time. The overall model takes the form dθ dt = 1 − cos θ + (1 + cos θ) I + κ where · · · t −3 < t −2 < t −1 < t 0 < 0 are the firing times in the past that enter with strength κ and after the delay time τ . System (1) for κ = 0 (without self-feedback) is simply the equation of a single theta neuron with constant input current I. The main advantage of the theta neuron is that its dynamics can be solved explicitly, because its right-hand side is of a particularly simple form. More specifically, it is the (angular) normal form of a saddle-node bifurcation on an invariant-circle (SNIC) bifurcation [2,1], which occurs at I = 0. One often visualizes the theta neuron as a point moving over the unit circle in the plane, as given by the periodic angle θ(t). For I < 0 the system has two equilibria θ ± = ±2 tan −1 ( √ −I), and it is excitable close to the SNIC bifurcation. The point θ − is an attractor, while θ + is unstable (a saddle on the attracting circle) and acts like a threshold. If an initial condition θ(0) (caused by some external perturbation) is smaller than θ + then θ(t) relaxes back to θ − . However, if θ(0) is (slightly) larger than θ + then θ(t) increases through π and, hence, the theta neuron fires before approaching θ − (in the absence of any further perturbations). For I > 0, on the other hand, there are no equilibria and θ(t) increases monotonically, and the theta neuron fires periodically with period π/ √ I. This discussion shows that the theta neuron (without feedback) is an excitable system in the parameter range before it bifurcates to producing sustained periodic oscillations. Excitability is a phenomenon that is observed in numerous natural and man-made systems. It is characterized by an all-or-none response (a pulse or spike) to an external input perturbation, depending on whether or not the perturbation exceeds the so-called excitability threshold. An output pulse corresponds to the sudden release of stored energy and it is followed by a refractory period during which energy is replenished and the excitable system is not able to react to a further purturbation [3]. Excitability has been found experimentally and in associated mathematical models of neurons and other cells, as well as different types of laser systems; see, for example, [4,5] as entry points to the literature.
An excitable unit or cell can react to several input perturbations, where the overall strength of all inputs and their timing determine whether an output pulse ensues. In a (neural) network there are necessarily communication delays between cells that need to be taken into account [6] in order to understand the ensuing dynamics. In the simplest case, a single excitable system receives its own feedback by its output being re-injected after a delay τ . This overall system is able to regenerate its own response after an external input, resulting in feedback-induced self-pulsations whose timing is controlled by the delay time τ . This general mechanism for sustaining pulses in the delay line or delay interval has been demonstrated in a number of laser systems [7,8,9,10,11,12] and, recently also in an experiment with an actual cell [13]. The connection of a neuron to itself is known as an autapse [14] and these are known to occur naturally [15]. Stable self-generation of spikes or pulses are important, for example, as a means of memory storage [16,17].
We study the theta neuron with self-feedback (1) as a prototypical example of such a system, which allows us to determine the observed dynamics explicitly. The underlying idea is that exact details of the excitability are not important when it comes to identifying underlying principles of self-sustained oscillations in the presence of delayed feedback. Figure 1 shows three periodic solutions of Eq. (1) that coexist and are characterized by one, two and three equidistant pulses being regenerated in the feedback loop over one delay interval. It is the existence and stability of these types of solutions that are the focus of this paper.
A number of other authors have studied oscillators with pulsatile interactions [18,19,20], exploiting the fact that it is often possible to explicitly calculate the effects of these interactions. One example is the leaky integrate-and-fire neuron [19], featuring a reset when the firing threshold is reached, which can be solved explicitly. However, that neuron model has the disadvantage of being nonsmooth and phenomenological, whereas the theta neuron is smooth and has been derived from a more complex Hodgkin-Huxley type model via the technique of phase reduction [2,1]. Note further that (1) is exactly equivalent under the transformation V = tan (θ/2) to the quadratic integrate-and-fire (QIF) neuron given by the equation for the voltage V , along with the rule that if V (t − ) = ∞ then V (t + ) = −∞ [21,22]. We will use this equivalence several times below to derive solutions of the theta neuron model. We consider here two important subcases of Eq. (1): the case of I < 0 when the theta neuron is excitable and has positive (excitatory) self-coupling for κ > 0; and the case of I > 0 when the theta neuron is intrinsically firing and receives either excitatory (for κ > 0) or inhibitory (for κ < 0) self-coupling. For ease of notation we define the input current I separately as for these two cases. In Sec. 2 we consider the case of −I 2 m = I < 0 and κ > 0 of the excitable theta neuron with excitatory self-feedback; we derive the existence of periodic orbits, determine their stability analytically, and give an explicit expression for the saddle-node bifurcations of these orbits. Section 3 concerns the case I 2 p = I > 0 of a theta neuron firing periodically even without any feedback. We derive existence and stability of periodic orbits for κ > 0 and then show that these are related to those for κ < 0 via a simple geometric transformation; again, curves of saddle-node bifurcations are given explicitly. In Sec. 4 we consider the theta neuron with a feedback term that is smooth, rather than a delta function. The resulting system needs to be studied numerically, and we find that for both I < 0 and I > 0 the dynamics for excitatory coupling with κ > 0 are qualitatively as those of (1). For I > 0 and κ < 0, on the other hand, we find additional bifurcations that may lead to chaotic behaviour; however, it can be argued that such a model may not be an appropriate description of a neuron operating in this regime of an oscillating neuron with inhibitory self-feedback. We conclude with a discussion in Sec. 5.

Excitable theta neuron for negative current
We first consider the case −I 2 m = I < 0. In order to have nontrivial solutions we consider excitatory coupling only. Our starting point is the derivation of the solution of (1) for κ = 0, using its equivalence to (2). The solution of (2) for −I m < V (0) < I m is Using the transformation V = tan (θ/2) it follows that if −I m < tan (θ(0)/2) < I m , i.e. the neuron's state is initially between the fixed points, then the solution of (1) is and θ approaches the lower fixed point from above as t → ∞. Conversely, if I m < tan (θ(0)/2), i.e. the neuron's state is initially above the upper fixed point, then the solution of (1) for κ = 0 is This means that θ increases through π (the neuron fires) and θ then approaches the lower fixed point from below as t → ∞. While these expressions are general, we will often use them in cases where they simplify. Suppose now that the delta function acts at time t and κ = 0. From (2) and the relationship V = tan (θ/2) we see that this has the effect of instantaneously changing V by the rule where superscripts indicate the state immediately before (−) or after (+) the action of the delta function.
2.1. Existence of periodic solutions for negative current. We now derive the existence of ocillations such as those in Fig. 1. Suppose the neuron fires at time 0 (i.e. θ(0) = π) and there are an additional n past firing times in the interval (−τ, 0) (here n may be zero); in other words, there are n + 1 spikes in the delay interval (−τ, 0]. Assume that the past firing times are evenly spaced with time T between them, meaning that T is the period of oscillation. It follows that τ n + 1 < T < τ n . The time until the neuron feels the next delta function is τ − nT , which is less than T . Since θ(0) = π we have from (6) that The delta function moves θ to the new value: Assuming that tan [θ(τ − nT + )/2] > I m , i.e. we are above the upper fixed point, the phase θ continues to increase and the neuron fires after a further time ∆. Again using (6), this happens when or, equivalently, when We know that the amount of time the neuron waits before the delta function acts, (τ − nT ) plus ∆, has to equal T , i.e. ∆ + τ − nT = T . Thus we have ∆ = (n + 1)T − τ and finally (n + 1)T = τ + 1 Recall that this expression is valid only for τ n+1 < T < τ n . For n = 0 (9) gives the period T explicitly in terms of the other parameters. The periodic solutions shown in Fig. 1 with one up to three equidistant spikes per delay interval of length τ = 20 in panels (a), (b) and (c) correspond to n = 0, 1 and 2, respectively.
Note from (9) that we can rescale parameters to scale one of I m , κ and τ to 1. From now on we set I m = 1 and we investigate the effect of varying τ and κ, the control parameters of the feedback term. For convenience we rewrite (9) as Since nT − τ < 0, we have coth [nT − τ ] < −1, and since (n + 1)T − τ > 0, we have coth [(n + 1)T − τ ] > 1. Thus in order to satisfy (10) we must have κ > 2. Another way to see this is that V has the values ±1 at the equilibria, and κ must be larger than the "gap" between them, which is of size 2.
2.2. Stability of periodic solutions for negative current. We now consider the stability of solutions satisfying (10). As above, assume that the neuron has just fired at time t 0 and there are n past firing times t −1 down to t −n in (t 0 − τ, t 0 ). We wait a time τ − (t 0 − t −n ) until the delta function acts, which maps θ from We then wait a time ∆ until the neuron fires at time t 1 , where Thus . This gives t 1 in terms of previous firing times, but the calculation is general. Hence, we have that is, a map giving t i in terms of previous firing times, which we write as F (t i , t i−1 , . . . , t i−n−1 ) = 0. Note that assuming t i = iT we recover (10).
To calculate stability we perturb t i to t i + η i about a periodic solution. To first order we have where Equation (12) is of the form a i = J a i−1 where a i ∈ R n+1 , so the growth or decay of its solutions (and thus the (in)stability of the periodic solution) is determined by the magnitude of the eigenvalues of the matrix J.
Differentiating (11) we find that ∂F/∂t i = −1, F k = 0 for 2 ≤ k ≤ n, and F n+1 = 1 − F 1 . Now For a periodic solution t i−1 − t i−n−1 = nT , so we can express F 1 in term of T and the n. Since this quantity determines the stability of periodic solutions, we denote it γ from now on for notational convenience, given as We see that γ is always positive: the numerator is clearly positive, and from (10) we have that κ + coth (nT − τ ) = κ − coth (τ − nT ) > 1, so the denominator is also positive. The parameter γ is a function of the other parameters in the model: the feedback parameters κ and τ , the integer n which specifies the form of the solution we consider (with n + 1 equidistant spikes), and also the period T which, in turn, is determined by the values of κ, τ and n.
The matrix J of Eq. (12) is thus To find out more about the eigenvalues of J we consider the determinant of J − λI, which can be expanded down the first column to obtain For even n we have |J − λI| = −λ n+1 + γλ n + 1 − γ, and for odd n we have |J − λI| = λ n+1 − γλ n − 1 + γ. In both cases the condition |J − λI| = 0 can be written as g(λ) = 0 where It is the roots of this polynomial that determine the stability of a periodic solution. We see from (15) that λ = 1 is always a root of g, reflecting the invariance with respect to uniform time translation of the original system. In particular, Note that h(λ) ≡ 1 for n = 0. The properties of the roots of h(λ) were studied in [23]; see also [24,18]. These authors studied a phase oscillator with delayed pulsatile self-feedback but in a model for which, in the absence of feedback, dθ/dt = 1. Their results are not equivalent to ours since in the absence of feedback our model approaches a fixed point. We now collect some useful results.
Proposition 1 (Properties of the eigenvalues at a periodic solution).
These superstable points on branches of solutions corresponding to different values of n all occur at the same value of T . To see this note from (19) that τ −nT = coth −1 (κ/2). Using the oddness of coth and substituting this into (10) we get and, hence, T = 2 coth −1 (κ/2), which is independent of n.
In summary, if γ varies monotonically along a branch of solutions (which we find to always occur), the solution is stable when 0 < γ < 1 and superstable when γ = 1. Moreover, it becomes unstable as γ increases through (n + 1)/n, where an eigenvalue increases through 1; generically, this is a saddle-node bifurcation. Figure 2 illustrates Proposition 1 by showing the typical behaviour of the roots of the function g(λ) in the complex plane as the parameter γ is varied; specifically, for the case n = 4 of a periodic solution with five spikes in the delay interval (−τ, 0]. The value of γ is represented by color in Fig. 2 and panel (a) shows the roots for 0 < γ < 1. For γ = 0 the roots of g(λ) are the fifth roots of unity. As γ increases, the four roots of the factor h(λ) decrease in magnitude and their arguments tend to those of the fourth roots of −1, that is, to π/4, 3π/4, 5π/4 and 7π/4; at γ = 1 there is a quadruple root 0 of h and of g and the periodic orbit is superstable. This behavior can be understood by considering the function h(λ). Since |λ| → 0 as γ → 1, we can approximate h(λ) by λ n + 1 − γ, showing that its roots λ tend to the nth roots of γ − 1, which is negative here. Figure 2(b) shows the behaviour of the roots of g(λ) as γ is increased from 1 in the interval 1 < γ ≤ 2. The four roots of h(λ) leave the origin with arguments close to those of the fourth roots of +1, i.e. 0, π/2, π and 3π/2. This can also be understood with the argument above, since now γ − 1 is positive. Moreover, we see that one real multiplier leaves the unit circle through 1 as γ increases through (n + 1)/n = 5/4. Note that the trivial root 1 of g(λ) remains unchanged throughout in Fig. 2. 2.3. Branches of periodic solutions for negative current. Figure 3 shows solution branches of Eq. (1) with κ = 5 (and I m = 1) for different values of n as indicated. Each branch represents a periodic solution with n + 1 spikes in the delay interval, represented by its period T as a function of τ , with stability indicated as determined in the previous section. The single-spike solution for n = 0 is stable throughout. Branches for n ≥ 1 emerge in pairs, one stable and one unstable, at saddle-node bifurcations as τ is increased; the different stable solutions coexist, leading to an increasing level of multistability. Similar figures for other excitable systems with delayed feedback appear in [13,25,26]. The solution branch for n = 0 is special: it is always stable (when it exists) and the period approaches infinity as τ is decreased sufficiently. This corresponds to a homoclinic bifurcation to the unstable fixed point θ = 2 tan −1 (1) at a finite value τ * . To understand this, imagine that θ(0) = π on this homoclinic orbit. Under backwards time, θ(t) will approach 2 tan −1 (1) as t → −∞. In forwards time, at time t = τ * the impulse from the delayed feedback will precisely make the neuron jump to θ = 2 tan −1 (1), where it will then stay forever. This condition is described by For κ = 5 we obtain the solution τ * = coth −1 (κ − 1) ≈ 0.25541, which is the position of the vertical asymptote in Fig. 3. The minimum in T of this curve for n = 0 occurs at (τ, To determine the branches for n ≥ 1 it is not necessary to solve (10) for additional values of n. Instead, they can be found and plotted by using the reappearance of periodic solutions of DDEs with fixed delay [26], which implies that the branches in Fig. 3 are images of one another under a similarity transformation. Suppose that a DDE (with fixed delay) has a periodic solution with period T 0 for a time delay τ = τ 0 . This same periodic solution is then also a solution for delay τ = τ 0 + mT 0 where m is an integer. Thus, a particular branch of periodic solutions of (1) reappears at higher (and lower) values of the delay. We can express all of these branches parametrically. The branch for n = 0, for which τ < T (τ ), is referred to as the primary branch. It is given explicitly as a function of τ by The branches for n ≥ 1 are referred to as secondary branches. Letting s be a parameter, where coth −1 (κ − 1) < s < ∞, the nth branch is of the form (τ, T ) = (s + nT (s), T (s)).
The reappearance of periodic solutions can also be seen explicitly from (10). Suppose for some n and τ = τ 0 there is a periodic solution with period T 0 , i.e.
Then it follows that where m is any integer; hence, there is also a solution of period T 0 at delay τ = τ 0 + mT 0 with n + m past firing times in the interval (−τ, 0). This reappearance also explains why all superstable points occur at the same value T of T , as described in Sec. 2.2: they are the images of the minimum at (τ, T ) = (T /2, T ) on the primary branch. It follows that the minimum on the n = 1 branch is at (τ, T ) = (3T /2, T ), on the n = 2 branch at (τ, T ) = (5T /2, T ), and on the nth branch at (τ, T ) = ((2n + 1)T /2, T ). Thus, we conclude that the minimum on the nth branch occurs at its (nontransverse) intersection point with the line T = 2τ /(2n + 1). We conclude this section by considering the stability of simultaneously stable periodic solutions. Figure 4 shows the sets of Floquet multipliers for the different stable solutions given by (10) with n = 0, 1, 2, 3, 4 at τ = 4 in Fig. 3. As discussed in Sec. 2.2, these multipliers are roots of the polynomial g(λ) from (15) and can, thus, easily be found numerically for any value of γ (see (13)) which in turn depends on κ and τ . For solutions such as those in Fig. 3, the parameter γ tends to zero from above as one moves along a stable branch away from the saddle-node bifurcation where it is created. According to Proposition 1(i), the respective Floquet multipliers approach the unit circle from within as γ → 0, namely at roots of unity; see also Fig. 2. This explains why the multipliers for smaller n are closer to the unit circle, which is certainly the case for n = 0, 1, 2, 3 in Fig. 4, meaning that the stability of the corresponding periodic solutions is already quite weak. This analytical result for the theta neuron (1) with Dirac delta function is in agreement with, and may serve as an explanation for, the observation in other contexts [10,27] of only weakly stable solutions with n spikes or pulses in the delay interval, whose Floquet multipliers are near roots of unity.

2.4.
Bifurcation curves in the (τ, κ)-plane for negative current. The bifurcation of the branches in Fig. 3 can be continued in the additional parameter κ. Figure 5 shows in the (τ, κ)-plane the curves of homoclinic bifurcation on the primary branch, given by (23), and of saddle-node bifurcations on the secondary branches for n = 1, 2, 3, 4, 5, 6. Generally, such loci of bifurcations need to be continued numerically with standard numerical algorithms [28,29,30]; see also Sec. 4. However, for the theta neuron with delta feedback (1) they can actually be found analytically. Writing the value of τ on the nth branch as τ n = τ 0 + nT (τ 0 ) where τ 0 is the value of τ on the primary branch, a saddle-node bifurcation occurs when dτ n /dτ 0 = 0 [25], i.e. when 0 = 1 + nT (τ 0 ), or T (τ 0 ) = −1/n. Differentiating (24) we have Setting this T (τ 0 ) = −1/n and then solving the resulting quadratic equation for coth τ 0 we obtain coth τ where we have taken the physically meaningful of the two possible solutions. Substituting (26) into (24) we obtain the value of T at which the saddle-node bifurcation on the nth branch occurs as In Fig. 6 the first 100 of the points T (n) for κ = 5 in (27) are plotted on the primary branch. This illustrates that T (n) coverges to the value T at the minimum of the primary branch. The value of τ at which the saddle-node bifurcation of the nth branch occurs is τ (n) = τ (n) 0 + nT (n) ; this expression depends on κ and the curves in Fig. 5 were obtained by plotting τ (n) as a function of κ for various ns indicated.

Oscillating theta neuron for positive current
We now consider (1) with a positive current I = I 2 p > 0, which means that the uncoupled theta neuron is producing regular spikes. We first consider the existence and stability of periodic orbits in the presence of delayed self-coupling. Subsequently, we consider the associated branches of periodic solutions and their bifurcations for the two sub-cases of excitatory selfcoupling for κ > 0, and of inhibitory self-coupling for κ < 0 -showing that they are related by an explicit transformation.
3.1. Existence of periodic solutions for positive current. Our starting point here is again the solution of (2) with κ = 0, which for with I = I 2 p > 0 is With the transformation V = tan (θ/2) we have θ(t) = 2 tan −1 I p tan I p t + tan −1 tan ( θ(0) 2 ) I p .
To find periodic solutions for κ = 0 we assume, as before, that the neuron fires at t = 0 and there are n past firing times in (−τ, 0). The delta function acts at time τ − nT , at which point The delta function acts to move the phase to θ(τ − nT + ) where The neuron fires (i.e. θ reaches π) after another time ∆ where (using (28)) i.e. when which gives Since τ − nT + ∆ = T we have As before, we can rescale either I p , κ or τ to be 1, so set I p = 1 and rewrite (29) as This is the equation relating the period T to the other parameters τ and κ, for a given integer n. As above, for n = 0 expression (30) gives T explicitly as which is valid for 0 ≤ τ ≤ π. Note that, when I = I p = 1, dθ/dt = 2 except at the times at which the feedback acts, which simplifies the derivations below.

3.2.
Stability of periodic solutions for positive current. Determining the stability for I = I 2 p > 0 is also similar to the excitable case. Assume again that the neuron has just fired at time t 0 and there are n past firing times in (−τ, 0). We wait τ − (t 0 − t −n ) until the delta function acts, which maps θ from We then wait a time ∆ until the neuron fires at time t 1 , where π = 2∆ + θ(τ − t 0 + t + −n ). Thus and, in general, This is a map giving the next firing time, t i , in terms of the previous ones, back to t i−n−1 .
Perturbing the firing times defined by this equation as in Sec. 2.2, we obtain the same matrix J as in (14). However, now for a periodic orbit with period T we have which is clearly positive for any value of κ.

3.3.
Branches of periodic orbits for positive current. As in Sec. 2.3, the primary branch of periodic solutions is given by (31) and the secondary branches are given parametrically by (25). To plot and discuss these branches we need to distinguish the sub-cases κ > 0 of excitatory and κ < 0 of inhibitory self-coupling.
3.3.1. The case of excitatory delayed self-coupling. Figure 7 shows the branches of periodic solutions for n = 0, 1, 2, 3, 4, 5 with I = I 2 p = 1 and κ = 2 so that the self-coupling is excitatory. Notice that all branches now connect to form a single curve of periodic solutions; similar plots appear in [31,26]. However, different parts of this curve still correspond to different values of n, as is indicated by colour. The primary branch for n = 0 is again given by a function T (τ ) and it has a finite period throughout for positive current I; that is, there is no longer a homoclinic bifurcation. More specifically, the maximum period now occurs at τ = 0 and is equal to that of an uncoupled oscillator; this is the case since at firing, θ = π and in (1) the term 1 + cos θ is equal to zero, so the feedback can have no effect. The period of this free oscillation is π/I p = π. For the chosen value of κ the primary branch is entirely stable. On the other hand, all secondary branches have two saddle-node bifurcations on them, with an unstable middle sub-branch. Notice also that the secondary branches are increasingly tilted to the right as n increases, leading again to an increasing level of multistability with increasing τ . A periodic orbit is superstable when dT /dτ = 0, and this again occurs at the minima of the period but now also at the maximum period (with τ > 0) where neighbouring branches meet. Differentiating (31) with respect to τ , setting this to zero, and substituting the expression for τ back into (31), we find that the minimum period on the primary branch occurs at (τ, T ) = (T /2, T ) where now T = 2 cot −1 (κ/2); the minima of the secondary branches are therefore at (τ, T ) = ((n + 1/2)T , T ). (For this value of κ we have T = π/2.) The primary branch exists for 0 ≤ τ ≤ π, so the transitions between branches occur when τ is a multiple of π; these are also the locations of the maxima where the periodic orbit is also superstable.
3.3.2. The case of inhibitory delayed self-coupling. Figure 8 shows the shows the branches of periodic solutions for n = 0, 1, 2, 3, 4, 5 with I = I 2 p = 1 and an inhibitory self-coupling with κ = −2. At τ = 0 the period is again π but this is now the minimum period possible; this makes sense since inhibitory feedback can only increase the period. The primary branch is again entirely stable, while the secondary branches each have two points of saddle-node bifurcation that delimit a sub-branch where the periodic orbit is unstable.
The branches in Fig. 8 have been plotted by evaluating (31) and (25) for κ < 0. However, they can also be obtained directly from those in Figs. 7 via the following interesting geometric relationship between the branches of periodic solutions for κ > 0 and for κ < 0.
Solving for τ and T and substituting back into (34), we see that ( τ , T ) satisfies which is exactly of the form (34) but now for κ = −K.
3.4. Bifurcation curves in the (τ, κ)-plane for positive current. To find the saddle-node bifurcations on the nth branch for I = I 2 p > 0 we follow the analysis in Sec. 2.3 of letting τ 0 be the value of τ on the primary branch and solving T (τ 0 ) = −1/n, where T (τ ) is given by (31). This now gives (generically) two values of τ 0 for each n, given by Substituting these into (31) gives the corresponding T -values We now consider first the case that κ is positive. Then the upper sign corresponds to the saddle-node bifurcation with the larger period and the lower sign to that with the smaller period; see Fig. 7. The saddle-node bifurcations with T (n) Note that for the values of τ (n) 0,± to be real we must have κ 2 (n 2 + n) > 1. So for a fixed κ only branches with 1 + 4/κ 2 − 1 2 < n (38) have saddle-node bifurcations on them. For κ = 2, as in Fig. 7, this is satisfied for all n ≥ 1, which means that all secondary branches have a pair of saddle-node bifurcations on them.
(39) Figure 9 shows the curves of saddle-node bifurcations for n = 1, 2, 3, 4, 5 of periodic orbits of (30) with I = I 2 p = 1. The case of excitatory self-feedback is shown in the upper half of the (τ, κ)-plane where κ > 0. Notice that the respective cusp points, where the two curves defined by (37) come together, are minima here. Figure 9 also shows the curves of saddle-node bifurcations for n = 1, 2, 3, 4, 5 for the case of inhibitory coupling, namely in the lower half of the (τ, κ)-plane where κ < 0. These curves can be obtained from Eq. (37) (via (36)), now for negative κ. However, Proposition 2 implies a transformation also of the loci of saddle-node bifurcations on the nth branch of periodic solutions when the sign of κ is changed. Specifically, these two loci are each other's images under rotation by π about the point (τ c , κ c ) = ((n + 1/2)π, 0) .

Theta neuron with smooth self-feedback
We now consider the model of a theta neuron with smooth feedback given by dθ dt where τ is the delay and P (θ) = (8/63)(1 − cos θ) 5 (43) is a pulsatile function centred at θ = π; here the factor 8/63 is a normalisation ensuring that 2π 0 P (θ)dθ = 2π.
The function P (θ) is shown in Fig. 10. Its pulse is quite broad, certainly compared to the Dirac delta function that models the instantaneous delay in system (1). Such forms of smooth coupling have been considered elsewhere [32,33,34], although in infinite networks and without delays. In contrast, the mathematically similar Kuramoto model of coupled phase ocillators with delayed coupling has been well-studied [35,36,37].
We study here the theta neuron with smooth feedback (42)-(43), or smooth theta neuron for short, to investigate the validity of the results we found for the case of instantaneous feedback in system (1). A significant difference between these two systems is that P is a function of θ,  Figure 10. The pulsatile function P (θ) given in (43) as used in the theta neuron with smooth self-feedback (42).
not of time. So if θ increases through π with a speed bounded away from zero, the function P (θ(t)) will be pulse-like in time with a maximum at the time at which θ = π. However, if θ decreases through π (as a result of inhibitory coupling, for example) the neuron will then emit a spurious pulse. The smooth theta neuron (42)-(43) is a delay differential equation (DDE) with a single fixed delay. As such, it is an infinite-dimensional dynamical system whose equilibria and periodic orbits must be expected to undergo standard bifurcations; see, for example, [38,39,40,41]. Periodic solutions of DDEs such as (42)- (43) are not known analytically but must be found with numerical methods [42,43]. To find (stable) periodic solutions we numerically integrate Eq. (42)-(43) with Matlab's dde23 integration routine. We then continue such periodic solutions in a parameter with the software DDE-BIFTOOL [30]. This allows us to compute branches of periodic solutions, regardless of whether they are stable or not, to identify their bifurcations, and to continue the respective bifurcation curves in a two-dimensional parameter space.

4.1.
Excitable smooth theta neuron. We again first consider the case of an excitable smooth theta neuron with I < 0, which means that we must have κ > 0 in (43) to obtain self-sustained solutions. Figure 11 shows typical coexisting stable periodic solutions of the smooth theta neuron (42)- (43). Comparison with those in Fig. 1 reveals that they have similar shapes, with 1, 2 and 3 spikes per delay interval, respectively. However, the influence of the self-feedback is now smooth, rather than acting instantaneously. Note, in particular, that each periodic solution is smooth across the discontinuty in θ, that is, when θ is seen as a point on the unit circle. Figure 12 shows the first four branches of the periodic solutions obtained by numerical continuation as τ is varied; observe that for τ = 4 there are indeed three simultaneously stable solutions as shown in Fig. 11. The branches of periodic orbits are distinct and, as before, can be indexed in Fig. 12 based on the number of times θ increases through π in one delay period; this corresponds to the number n of additional spikes per delay interval. We observe qualitatively the same behaviour as in Fig. 3 for instantaneous self-feedback: the single-spike basic branch for n = 0 is a function over θ, while for n ≥ 1 there is a stable and an unstable branch that meet at a saddle-node bifurcation. Figure 13 shows the bifurcations for n = 0, 1, 2, 3, 4, 5, 6 as curves in the (τ, κ)-plane, which must be found by numerical continuation. The primary branch approaches a homoclinic bifurcation, which was identified, and then continued, as a single-spike orbit of sufficiently large period. This resulted in the curve for n = 0 of Fig. 13, to the left and below of which no periodic solutions exist. The curves for n = 1, 2, 3, 4, 5, 6 are curves of saddle-node bifurcations, and they have been continued from the identified points of saddle-node bifurcations on the secondary branches in Fig. 12. We observe that the bifurcation curves in Fig. 13 are also in perfect qualitative agreement with the corresponding analytically obtained curves in Fig. 5 for instantaneous self-feedback. We conclude that system (1) with Dirac delta function predicts correctly the observed behaviour of the smooth DDE model (42)- (43) for the case that the theta neuron is excitable. This is quite remarkable given that the feedback spike P (θ) in (43) is of a considerable width and not close to a Dirac delta function.

4.2.
Intrinsically oscillating smooth theta neuron. We now consider the smooth theta neuron (42)- (43) with positive I when, in the absence of coupling, the neuron fires periodically with period π/ √ I. When κ is positive the self-coupling is excitatory and numerical continuation started from a stable periodic solution results in the single branch of periodic solutions shown in Fig. 14. While the periodic solution varies smoothly along the branch, successive segments of it can still be associated with an increasing number of firings within a single delay period. We find that the transitions between consecutive numbers of firings take place at the maxima, and we distinguish the respective segments again by colour in Fig. 14 to represent the associated value of additional spikes n = 0, 1, 2, 3, 4, 5. This highlights the very good agreement with Fig. 7 for (1) with Dirac delta function. Excellent qualitative agreement is also observed when we continue the saddle-node bifurcations identified in Fig. 14. The resulting bifurcation curves are shown in Fig. 15 for n = 1, 2, 3, 4, 5, and should be compared with the top half of Fig. 9.
We finally consider the smooth theta neuron (42)- (43) with positive I for the case κ < 0 when the self-coupling is inhibitory. As was expected from the equivalent case of the theta neuron (1) with Dirac delta function in Fig. 8, there is again a single branch of periodic solutions, which is    Compare with Fig. 9 for κ > 0. shown in Fig. 16. More precisely, on this branch one finds orbits with a different number n of additional spikes (this is not shown by colour here), and these all feature a pair of saddle-node bifurcations. We remark that the corresponding curves of saddle-node bifurcations in the (κ, τ )plane (not reproduced here) are qualitatively as those in the bottom half of Fig. 9. However, additional points of bifurcation are now also detected during the continuation of the branch in Fig. 16, specifically, period-doubling and Neimark-Sacker bifurcations. This leads to additional intervals, between pairs of these two types of bifurcations, where the respective periodic orbit is unstable. In particular, there are now ranges of τ for which neither the primary orbit nor the secondary periodic orbit is stable. This opens up the possibility of more complicated periodic and even chaotic dynamics. This is demonstrated in Fig. 17 with three time series near the first period-doubling bifurcation at τ ≈ 2.99, showing periodic, period-doubled and chaotic spiking, respectively. The statement that the solution in panel (c) is chaotic has been checked by verifying that it features a positive Lyapunov exponent.
Overall, we find a more diverse picture for the case that the theta neuron is intrinsically oscillating. For excitatory self-feedback there is still excellent agreement between the theta neuron (1) with Dirac delta function and the smooth theta neuron (42)- (43). For inhibitory self-coupling, on the other hand, we still find a single branch with pairs of saddle-node bifurcations as for system (1), but the smooth DDE (42)-(43) now features additional bifurcations that lead to more complicated dynamics in certain ranges of the delay τ , including chaotic spiking. It might be argued that this more complex behaviour results from the somewhat non-physical nature of the model in this parameter regime of positive current I and negative self-feedback κ, which is why we do not investigate it further here.

Discussion
We have studied the dynamics of a single theta neuron with delayed self-coupling in the form of a delta function. Because the dynamics can be solved explicitly between times at which the feedback occurs, we were able to analytically describe periodic orbits and determine their stability in terms of the roots of a finite-order polynomial. The only possible bifurcations of periodic orbits are saddle-node bifurcations, and we gave explicit expressions for curves in two-dimensional parameter space along which they are found. In this way, we provided a complete description of the types of spiking solutions, where they occur and their stability.
The theta neuron with delayed delta-function self-coupling can be thought of as a "normal form" of an excitable system subject to delayed self-feedback. By this we mean that the system is solvable explicitly, while still capturing the essentials of the behaviour of other excitable systems with non-instantaneous feedback that do not have analytical solutions. Indeed, the kind of dynamics and bifurcation structures presented here have been found (by means of numerical techniques) in different contexts as well [13,7]. To test the predictive power of our results for delta-function feedback, we investigated a single theta neuron with smooth delayed self-coupling. This showed that for excitatory feedback the dynamics are qualitatively the same for both types of intrinsic dynamics (excitable and intrinsically oscillating). In the intrinsically oscillating regime with inhibitory feedback the basic structure of periodic solutions is still as predicted, but we found additional, more complex behaviour generated by period-doubling and Neimark-Sacker bifurcations. These bifurcations have also been found in [27] on branches of periodic spiking in the context of pulse-timing symmetry breaking in a nanolaser with optical feedback.
We now discuss related work that concerns networks of neurons. Several groups have studied infinite networks of QIF (or equivalently, theta) neurons with delayed feedback [21,22,44]. Devalle et al. [21] considered an infinite network of excitable QIF neurons with delayed delta function coupling. The synchronous state in that network is described by the dynamics of one neuron, as we study here. These authors considered the case of a single spike in the delay interval and derived an expression equivalent to (23) giving the minimum value of the feedback strength for which a synchronous periodic solution can exist. They also analysed the stability of such a solution but obtained different results from us, as there are instabilities in infinite networks that cannot occur for a single self-coupled neuron. In similar work, Pazó and Montbrió [22] considered the same network but in the intrinsically firing regime. These authors derived an expression describing the existence of a synchronous state equivalent to (30). They analysed the stability of such a solution and again obtained different results from us, due to the instabilities mentioned above. In relation to the occurence of more complex spiking behaviour dynamics, chaotic dynamics was found, also in [22], in an intrinsically oscillating infinite network of identical QIF neurons with delayed inhibitory delta function feedback. However, this chaotic behaviour required that the neurons are not synchronised, meaning that it is not equivalent to the dynamics of the single-neuron model we studied here. Chaotic dynamics were also found in a network of three intrinsically oscillating theta neurons with nondelayed smooth inhibitory feedback [45], but this is due to the reversibility of the network's dynamics.
Possible generalisations of our work include the study of two coupled neurons [46] or a ring of unidirectionally coupled neurons [18]. For more complex networks, the fact that we can explicitly solve for the dynamics of an uncoupled neuron means that we could efficiently implement event-based simulations: jumping straight from one firing event to the next without having to numerically integrate differential equations between them [47]. Such schemes are very efficient and have been implemented in [48,49], for example. Another possibility is to consider an infinite network of heterogeneous theta neurons for which the Ott/Antonsen ansatz [32,50] could be used to derive a single complex delay differential equation for the network's order parameter, as has been done for delayed Kuramoto oscillators [37,51].