M-current induced Bogdanov–Takens bifurcation and switching of neuron excitability class

In this work, we consider a general conductance-based neuron model with the inclusion of the acetycholine sensitive, M-current. We study bifurcations in the parameter space consisting of the applied current Iapp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$I_{app}$\end{document}, the maximal conductance of the M-current gM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$g_{M}$\end{document} and the conductance of the leak current gL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$g_{L}$\end{document}. We give precise conditions for the model that ensure the existence of a Bogdanov–Takens (BT) point and show that such a point can occur by varying Iapp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$I_{app}$\end{document} and gM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$g_{M}$\end{document}. We discuss the case when the BT point becomes a Bogdanov–Takens–cusp (BTC) point and show that such a point can occur in the three-dimensional parameter space. The results of the bifurcation analysis are applied to different neuronal models and are verified and supplemented by numerical bifurcation diagrams generated using the package MATCONT. We conclude that there is a transition in the neuronal excitability type organised by the BT point and the neuron switches from Class-I to Class-II as conductance of the M-current increases.


Introduction
Neuromodulators are chemicals released by neurons that can alter the behaviour of individual neurons and large populations of neurons. Examples include dopamine, seratonin and acetylcholine. These chemicals occur widely in the brain and can affect many types of neurons. The effect of neuromodulators ranges from altering the membrane properties of individual neurons to altering synaptic transmission.
The M-current is a voltage-dependent, non-inactivating potassium current, which has been shown to occur in many neural types including excitatory neurons in the cortex [1] and inhibitory neurons in the hippocampus [2]. Its name arises from the fact that this current is down-regulated by the presence of the neuromodulator acetylcholine through its action on the muscarinic receptor. At the simplest level, this current reduces firing activity since it is a potassium current [2,3]. However, this current has been implicated in many aspects of both individual cell and network activity.
Before reviewing the literature on the M-current, we first recall some terminology. Neurons and neural models are often classified by their membrane excitability class. As described by Hodgkin [4], neurons with the Class-I excitability have a continuous frequencycurrent (F/I) curve because they begin repetitive firing with zero frequency from the resting state. On the other hand, the frequency-current curve of Class-II neurons is discontinuous because they start firing with non-zero frequency from the resting state [5]. The phase resetting curve (PRC) describes the effect of stimulation on the phase of an oscillator as a function of the phase at which the stimulus is delivered. At phases where the PRC is positive the phase is advanced, meaning the period of the oscillator is increased by the perturbation. At phases where the PRC is negative the phase is delayed, corresponding to a decrease in the period of the oscillator [6,7]. As introduced by Hansel [8], a Type-I PRC is one where an excitatory stimulus produces only phase advances, while in a Type-II PRC either phase advance or phase delay can occur, depending on the phase of the stimulus. For two oscillators with reciprocal excitatory coupling, a Type-I PRC means the coupling cannot synchronise the oscillators, while a Type-I PRC means that the coupling can synchronise the oscillators. For inhibitory coupling, the opposite occurs [6,7]. Another important classification of neurons is whether or not they exhibit subthreshold oscillations. Neurons that do exhibit subthreshold oscillations are called resonators, while neurons that do not are called integrators.
At the single cell level, the M-current has been shown to affect the neuronal excitability [1,9] and resonant properties [10][11][12]. For example, in [1], the authors recorded from layer II/III pyramidal neurons and determined PRCs. Stiefel et al. [1] found that down-regulation of slow voltage-dependent potassium currents such as the M-current can switch the PRC from Type-II to Type-I, thus changing the expected synchronisation of pairs of coupled neurons. In a follow-up paper [1], they showed for that the M-current could produce the same effect in several different neural models. The work of [13] showed that these differences in PRC type due to M-current modulation translate into differences in synchronisation properties in networks of model neurons. The experimental work of [11,12] showed that increased membrane conductance (shunting) could switch a hippocampal pyramidal neuron from an integrator to a resonator. Using a simple model, they attributed this change to the combined effect of shunting (modelled as a leak current) and the M-current. Interpreting the shunt as representing the effect of background synaptic on a neuron, Prescott et al. [12] concluded that neurons that present as integrators in vitro may act as resonators in vivo. At the network level the M-current has also been implicated in the organisation of rhythms in striatal microcircuits. In [14], the authors studied an inhibitory neuron model with M-current under forcing from gamma pulses and a sinusoidal current of theta frequency. They found that the M-current expands the phase-locking frequency range of the network, counteracts the slow theta forcing and admits bistability in some parameter range. In [15], the effects of the M-correct on β oscillations was studied.
In all the studies cited above, the effect of acetylcholine, through the M-current, was explored in models for specific cells. While this is important for understanding the behaviour of specific cells and brain networks, it can be difficult to extract the essential effects of the M-current from its interplay with other specific currents in the models. Here we take a different approach and consider the effect of the M-current in a general conductance-based model. We study the bifurcations of the model in the parameter space of two parameters common to any conductance-based model with an M-current: the applied current and the maximal conductance of the M-current. We derive necessary and sufficient conditions for the existence of two codimension two bifurcations of the resting equilibrium point: the Bogdanov-Takens (BT) bifurcation and the cusp (CP) bifurcation.
The Bogdanov-Takens (BT) bifurcation is associated with an equilibrium point that has a zero eigenvalue with algebraic multiplicity two and geometric multiplicity one. The cusp bifurcation occurs when three equilibrium points coalesce into one, and can be thought of as the simultaneous occurrence of two fold bifurcations. When an equilibrium point simultaneously undergoes BT and cusp bifurcations, a Bogdanov-Takens-cusp (BTC) occurs, which is a codimension three bifurcation. We show that variation of a third parameter, the leak conductance, can lead to a BTC bifurcation point.
In the literature, there are many instances where the presence of the BT, cusp and to a lesser extent the BTC, bifurcations has been shown to occur in particular conductancebased models. For example, the presence of BT and cusp bifurcations [16] and BTC bifurcation [17] has been shown in the Hodgkin-Huxley model. In [18], the author showed the existence of BT and cusp bifurcations in Morris-Lecar model [19]. While in [7] the BT and cusp bifurcations were shown in the Wang-Buzsáki interneuron model [20]. The majority of these studies used numerical bifurcation analysis to show that these bifurcations occur as particular parameters are varied with all other parameters fixed at some specific, biologically relevant values. The prevalence of these codimension two bifurcations in particular studies would seem to indicate that these bifurcations are associated with some underlying structure in conductance-based models in general. Indeed, two recent papers give support to this hypothesis. The authors in [21] considered a general conductancebased neuron model and studied the existence of the BTC point in the parameter space of the applied current, leak conductance and capacitance. In [22], the authors give general conditions for the existence of the BT bifurcation in any conductance-based model. Our work builds on these latter two papers and extends them to the situation where an M-current is present in the model.
To understand the implications of the codimension two bifurcations, we related them to the neural behaviours described above. The resonance property of neural models is quite simply related to the bifurcation that causes the loss of stability of the resting state when the applied current is increased. If this bifurcation is a Hopf bifurcation the model is a resonator, otherwise it is an integrator. As pointed out by Izhikevich, a Bogdanov-Takens bifurcation can switch the resonator type of a neuron [5]. Class I/II excitability was first linked to bifurcations in neural models by Rinzel and Ermentrout [23]. Rinzel and Ermentrout showed that neuronal models where the onset of repetitive firing occurs via a saddle-node bifurcation on an invariant circle are Class-I, while models where the onset occurs via a subcritical Andronov-Hopf bifurcations are Class-II. This link can be extended to other types of bifurcations by studying the associated F/I curves. The excitability class of individual neurons has been linked to the synchronisation properties of the neuron in a network through the phase resetting curve (PRC). In particular, it has been shown in certain circumstances that Class-I neurons have Type-I PRCs [24]. No conclusive link between Class-II neurons and a particular PRC type was found in that paper. More recently, Izhikevich has made a subtly different classification of excitability based on ramped current inputs as opposed to step current inputs. Izhikevich defines Class I/II excitability based on the bifurcation that causes the loss of stability of the resting state when the current is increased. Further, Izhikevich defines Class I/II spiking by the bifurcation that destroys the stable oscillations as the current is decreased [5]. A focus for this paper will be on how the presence of a BT point is linked to the emergence of a Hopf bifurcation and thus could be associated with a change of oscillation class for a conductance-based neural model.
The paper is organised as follows. In the next section, we provide a general conductancebased neuron model with the inclusion of the M-current and study the existence of the steady-state solutions. In Sects. 3 and 4, we give a complete characterisation of the BT bifurcation, provide a condition for the cusp bifurcation and discuss the existence of Bogdanov-Takens-cusp (BTC) bifurcation. In Sect. 5, we consider three example models and show that all three models exhibit the BT, CP and BTC bifurcation points. We construct bifurcation diagrams using MATCONT to explain possible behaviour of each example and use the numerical solution of each model to construct the frequency-current curves. In Sect. 6, we use numerical simulations to study the influence of varying of g M on the neurons synchronisation in two coupled neurons model with synaptic coupling. In Sect. 7, we discuss our results.

General model
In non-dimensional variables, a general conductance-based neuron model with the inclusion of the M-current can be written as follows: where a = (a 3 , . . . , a N ) T , where I app is the applied current and φ i is the set of indexes that represents the identities of the gating variables present in a given ionic current. In the rest of the manuscript, we assume that all conductances g j are positive, and the steady state activations w ∞ and a j,∞ , j = 3, . . . , N , are non-negative bounded functions (0 ≤ f (V ) ≤ 1), monotonic C 3 (R, R) and become sufficiently flat in the limits V → ±∞.

Equilibria
By applying the scaling t → t C m , system (1) can be written as where f 3 (V , a) = (f 33 (V , a 3 ), . . . , f 3N (V , a N )) T . Assume that (2) has an equilibrium point E * = (V * , w * , a * 0 ). From the equations above it follows that where V * satisfies Here, I ∞ is the steady-state I -V curve [5,23] defined by where I ion,∞ (V ) = I ion (V , a ∞ (V )) is the stationary ionic current. Notice that (3) can be written as

Bogdanov-Takens bifurcation
In the following we discuss Bogdanov-Takens point (BT point) of codimension two in (I app , g M )-plane, when all other parameters in the model are fixed. Assume that V * is a solution of (3), then there exist parameters (I * app , g * M ) such that It is well known [25][26][27] that the equilibrium point V * is a BT point if the zero eigenvalue has algebraic multiplicity two and geometric multiplicity one. Using an approach similar to [21,22], we obtain the following.
Theorem 3.1 Let V * be a solution of (3) at (I * app , g * M ) and define Then E * is an ordinary BT point of codimension two.
Then the Jacobian of (2) is Consequently, we have Notice that Thus, at E * , we have that the equation is equivalent to d dV I ∞ (V )| V * = 0. Thus, (0) = 0 when (6) holds.
It easy to check that Thus, at E * , (0) = 0 when (6) and (7) hold. Hence, λ = 0 is a double root. Now, we show that a Jordan block arises when λ = 0 is a double multiplicity root. In other words, when (6) and (7) hold, we demand the existence of four generalised eigenvectors q 0 , q 1 , p 0 , p 1 of A such that Then we obtain from Aq 0 = 0 the following equations: From (11) and (12), we have respectively. Hence, and it follows from (10) that Similarly, from Aq 0 = q 1 , A T p 1 = 0 and A T p 1 = p 0 , we have where I N-3 is the identity matrix of size N -3, and As the generalised eigenvectors must be non-zero, we let q 01 and p 11 be non-zero arbitrary constants. Thus, when (6) and (7) hold, equations (13)-(16) hold. Thus, four generalised eigenvectors exist. Hence, V * is an ordinary Bogdanov-Takens point.

Remark 3.1 With the additional condition
we can guarantee the uniqueness of the generalised eigenvectors q 0 , q 1 , p 0 , p 1 of A.
When V * is a BT point, system (2) has a two-dimensional centre manifold, with normal form given by (see, e.g. [25-27, 29, 30]) where where h 20 is the solution of the equation and the function G is defined as Here, ,j≤N is the Hessian matrix of a quadratic form at V * .
Proof From f m and f a , we have Hence, the components of G are Recall that all of these derivatives are calculated at V * . It follows from (8) that Thus, α 2 = 0 if and only if Consequently, from (19), we have Ah 20 = -G(q 0 , q 0 ), which has infinite solutions. This system is consistent due to the Fredholm solvability condition [26]. Hence, h 20 can be chosen such that β 2 = 0 in (18). This completes the proof.

Existence of bifurcations
Theorems 3.1 and 3.2 imply three bifurcations: BT, CP and BTC which are characterised by equations (5)- (7) and (20). In the following we discuss the solution of these equations. Recall that equation (5) relates the equilibrium point voltage value V * to I app and the other parameters.
Rearranging (6), we obtain where Similarly, (7) leads to where It is easy to check that the second derivative of I ∞ (V ) is Thus, (20) holds when Bogdanov-Takens bifurcation Suppose that there is V * that satisfies and I * app is given by (5).
Cusp bifurcation Suppose that there is V * that satisfies and I * app is given by (5).

Bogdanov-Takens-cusp bifurcation
Suppose that there is V * that satisfies Then there is an equilibrium E * = (V * , w * , a * ) that undergoes a BTC bifurcation at (I app , g M , and I * app is given by (5).
Remark 4. 1 We have explicitly included the leak current in our formulation. The leak current is not necessary for the occurrence of the BT and CP bifurcations. If g L = 0, then equation (21) becomes and the solution will go through as above. However, for the BTC bifurcation to occur, we must have another parameter to vary. We have shown that this third parameter can be the leak conductance g L . Solving the equations in a different way shows that the capacitance C M could also be used.

Implications
In the previous section we gave conditions which guarantee that a BT bifurcation can be induced by the variation of two parameters found in our general model: the applied current I app and the maximal conductance of the M-current g M . Further, if the conditions are met, we gave explicit expressions for the bifurcation point in terms of g M and I app . Near this bifurcation point the behaviour of the system will be described by the unfolding of the normal form (17) in terms of two parameters. The normal form and unfolding were first studied by [29,30]. The details can also be found in [25,27]. A key point for our work is that emanating out of the BT point are three codimension one bifurcation curves: Hopf bifurcation, saddle homoclinic bifurcation and saddle node (fold) of equilibria. A periodic orbit exists between the Hopf and homoclinic bifurcation curves, the stability of which depends on the sign of the coefficients α 2 , β 2 in (17). Thus the emergence of periodic solutions via a Hopf bifurcation can be linked to the presence of the BT point.
In the previous section, we also gave conditions which guarantee that a BTC bifurcation can be induced by I app , g M and the conductance of the leak current g L . The normal form and unfolding for the case considered in Theorem 3.2 was first studied in [31]; see also [21,26]. There are various possibilities for the bifurcations in the unfolding which are determined by the higher order terms in the normal form. The key results for our analysis are that in the three-dimensional parameter space there are two curves of cusp bifurcations and two curves of BT bifurcations with a surface of Hopf bifurcation starting at one BT curve and ending at the other. Near one BT bifurcation the Hopf bifurcation is supercritical (produces an asymptotically stable periodic orbit), while at the other it is subcritical. There is a saddle-node (fold) of limit cycles bifurcation associated with the change in criticality of the Hopf bifurcation. Fixing the value of one parameter (such as the leak conductance g L ) amounts to taking a two-dimensional slice in the three-dimensional parameter space. Thus, in general one should expect to see some subset of the bifurcations we have just described.

Numerical examples
In this section, we implement three examples with different ranges of g M corresponding to the range between the BT and cusp points. We apply our theoretical results and compare them with computations carried out with MATCONT [32]. We also construct bifurcation diagrams using MATCONT to explain the possible behaviour of each example. The labels used in these bifurcation diagrams are summarised in Table 1 supplemented by the dynamics for the gating variables h and n as in (1). Parameter values and other details of the model are given in the Appendix. Figure 1 shows the contour plot of equations (6), (7) and (20). In Fig. 1a, there are two intersections of equations (6) and (7)  The analysis of Sect. 4 shows that of the three curves, only the one defined by equation (6) depends on g L . Further, the representation (21) of this equation and the properties of the M-current show that increasing g L will move this curve downward. Given the shape of the curves in Fig. 1a, it is clear that increasing g L will move the BT and CP points closer together, and for sufficiently large g L we should obtain a single intersection point of all three curves corresponding to a BTC point. Figure 1b confirms that when we increase g L to 0.7507, we find the (approximate) BTC point (-46.6416, 7.75907, -0.0166046).
We use the MATLAB numerical continuation package MATCONT [32] to verify the theoretical results and supplement them with numerical bifurcation diagrams.  Table 3. This is consistent with our results in Fig. 1.  Table 3. (a) The conditions given by equations (6), (7) and (20) (26). Green curves are limit point (fold/saddle-node) bifurcations of equilibria, blue are Andronov-Hopf bifurcations, magenta are homoclinic bifurcations and red are limit point (fold) bifurcations of limit cycles. Codimension two bifurcation point labels are described in Table 1 Now, we discuss the switch in the model neuronal excitability class as g M increases. We plot a bifurcation diagram in the I app , g M parameter space for Wang-Buzśaki model (26) in Fig. 2. As expected from the normal form analysis, there is a curve of homoclinic bifurcations, a curve of Hopf bifurcation and a curve of saddle-node of equilibria emanating from the BT point. The Hopf is subcritical and thus an unstable periodic orbit exists for any parameters between the homoclinic and Hopf curves. See Fig. 2b. These curves are associated with the transition in the neuronal excitability class and show three cases.
• g M < g * M : In Fig. 3a, when g M = 0 and I app < 0.16, there exists a stable equilibrium point that determines the resting state and two unstable equilibria. As the applied current increases, the stable and one unstable fixed points collide in a saddle-node bifurcation point ("LP"). Consequently, a limit cycle is born simultaneously and emanates from the LP, that is, the limit cycle is created via a saddle-node on invariant circle bifurcation ("SNIC"). As expected, the oscillations on the limit cycle appear with arbitrarily slow frequency (see Fig. 4a), indicating Class-I excitability [5,23]. • g M > g M : For large enough g M , a different sequence of bifurcations is observed. In Fig. 3c, when g M = 3, at I app = 1, a limit point bifurcation of cycles "LPC" occurs giving rise to one unstable and one stable periodic orbit. Then, at I app = 1.1416, the unstable periodic orbit disappears in a subcritical Hopf bifurcation (subHopf ) of the lone equilibrium point, destabilising it. Consequently, firing with a positive frequency appears via LPC, and hence neuronal excitability Class-II occurs [5,23] Fig. 3b, an unstable limit cycle is created via a homoclinic bifurcation (the magenta curve in Fig. 2b) and disappears in the subHopf. In this case, the stable limit cycle appears via an LPC with a different unstable limit cycle which disappears via homoclinic orbit bifurcation (not shown in Fig. 2b). (ii) When 2.1 g M < g M , the sequence of bifurcations is very similar to that for g M > g M . An LPC bifurcation creates both unstable and stable periodic orbits. The former is lost in the subHopf, see Fig. 5c. For all g M ∈ (g * M , g M ), there is a region of bistability between a stable limit cycle and a stable equilibrium point, between the LPC and subHopf bifurcations. Consequently, when g * M < g M < g M , a neuronal excitability Class-II occurs [5,23].  Table 1. Bottom row: corresponding F/I curves Therefore, the model neuronal excitability type switches from Class-I to Class-II when the conductance of the M-current g M passes through the BT point. Now let us consider the effect of the leak conductance g L . As shown above, increasing g L monotonically decreases the g M value at the bio-physically permissible BT point. This means that the range of values of g M where the model has Class-I excitability will be decreased. Equivalently, smaller changes of g M are needed to switch the model from Class I to Class II. If g L is increased enough, then g * M may become negative, in which case the model will exhibit Class-II excitability regardless of the value of g M .
Parameter values and other details can be found in the Appendix. Solving equations (6), (7) and (20)  Applying the analysis of Sect. 4 to this model also shows that increasing g L should lead to a BTC point. This is confirmed in Fig. 6b. We find the BTC point  Table 4. (a) The conditions given by equations (6), (7) and (20) (27). Green curves are limit point (fold/saddle-node) bifurcations of equilibria, blue are Andronov-Hopf bifurcations, magenta are homoclinic bifurcations and red are limit point (fold) bifurcations of limit cycles (LPC). Codimension two bifurcation point labels are described in Table 1 (-43.1385, 2.9461, 0.0008) when we increase g L to 0.3785. As in the previous example, the neuronal excitability type switches from Class I to II as the conductances of the M-current increase, Class I when g M < g * M and Class II otherwise, see Figs. 7, 8 and 9. Although the range (g * M , g M ) is much smaller than for Example 1, model (27) exhibits a similar behaviour in this range, see Fig. 5.
Example 3 The reduced Traub-Miles (RTM) model is a substantial simplification of a model of a pyramidal excitatory cell in rat hippocampus due to Traub and Miles [34]. The RTM model with the M-current can be written as [35] C m dV dt  Fig. 10a. Applying the analysis of Sect. 4 again shows that increasing g L should lead to a BTC point. This is confirmed in Fig. 10b. When we increase g L to 13.79, the BT and CP points collide producing the BTC point (-49.8762, 166.25, -0.6745). In this example, we notice that the range (g * M , g M ) is bigger than those in Example 1 and 2, but the transition in the neuronal excitability type is consistent with previous examples: Class I when g M < g * M and Class II otherwise, see Figs. 11 and 12.

Implications for synchronisation
In Sect. 4 we showed that the M-current will give rise to a BT bifurcation in any conductance-based neural model, when certain conditions are met. In Sect. 5 we showed in three examples that these conditions are met and a BT bifurcation occurs. Further, we showed that this BT bifurcation induces a transition from Class-II to Class-I excitability in these models as the conductance of the M-current is decreased (as would be the case in the presence of acetylcholine). In this section we explore one implication of this transition. There are many studies in the literature describing the relationship between the synchronisation of coupled neurons and their neuronal excitability type, see, e.g. [8,24]. The classic result is that the in-phase solution of a pair of weakly coupled Class-I oscillators Figure 10 Existence of codimension two and three bifurcation points in Wang-Buzśaki model (28) with the parameter values given in Table 5. (a) The conditions given by equations (6), (7) and (20)  Bifurcation diagram in the I app , g M parameter space for RTM model (28). Green curves are limit point (fold/saddle-node) bifurcations of equilibria, blue are Andronov-Hopf bifurcations, magenta are homoclinic bifurcations and red are limit point (fold) bifurcations of limit cycles (LPC). Codimension two bifurcation point labels are described in Table 1 Figure 12 One-parameter bifurcation diagrams for RTM model (28), showing the change in the bifurcation structure as g M is varied. (a) g M < g * M (the value at the BT point); (b) g * M < g M < g M ; (c) g M > g M (the value at the CP point). Green/blue curves show stable/unstable equilibria. Pink curves show maxima/minima of periodic orbits. Codimension one bifurcation point labels are described in Table 1 model with synaptic coupling is stable when there are inhibitory coupling and unstable for excitatory coupling, while the anti-phase solution exhibits the opposite stability [24]. The synchronisation of Class-II oscillators is less clear, and other factors such as the synaptic time constants and firing frequency may affect these conclusions [24,33]. By in-phase solution, we mean both oscillators reach their highest peak at the same time, whereas an anti-phase solution means one oscillator reaches its highest peak one half-period after the other oscillator.
To study the stability of phase-locked solutions and the correspondence with the neuronal excitability type as g M varies, we write two coupled neurons with synaptic coupling as follows: where I ion are ionic currents in Examples 1-3. The synaptic coupling function and parameters are given in Table 2.
To determine the stable phase-locked solution(s), first we solve (29) numerically with ten random initial conditions at each step of g M , then we calculate the period of the oscillators (T 1 and T 2 ) in the numerical solution. Finally, we approximate the phase shift as where · is the floor function, T = (T 1 + T 2 )/2 and τ is the argument shift satisfying V 1 (t) = V 2 (t + τ ) for all t. Figure 13 shows bifurcation diagrams for (29) with excitatory and inhibitory synaptic coupling in Examples 1-3. For instance, for coupled Wang-Buzsaki model (Example 1), we notice in Fig. 13a that when g M < g * M (Class-I dynamics in (26)), the in-phase solution is unstable and the anti-phase solution is stable with excitatory coupling V syn = 0. The reverse is true for inhibitory coupling V syn = -75. This is consistent with [24]. When there is an excitatory synaptic connection, as the M-current reaches g M ≈ 0.5, the anti-phase solution loses its stability and two stable out-of-phase solutions (neither in-phase nor anti-phase) appear. As the conductance of the M-current is increased any further, a stable in-phase solution appears. Hence, there is a transition from stable antiphase solution to stable in-phase solution via stable out-of-phase solutions. The transition also occurs at g M ≈ 0.5 when there is the coupling is inhibitory. We observe a similar dynamical behaviour in Examples 2 and 3, see Fig. 13b-13c, although the transition is not as Table 2 Synaptic coupling function and parameters in (29)  phase-locking in anti-phase (phase difference π ), while inhibitory coupling leads to in-phase (phase difference 0). In all cases g M has to be increased significantly past the BT value before the solution switches to in-phase for excitatory coupling and anti-phase for inhibitory coupling clear in all cases. Although the relationship of the transition point to the codimension two bifurcations varies with the different models, in all cases it occurs at some g M ∈ (g * M , g M ), that is, when the model has Class-II dynamics.
As indicated above, other factors may affect the synchronisation of neurons. We focus here on the firing frequency of the neuron. In [33] it was shown that increasing the firing frequency by increasing the applied current could switch the PRC of a model neuron with an M-current from Type-II to closer to Type-I. In [13], the authors reproduced this result for other neural models and studied how changing the firing frequency modulates the synchronisation properties induced by the M-current. They found that synchrony in excitatory networks of neurons with a Type I PRC (low g M ) was largely unaffected by frequency modulation, whereas networks of Type II PRC neurons (high g M ) synchronised much better at lower frequencies. In [7], the authors studied how the stability of in-phase and anti-phase phase-locked solutions in Wang-Buzśaki model (with no M-current) varied with firing frequency. At low frequencies with inhibitory coupling, they showed that both in-phase and anti-phase phase-locked solutions were stable. However, at higher frequencies only the in-phase solution was stable. In contrast, with excitatory coupling, they showed that the in-phase solution was unstable for both high and low frequencies. Recalling that the Wang-Buzśaki model is a Class-I oscillator, this latter result is consistent with that of [13].
To consider if firing frequency has an effect in our results, we determined the variation of firing frequency with the conductance of the M-current g M for our example models, see Fig. 14. In all cases the firing frequency decreases rapidly as g M increases. When the models are in the Class-I excitability regime (below the BT point), the frequency change does not affect the synchronisation properties. This is consistent with the results described above [13,33], given that neurons with Class-I excitability typically have Type-I PRCs [24]. Recalling that the main switch in synchronisation behaviour in all cases occurs within the Class-II regime, we conclude that this switch is likely due to the decrease in the frequency as g M increases.  Fig. 13 In summary, while the excitability class of the model changes exactly at the BT point, the synchronisation property of the models switches at g M value larger than the BT point, when the frequency of the intrinsic oscillations is small enough.

Discussion
In this paper, we studied Bogdanov-Takens (BT) bifurcation in a general conductancebased neuron model with the inclusion of the M-current. We started by showing the existence of equilibrium points. Then we derived the necessary and sufficient conditions for the equilibrium point to become a BT point. A degenerate Bogdanov-Takens (BTC) point appears when BT and cusp points merge. To discuss the occurrence of such a point, we provided the condition for a cusp bifurcation. We then showed that the conditions for the BT and cusp bifurcation may be satisfied by varying the applied current and the maximal conductance of the M-current and that for the BTC point by additionally varying the conductance of the leak current.
As previously noted, our theoretical work was inspired by two recent papers. In [21] they show that the BTC point can occur in any conductance-based model in the parameter space of the applied current, leak conductance and capacitance. They use this to study the effect of the leak current on the excitability properties of models for single neurons and synchronisation properties for networks of neurons. In [22] they study a general conductance-based neural model. They show that if the model has an equilibrium point with a double zero eigenvalue for some parameter values, then it is a BT point. Further, they give conditions on the gating variables and time constants for a BT bifurcation to occur. They propose the BT normal form as a generic minimal model for a single neuron.
Numerically, we applied our analytical results to three examples and compared them with the computations of MATCONT, a numerical bifurcation analysis toolbox in Matlab. Furthermore, we constructed bifurcation diagrams using MATCONT to explain the possible behaviour of each example and discuss the switches in the neuronal excitability class with respect to the M-current g M . As predicted by normal form theory [25,27,29,30] in all examples a curve of homoclinic bifurcation, a curve of Hopf bifurcation and a curve of saddle-node of equilibria emanate from the BT point. These latter two curves particularly affect the neuronal excitability class. We found that a transition is determined by the BT point which occurs at (g M , I app ) = (g * M , I * app ). The model is a Class-I oscillator when g M < g * M and Class-II when g M > g * M . More precisely, when g M < g * M as I app is increased, oscillations with arbitrarily slow frequency appear via a saddle-node on invariant circle bifurcation, while when g M > g * M oscillations with a positive frequency appear via a fold bifurcation of cycles, followed by a subcritical Hopf bifurcation.
Using systems of two synaptically coupled cells, we explored how the change in excitability class with the variation of g M affects synchronisation in the example models. We found that while the excitability class of the model changes exactly at the BT point, the synchronisation property of all the models switches at g M value larger than the BT point. We attributed this change to the fact that the M-current also affects the frequency of the intrinsic oscillations and that the synchronisation of Class-II oscillators has been shown to be sensitive to intrinsic frequency. Thus the necessary condition for the switch of synchronisation, we observed, is that the system be Class-II and the frequency be sufficiently small.
We also considered the effect of the leak conductance g L showing that, in the examples we considered, increasing g L decreases the g M value of the BT point. This means that the range of values of g M where the model has Class-I excitability will be decreased. Equivalently, smaller changes of g M are needed to switch the model from Class I to Class II. If g L is increased enough, then g * M may become negative, in which case the model will exhibit Class II excitability regardless of the value of g M . Since the switch of synchronisation occurs at a higher value of g M than the BT point, this does not necessarily mean that the system will not exhibit changes in synchronisation associated with a change in g M , it just means that smaller changes in g M are needed to switch the synchronisation property. We note that Prescott et al. [11,12] represented the increase in membrane conductance due to background synaptic input using a leak current with a reversal potential near rest in a Morris-Lecar model with an M-current. The one-parameter bifurcation diagrams in [12] are consistent with what we have seen in our analysis.
Our analysis of the effect of g L on the BT point relies on understanding how the intersection points of two curves vary with g L . Only one curve depends on g L , and we can show in general (i.e. for any model) that the curve will move downward as g L increases. This effect depends on two aspects of the M-current: the reversal potential is a large negative value (since it is a potassium current) and the current is non-inactivating, see equation (24).
The implications of these results for the action of acetylcholine are as follows. If the neuron is of Class-II in the absence of acetylcholine (corresponding to high g M ), then the presence of acetylcholine may push the system past the BT bifurcation point and change the neural excitability type to Class-I. The expected synchronisation in the presence of sufficient acetylcholine is then clear: neurons with excitable coupling will likely desynchronise, while those with inhibitory coupling will synchronise. This is consistent with the changes to the PRCs induced by acetylcholine observed in [1]. Whether or not acetylcholine induces a change in synchronisation may depend on intrinsic firing frequencies of cells. Expanding on the idea of Prescott et al. [11,12], an increase in membrane input conductance would make the system more sensitive to the effects of acetylcholine, so that switches of synchronisation could occur more easily.
These conclusions, of course, assume that the only effect of acetylcholine is to downregulate the M-current. However, acetycholine has been observed to have other effects, including down-regulating an afterhyperpolarization current I AHP [3,37] and the leak current [1]. As indicated above, our work indicates that decreasing g L will increase the value of g * M . Thus the simultaneous down-regulation of the leak and M-currents would cause the switch of excitability class at higher values of g M . The net effect would be to increase     (28) Conductance (mS/cm 2 ) Reversal potential (mV) Capacitance (μF/cm) g L = 0.1 V L = -67 C M = 1 g Na = 100 V Na = 50 g K = 80 V K = -100 • Example 3: Reduced Traub-Miles model. The infinity and τ σ , σ ∈ {m, h, n, w}, functions are: