Mode locking in a periodically forced "ghostbursting" neuron model

We study a minimal integrate-and-fire based model of a "ghostbursting" neuron under periodic stimulation. These neurons are involved in sensory processing in weakly electric fish. There exist regions in parameter space in which the model neuron is mode-locked to the stimulation. We analyse this locked behavior and examine the bifurcations that define the boundaries of these regions. Due to the discontinuous nature of the flow, some of these bifurcations are nonsmooth. This exact analysis is in excellent agreement with numerical simulations, and can be used to understand the response of such a model neuron to biologically realistic input.


Introduction
"Ghostbursting" is a novel type of neuronal bursting that has recently been characterised [Doiron et al., 2002]. It occurs in pyramidal cells in the electrosensory system of the weakly electric fish Apteronotus leptorhynchus. These fish communicate, navigate and find prey using their electrosense. They generate a weak electric field with their electric organ (located in the tail) and sinusoidally modulate this in time at frequencies between approximately 800 and 1200 Hz. This electric field then permeates their local environment. Modulations in the amplitude and phase of the oscillation caused by nearby objects are detected by the fish's electroreceptors, which cover the fish's skin. These signals are then passed on to the electrosensory lateral line lobe (ELL) of the fish, and from there to higher brain centers.
In vitro studies of pyramidal cells from the ELL of A. leptorhynchus show that if a constant current (chosen from some appropriate interval) is injected into a cell, the cell's action potentials are clustered together in bursts [Lemon & Turner, 2000]. These bursts are generated in an approximately periodic fashion, provided the injected current is kept fixed. Within a burst, the times between successive action potentials (spikes) monotonically decrease, i.e. the action potentials "speed up." The end of a burst is marked by a very short interspike interval (ISI), which is followed by a long ISI, before the next burst commences. An example is shown in Fig. 1. This type of bursting (dubbed "ghostbursting") was recently reproduced in a biophysically detailed model neuron [Doiron et al., 2001] and further analysed in a simpler yet still realistic model [Doiron et al., 2002].
One novel aspect of this type of bursting is that it seems necessary for the neuron undergoing bursting to be spatially extended, i.e. ghostbursting cannot be seen in a "point" neuron. This is apparent once the bursting mechanism is understood. In summary, both the pyramidal cell soma and dendrite are capable of producing spikes. Just after all but the last somatic spike in a burst, the dendrite fires a spike. During a burst there is a slow decrease in the dendritic potassium inactivation variable, causing the dendritic spikes to become wider in time, and causing the slow speed up of spikes -see Fig. 1. This speed up continues until the dendrite can no longer respond to a somatic spike (the dendritic refractory period is longer than that of the somatic). The lack of a dendritic spike then causes a long gap between two somatic spikes, during which the dendritic potassium inactivation variable increases, and another burst starts. See Doiron et al. [2002] for more details.
The model of Doiron et al. [2002] involves six nonlinear ODEs, which reproduce qualitatively, and to a large extent quantitatively, the experimentally observed bursting [Lemon & Turner, 2000]. Laing & Longtin [2002] introduced a minimal model of a ghostbursting neuron that had two variables, one "voltage" variable whose dynamics were governed by a modified integrate-and-fire mechanism, and a second variable that modulated the spike patterning so as to qualitatively reproduce ghostbursting. It is that model that we study here, under periodic modulation of its input current.
Investigating the effects of periodic forcing on such a model sensory neuron is important for the following reasons. The frequency of electric organ discharge (EOD) for a particular fish is essentially constant. Males have higher frequencies than females, and when two fish are physically close the difference in EOD frequencies gives rise to a beat frequency, equal to that difference. Thus there are at least two types of periodic signals (the fish's own EOD and the beat signal) that could be detected by electroreceptors on the fish's skin and passed directly to the pyramidal cells that we are modelling. It is well known that weakly electric fish can detect beat frequencies, and change their own EOD frequency so as to increase the beat frequency [Zakon et al. 2002]. During such a jamming avoidance response the fish with the lower initial EOD frequency lowers its frequency while the other does the opposite, thereby increasing the beat rate to values that have little effect on electrolocation. This work can also be seen as a generalisation of work on periodically forced oscillators [Glass et al., 1984, Glass & Mackey, 1988.
Some work on the periodically forced ghostburster has been done previously [Laing & Longtin, 2002, Laing & Longtin, 2003. Laing & Longtin [2002] studied the same model that we study here, but their analysis only allowed one to study mode-locked solutions in which the neuron fired once during a whole number of forcing periods, and they also ignored non-smooth bifurcations, whose study is essential for a complete understanding. Here we perform a full analysis of arbitrary mode-locked states for arbitrary periodic forcing, and investigate non-smooth bifurcations of the system. Laing & Longtin [2003] studied the periodically forced neuron model consisting of six nonlinear ODEs presented by Doiron et al. [2002], but did not analyse mode-locked solutions in any detail.
The work presented here is similar to that in Coombes et al. [2001], where the authors analysed a periodically forced integrate-and-fire-or-burst neuron model. That model had a second variable to model a Ca 2 ion current, and was capable of bursting under appropriate stimuli.
The structure of the rest of the paper is as follows. In Sec. 2 we present the discontinuous differential equations describing the minimal ghostbursting neuron and derive an implicit map for successive firing times. In Sec. 3 we describe mode-locked solutions and formulate a set of nonlinear simultaneous equations that any mode-locked solution must satisfy. Sec. 4 discusses the stability of these solutions and Sec. 5 shows some numerical results, comparing direct simulations with numerically determined resonance tongue boundaries. We summarize our main results in Sec. 6.

Model Equations
The differential equations describing our model neuron are where H is the Heaviside step function and ¦ is the Dirac delta function. A(t) is a periodic function of time, with period ∆.
This model was first introduced by Laing & Longtin [2002], although a redundant parameter in that formulation has been scaled out here. Equations (1)-(2) define a modified integrate-and-fire (IF) neuron subject to periodic forcing. Periodically forced IF neuron models have been studied previously [Keener et al., 1981, Coombes, 1999, Coombes & Bressloff, 1999, as have IF models with periodically varying [Glass, 1991] and dynamic [Chacron et al., 2003] thresholds.
The times t n ¡ are referred to as the firing times of the neuron. Between firing times, c decays exponentially to zero with time constant Assuming that the last ISI (t n ¢ t n ¥ 1 ) was greater than r, V is incremented by the value of c at a time t n £ § . If the last ISI was less than r, V is not incremented (this mimics the failure of the dendrite to respond to a somatic action potential). Typical "bursting" behavior of the system in shown in Fig. 2. The current is given by A(t) Note that when I 1 ¡ 0, the system progressively moves from bursting to periodic firing and then quiescence as I 0 is decreased, in agreement with experimental observations [Lemon & Turner, 2000] and the behavior of more sophisticated models [Doiron et al., 2001, Doiron et al., 2002. For more analysis of the model (1)-(2), see Laing & Longtin [2002].
Because (1)-(2) is piecewise linear, we can integrate over t n t t n 1 to obtain The specific case we will study from now on is A(t) The next firing time, t n 1 , is defined implicitly by There are two Heaviside functions in (3), one resulting from the "hard" cutoff in (1) -meant to represent the failure of the neuron's dendrite to respond if the previous ISI was shorter than the dendritic refractory period, r -and one resulting from the integration of the delta function in (1). The presence of these Heaviside functions can give rise to discontinuities in the solution, (3). We prefer to work with a smooth flow, so we can determine firing times by threshold crossings. We achieve this by replacing both Heaviside functions in (3) by Together, (4) and (5) provide a means for calculating the sequence of firing times t n ¡ , provided initial conditions are specified. We now discuss the existence of mode-locked solutions.

Mode-locked Solutions
Mode-locked solutions are those for which the model neuron fires p times for every q periods of the forcing function, where p © q ¤ ¦ ¥ . We refer to these solutions as p : q mode-locked. For such solutions, the firing times are given by ) are a collection of firing phases, and [ ] denotes the integer part. Let us also define p £ 2 firing times (These are just p £ 2 consecutive firing times, shifted so that all but the first and the last lie in the interval (0© q∆).) Substituting the firing times, (7), into (4) we see that the p firing phases are determined by the p equations The c n are given by (5)-(6), and c 0 satisfies Thus we can find p : q mode-locked solutions by simultaneously solving the p equations (11) and the single equation (12).
For periodically forced systems, there typically exist regions in parameter space in which mode-locked solutions exist, i.e. they are codimension-zero phenomena. These regions are called Arnol'd tongues and are often wedge-shaped, with one point of the wedge touching a curve of zero forcing amplitude [Coombes & Bressloff, 1999, Glass, 1991. The boundaries of these tongues are defined by instabilities of the relevant mode-locked orbit. By mapping out these boundaries we can attempt to partition parameter space in terms of the qualitative behavior of the system.

Stability
n about a modelocked solution gives a linearised system of equations: where all partial derivatives are evaluated at the mode-locked solution. We can rewrite (15)-(16) as where the coefficients are We have used the fact that G (t) Here A n ¡ A(T n ), and partial derivatives are given explicitly as and The stability of a p : q mode-locked state thus depends on the behavior of the map If all of the eigenvalues of (p) are less than one in magnitude, the corresponding mode-locked solution will be stable. Various bifurcations of the solution can occur when an eigenvalue of (p) passes through the unit circle in the complex plane under variation of a parameter. Specifically, we expect a saddle-node bifurcation when det( where I is the 3 3 identity matrix, a period-doubling bifurcation when det( and a Neimark-Sacker bifurcation when det( provided the rank of (p) is maximal [Wiggins, 1990]. As an example, it can be shown that det( and thus this equation can be used to find saddle-node bifurcations of 1 : q modelocked solutions. It is interesting to note that the same equation would have been obtained if we had ignored perturbations in the c n when deriving n .
This linear stability analysis of the firing map will not detect bifurcations that occur when mode-locked solutions interact with discontinuities of the firing map. Indeed, we expect non-smooth bifurcations in such an integrate-and-fire system whenever a tangential crossing of the firing threshold occurs [Coombes et al., 2001]. We have observed a "grazing" bifurcation due to this phenomenon in the system under study, in which varying a parameter causes a subthreshold local maximum to increase through the firing theshold, leading to a new firing event that occurs at some time earlier than usual. For mode-locked solutions, these non-smooth bifurcations are defined by In principle there could also be a type of grazing bifurcation in which varying a parameter causes the voltage to reach the firing threshold but with dV d t ¡ 0 at that time, rather than dV d t ¡ 0 [Coombes et al., 2001]. This would also cause the solution to be lost in a non-smooth fashion and could be detected by appending the condition V (t n ) ¡ 0 to the p £ 1 equations defining the mode-locked orbit (11)-(12). However, we have not observed this type of bifurction in the system under study, so any subsequent discussion of grazing bifurcations will refer to the type described in the previous paragraph.
By considering both smooth bifurcations of the firing map and non-smooth bifurcations of the underlying discontinuous flow, we can construct the Arnol'd tongue structure of the periodically forced ghostbursting model neuron. For instabilities of the firing map (saddle-node, period-doubling and Neimark-Sacker bifurcations) we need to solve p £ 2 nonlinear algebraic equations (p £ 1 defining the mode-locked solution and one of (32)-(34)). For the grazing bifurcation described above we need to solve p £ 3 equations, the extra two being (35) and (36). For this type of bifurcation we need to take care with specifying between which two existing firing times the graze occurs.
In practice, numerical integration of (1)-(2) is used to determine approximate solutions of the simultaneous equations we wish to solve and these are used as initial guesses for Newton's method. Once Newton's method has converged we continue the solution in the appropriate system parameters using pseudoarclength continuation [Doedel et al., 1991].

Numerical Results
In this section we show some numerical results, both from direct numerical simulation of (1)-(2) and from numerical continuation of solutions of the simultaneous equations defining the edges of various Arnol'd tongues. Figures 3, 6, 8 and 10 show the most positive Lyapunov exponent associated with the linearised map (17). This is one way of indicating the asymptotic behavior of the system at a particular point in parameter space. The maximal exponent is defined in the usual way: where v is randomly chosen vector in ¦ 3 and (n) is given in (31). At each point in these Figures, only one vector v was chosen to calculate , i.e. we did not average several values of calculated using (38) with different v's. In some regions of these Figures there are "speckles" -this reflects multistability of the underlying dynamical system in these regions. At each point in the parameter plane the system (1)-(2) was integrated for 300 time units, and the first 10 firing times were ignored as transients. Figure 3 shows the maximal Lyapunov exponent as a function of both I 0 and I 1 for ¤ ¡ 3. Other parameters are as in the caption of Fig. 2. We see that when I 1 ¡ 0, i.e. there is no temporal modulation of the input current, the model moves from periodic ( ¡ 0) to chaotic ( ¡ 0) behavior as I 0 increases through approximately 1 2 2. In fact, the model moves from periodic firing to bursting, which is known to be chaotic [Laing & Longtin, 2002]. It can also be seen from Fig. 3 that when I 1 ¡ 0 and I 0 is greater than about 1 38, the system behaves periodically. In this parameter region the model neuron is firing "doublets" -alternating long and short ISIs -behavior seen in other models of this type of neuron [Doiron et al., 2002, Laing & Longtin, 2002. Also, when I 1 ¡ 0, the frequency of firing for I 0 1 2 2 increases as I 0 is increased. Thus we see Arnol'd tongues emanating from the I 0 axis, the largest being the 1 : 1 tongue.
In Fig. 4 we see boundaries of several tongues, calculated using the ideas presented in Sec. 4. Here, as in other Figures, the Arnol'd tongues are constructed from the union of smooth bifurcations of the firing map and non-smooth bifurcations of the underlying discontinuous flow. Neither Neimark-Sacker nor period-doubling bifurcations seem to play a large role in the determination of the tongue structure, although both are present. For example, there is period-doubling within the 2 : 1 tongue shown in Fig. 4, and a horizontal slice through this Figure at I 0 ¡ 1 2 2 is shown in Fig. 5. The first few period-doublings can be seen, which strongly suggests that the system may pass through a period-doubling cascade to chaotic behavior as parameters are varied. Neimark-Sacker bifurcations have previously been found for this model and are discussed by Laing & Longtin [2002]. It is interesting to note that the 2 : 1 border of Fig. 4 has two distinct pieces. Moreover, within the upper 2 : 1 border we find 2(q £ 1) : (q £ 1) solutions that bifurcate in turn from 2q : q solutions as one moves deeper into the interior of the 2 : 1 tongue. In this, and indeed all Arnol'd tongue diagrams, we find excellent agreement between theory and direct numerical simulations.
In Fig. 6 we show the maximal Lyapunov exponent as a function of ¤ and I 1 for I 0 ¡ 1 2 5. From Fig. 3 we see that when I 1 ¡ 0 and I 0 ¡ 1 2 5, the system is chaotically bursting. Thus Fig. 6 shows that in this situation, increasing the amplitude of the forcing can move the system from chaotic behavior to periodic, by forcing it to become mode locked to the forcing. The minimum amplitude required for this depends on the frequency of forcing. In Fig. 7 we plot the corresponding boundaries of several tongues shown in Fig. 6. Figure 8 also shows the maximal Lyapunov exponent in the (¤ © I 1 ) plane, but for I 0 ¡ 1 1 5. For this value of I 0 , the system fires periodically when I 1 ¡ 0, in contrast to the case shown in Fig 6. Here we clearly see Arnol'd tongues emanating from the zero forcing amplitude line. The corresponding theoretical Arnol'd tongue structure is shown in Fig. 9, and nicely illustrates another organizing feature of the Arnol'd tongues; namely that those emanating from the zero forcing amplitude line are arranged such that between a p : q and a p : q tongue one can expect to see a p £ p : q £ q tongue. Figure 10 shows the maximal Lyapunov exponent as a function of both ¤ and I 0 for I 1 ¡ 0 4 . Recalling that when I 1 ¡ 0, the system switches from periodic to chaotic firing at I 0 1 2 2, we see that the nonzero forcing amplitude can either delay the onset of chaotic behavior as I 0 is increased (for example, when ¤ ¡ 5) or hasten it (for example, when ¤ 7). Finally in Fig. 11 we plot the corresponding Arnol'd tongue structure for Fig. 10. For these parameter values we note the extensive overlap of tongue borders, highlighting the multi-stable nature of modelocked solutions, thus accounting for the speckled nature of Fig. 10.

Discussion
We have exactly analysed the response of a model ghostbursting neuron to an arbitrary periodic modulation of its input current, and compared these results with those from numerical integration of the model undergoing sinusoidal forcing. The agreement is very good. Despite the fact that the unforced system undergoes a bifurcation from periodic firing to chaotic bursting as the injected current is increased, its behavior when periodically forced can be largely understood by tracing boundaries of Arnol'd tongues. The boundaries of these tongues correspond to either local bifurcations of the firing time map, or non-smooth bifurcations of the underlying flow, caused by its discontinuous nature.
Unlike many other bursting systems [Izhikevich, 2000] the ghostburster does not have a wide range of timescales, and all ISIs lie in a relatively small interval. We have only examined forcing frequencies for which the period of forcing is of the same order of magnitude as typical ISIs of the unforced system. In biological terms, this would correspond to frequencies around the range 50-500 Hz. Much lower frequencies are presumably also of biological interest, but for these frequencies there would typically be many spikes per period. The practical analysis of mode-locked solutions in this situation would be hampered by the need to solve large numbers of simultaneous nonlinear algebraic equations, and the corresponding Arnol'd tongues would likely be very small. In this case it is likely that the development of an alternative firing rate approach may be more appropriate. An example of "ghostbursting" in the model of Doiron et al. [2002].

List of Figures
The somatic voltage is shown as a function of time. Note the progressive shortening of successive ISIs during each burst, and the large ISIs that separate bursts. The injected current is I ¡ 10; other parameters are as in Doiron et al. [2002]. . . . . . . . . . . . . . . . . . 14 2 Typical behavior of the model (1)