The Bogdanov–takens Normal Form: a Minimal Model for Single Neuron Dynamics

Conductance-based (CB) models are a class of high dimensional dynamical systems derived from biophysical principles to describe in detail the electrical dynamics of single neurons. Despite the high dimensionality of these models, the dynamics observed for realistic parameter values is generically planar and can be minimally described by two equations. In this work, we derive the conditions to have a Bogdanov–Takens (BT) bifurcation in CB models, and we argue that it is plausible that these conditions are verified for experimentally-sensible values of the parameters. We show numerically that the cubic BT normal form, a two-variable dynamical system, exhibits all of the diversity of bifurcations generically observed in single neuron models. We show that the Morris–Lecar model is approximately equivalent to the cubic Bogdanov–Takens normal form for realistic values of parameters. Furthermore, we explicitly calculate the quadratic coefficient of the BT normal form for a generic CB model, obtaining that by constraining the theoretical I-V curve's curvature to match experimental observations, the normal form appears to be naturally cubic. We propose the cubic BT normal form as a robust minimal model for single neuron dynamics that can be derived from biophysically-realistic CB models.


Introduction
Our purpose here is to present a robust minimal model for single neuron dynamics.We have a dynamical systems approach, and the detailed calculations can be found in [1].
Near a bifurcation, the behavior of any physical system is universal (i.e., not dependent on its specific details) and is described by a universal equation called the normal form [2,3].These equations are universality classes, and very different dynamical systems (near a given bifurcation) will be described by the same normal form.Although this is mathematically strictly true infinitesimally near the bifurcation, it occurs often that the qualitative aspects of the behavior of the system are still given by the normal form, even outside of the infinitesimal neighborhood of the bifurcation point.In this work, we will show that the electrical dynamics of a single neuron can be an example of this fact.
At the present time, to describe neuronal dynamics for different types of neurons, there is an overwhelming diversity of thousands of high dimensional nonlinear models.This class of models is called conductance-based models (CB models), and they are considered the most biophysically realistic models [4][5][6][7].In the last few decades, the detailed electrophysiological and biophysical description of the wide diversity of conductances [8] led scientists to the possibility of building a detailed CB model of many of the neuronal types studied [5].Each different CB model, built using adjustments of parameters with the experimental data, can have more than twenty nonlinear equations.This has generated in the last few decades a massive production of different CB models for different types of neurons or experimental scenarios and the development of databases of thousands of different CB models [9].
Although there is a huge diversity of types of neurons (and CB models), there is a robust dynamics observed experimentally (and numerically) in almost all of them.This distinctive electrical dynamics observed in neurons is called the excitable behavior of neurons.In 1948, it was Hodgkin who classified this excitable behavior observed in neurons into two major groups [10]: (1) Class I neural excitability: action potentials can be generated with arbitrarily low frequency, depending on the strength of the applied current; (2) Class II neural excitability: action potentials are generated in a certain frequency band that is relatively insensitive to changes in the strength of the applied current.
In the late 1990s, a connection between this classification and bifurcation theory was noticed [11].It was realized that Class I neural excitability is consistent with a saddle-node homoclinic bifurcation.On the other hand, Class II neural excitability was consistent with a local subcritical Andronov-Hopf bifurcation with a narrow range of bistability.However, it has been shown that this classification is consistent with a broader set of global bifurcation, such as saddle-node homoclinic bifurcation, the saddle-homoclinic bifurcation and the big homoclinic bifurcation, which can lead to secondary local bifurcations.Therefore, despite the huge diversity of types of neurons and models, they display a universal dynamics generically captured by a two-variable dynamical system [7,12,13].These minimal models can be phenomenological models or reduced CB models.Procedures have been developed to reduce the number of variables in CB models [7,14].However, the mathematical mechanisms by which single neurons display this universal planar dynamics are still not fully clear.
We propose here a solution to this problem adopting a dynamical systems approach, with the idea that we have to look for a robust behavior related to the experimental data and construct from its analysis the simplest possible model satisfying the requirements.
We first show that due to the very special mathematical structure of CB models if we impose the conditions to have in some neighborhood of the space of parameters both the saddle-node and the Andronov-Hopf bifurcations, we naturally arrive at the Bogdanov-Takens (BT) bifurcation point.Furthermore, we find that these conditions are equivalent to having both positive and negative feedback and a time scale separation, which are precisely two conditions that have been reported as a key dynamical mechanism for neurons to exhibit excitable behavior [15,16].However, we still have to impose that the system should be in the neighborhood of the transition between having one or three fixed points, since both situations are generically observed experimentally [7].We show that this last requirement is satisfied imposing a third condition, which implies that the coefficient of a quadratic term in the BT normal form at the critical point vanishes.We have then to include cubic nonlinearities at the critical point, and the quadratic term will appear with a small coefficient as an unfolding term (as happens in a subcritical bifurcation).The calculation of the cubic BT normal form gives us explicit formulas for the coefficients in terms of the parameters of the CB model under consideration.Additionally, all of the original variables of the CB model are expressed in terms of the critical variables, which we can choose to be (u, u), where u is the membrane potential and u its time derivative.The coefficients and the original physical variables are given by analytical expressions, which result from a standard calculation using the results in [1,17,18].We present here the BT cubic normal form in Arnold's choice, which leads to a second order in the time differential equation for the variable u, which has a very appealing mechanical interpretation as a Hamiltonian system with a nonlinear friction.We investigate numerically this equation, and we find that it robustly exhibits the local and global bifurcations observed in single neurons.

Mathematical Structure of Conductance-Based Models
CB models are a class of ODEs with a special structure that can be completely characterized.In any CB model, there is one variable that in this paper we will call u that represents the membrane potential and other set of variables that we will call m j that represents the gating probability of a ion channel sub-unit over time.In nondimensional variables, a general form of a CB model can be written as follows: The functions m ∞ j (u; σ j ) represent the stationary probability of a gating variable m j to be activated, and τ j (u; σ j ) is its decay time.Typically, m ∞ j (u; σ j ) and τ j (u; σ j ) are functions of the membrane potential u and depend on the same set of parameters σ j .These parameters are shaped by the biophysical structure of the ion channel and its kinetic properties.We will call the vector formed by the union of all these parameters σ T (i.e., σ T = ∪ k j=1 σ j ).The symbol I represents the input current to the neuron, and the function I ionic (u, m; η) represents the overall ionic current that goes through the membrane.This last quantity depends on parameters, such as the reversal potential and the maximal conductance of a specific ionic current.We will use the vector η to collect together all of these parameters of different ion channels (i.e., η = (u 1 , . . . ,u n , g 1 . . . ,g n )).In general, for any conductance-based model, the overall ionic current can be written as follows: Each I ionic j represents an ionic current associated with an ionic conductance j (or ion channel j) and, in general, can be written as: where g j is the maximal conductance and u j is the reversal potential associated with the j-th ionic current.For a given current I ionic j , there is a specific set of gating variables ( m l : l ∈ w j ) that dynamically controls the current (see Equation ( 3)).We will define w j as the set of indexes that represents the identities of the gating variables present in a given ionic current j (e.g., for the Hodgkin and Huxley model, the sodium current has the gating variables m and h).Because each gating variable m l represents the probability of a specific sub-unit of a given ionic channel j to be open, the same gating variable cannot belong to two different ionic currents.This restriction is based on the nature of the ionic channel's biophysics: the activity of one ionic channel does not directly affect the activity of another ionic channel (indirectly, it does through the coupling with the potential u), since their mechanisms of opening and closing are not physically connected [19].Then, all of the sets of indexes w j are disjoint.
Using Definition (2), the total ionic current reads: If we consider long times, such that all of the gating variables have relaxed ( x = 0), the ionic current written as (5) becomes the well-known stationary I-V curve widely used in electrophysiological studies [16,19].We will call the stationary ionic current I ∞ (u; σ T , η) the I-V curve of our model.It is important to notice that this function depends only on the membrane potential (u) and the sets of parameters σ T and η, but not the gating variables ( x).The stationary ionic current is given by: We will define the transient current as I T (u, x; σ T ), as the part of the total current that depends on both u and x and disappears for long times, such that all of the gating variables have relaxed, i.e.: Note that I T (u, x; σ T , η) is at least linearly proportional to any of the variables x j .This means that when all of the gating variables decay to the stationary value (i.e., m = m ∞ (u; σ j )), then the transient current is zero I T (u, x; σ T , η) = 0, which is the definition of a variable that acts transiently in this context.Now, we will define: Note that due to the form of m ∞ j (u; σ j ), which is a strictly increasing (or decreasing) function of u, one has that β j (u; σ j ) is strictly negative (or positive).Then, any CB model can be re-written in the following form: The dynamics of each gating variable m j is a linear decay to its stationary value (which corresponds to x j = 0) plus a perturbation driven by the difference between the input current and the ionic current and weighted by the derivative of their corresponding stationary gating variable function.The conditions to have a fixed point are: A fixed point occurs when all of the gating variables have relaxed to their stationary value (i.e., m j = m ∞ j (u; σ k )) and when there is no net flux of current through the neuron membrane.
The fixed points must satisfy x = 0 and I = I ∞ (u; σ T , η).This equation has at least one root, since the function I ∞ (u; σ T , η) for u → ±∞ can only have a non-negative and constant slope (because the conductances g j are always positive quantities), but generically has one or three [4,6,7].Then, the existence or disappearance of fixed points in this N + 1 dimensional dynamical system depends on a one-dimensional equation I = I ∞ (u; σ T , η).We can choose a root u * of this equation as a continuous parameter of bifurcation instead of the parameter I, which is the one that is experimentally tuned, as long as we stay in the same branch of the equation.Thus, by doing the standard translations in the neighborhood of the fixed point ( x = 0 + δ x and u = u * + δu) and after some algebra using Equation ( 9), we find that the linear matrix has the form: where: and: and: The structure of the matrix in (11) is known as a rank-1 updated matrix.As the name states, the matrix L is the sum of the diagonal matrix A plus a rank-1 perturbation, which corresponds to β M T .We will use this fact to provide a closed form to the characteristic polynomial of the linear approximation for CB models.For simplicity, throughout the rest of the text, we will not write explicitly the dependence on the parameters σ T and η.

Characteristic Polynomial
Using known expressions [20,21] for the determinant of a rank-1 updated matrix, we have that the characteristic polynomial of any conductance-based model can be written as: Providing that λ = −1/τ j ∀j, which is equivalent to: The last expression is the general form for the characteristic polynomial of any conductance-based model.Then, this expression can be used to find analytic formulas for local bifurcations.To our knowledge, this is the first time that an analytic expression is given for the characteristic polynomial of a conductance-based model with an arbitrary number of conductances.For simplicity, in this section, we did not include explicitly the dependence on the fixed point u * .

Bogdanov-Takens Bifurcation in Conductance-Based Models
We are looking for a critical point that encloses most of the diverse, but planar (two variables) dynamics observed in single neurons.In electrophysiological experiments, when the transmembrane current I is increased, neuron dynamics abruptly transit from a fixed point (rest potential) to a limit cycle (spiking).This transition to relaxation oscillations is the landmark of single neuron dynamics.Often, different classes of neurons have different types of spiking dynamics.The rest potential can exhibit an exponential decay after being perturbed (node) or damped oscillations (focus).In addition, the transition to spiking can occur via many dynamical mechanisms as saddle-node homoclinic bifurcation, bistability between a focus and limit cycle, subcritical Andronov-Hopf bifurcation between others.These dynamical mechanisms are seldom controlled by the experimenter, but they depend on the specific class of neuron, that is its specific set of ionic conductances and the values of its parameters.When the current I is increased, depending on the type of neuron, the local bifurcations observed are usually saddle-node or Andronov-Hopf.Then, the first condition that we will impose for a critical point that encloses the diversity observed in single neurons dynamics is to have both local saddle-node and Andronov-Hopf bifurcation in its neighborhood.This must be a bifurcation point of at least codimension two, since we need one parameter (usually I) to reach the bifurcation and another to transit from saddle-node to Andronov-Hopf (or vice versa).We exclude the bifurcation arising from the eigenvalues (0, +iω, −iω), i.e., a simultaneous saddle-node and Andronov-Hopf bifurcation, since it leads to a three-variable dynamics.A low codimension candidate for this is a double zero bifurcation (two eigenvalues equal to zero).Slight departures in the parameter space from this point will transform this double zero in a saddle-node or an Andronov-Hopf bifurcation locally.Using the characteristic polynomial (Equation ( 15)), we can find analytically the conditions to have double zero eigenvalues: Since the relaxation times of the gating variables (τ i (u * )) are always positive quantities, then these two conditions are equivalent to: Using Expressions ( 6) and ( 7), we obtain that the previous expression is equivalent to: At this point, it is necessary to define the following quantity: The function τ 0 (u * ) can be interpreted as the relaxation time of the ohmic part of the I-V curve.The function G j (u * ), called the gating function, will quantify the linear contribution of a given gating variable j to the overall CB model dynamics.The sign of G j (u) will determine if a gating variable is an amplifier (i.e., contributes to the dynamics as a positive feedback, is negative) or resonant (i.e., contribute to the dynamics as negative feedback, is positive) according to the classification used in computational neuroscience [7].This naturally leads us to a sign criteria to designate an amplifier or a resonant gating variable.Using these definitions, we can write the double zero condition in terms of the gating functions and relaxation times: Since the relaxation times (τ 0 (u * ) and τ j (u * )) are positive quantities, then after some simple algebra, we can conclude that the necessary conditions for having a double zero bifurcation are: (1) at least one amplifier and one resonant gating variable (i.e., at least one G j (u * ) must have a different sign); (2) at least one gating variables relaxation time is different.In other words, a conductance-based model must have gating variables that produce both positive and negative feedback, and there must be a time scale separation between at least one of them.Furthermore, we showed that if the system has two time scales, the fastest gating variables must be an amplifier and the slower resonant (Supplementary Materials, Section 2).This result is surprising, since the conditions mentioned have been reported as a key dynamical mechanism for neurons to exhibit excitable behavior [15,16].In addition, we showed that in CB models, a double zero bifurcation is always a BT bifurcation (Supplementary Materials, Section 1.1).Additionally, we calculated the complete Jordan base (Supplementary Materials, Section 1.2) and showed that generically, a Jordan block will arise.An example of the BT bifurcation in the widely-known Hodgkin and Huxley CB model [15] is shown in Figure 1B.

Additional Condition to the Bogdanov-Takens Bifurcation
We already identified BT as a critical point of interest.However, there is an important feature that is commonly observed in neurons and neuron models that was not included in the later discussion.It has been observed in experiments and simulations that under realistic parameter values, the I-V curve is usually either monotonically increasing or has two changes in its concavity (see Figure 1A).As we discussed in Section 2, the zeros of the I-V curve determine the fixed points of the dynamics.Then, if the I-V curve is monotonically increasing, the dynamics will have just one fixed point, and if its concavity changes just twice, it has three.Nevertheless, the stability of the fixed point depends on the eigenvalues of the linear matrix, which can be calculated using the expression for the characteristic polynomial in Equation ( 15), and they do not depend explicitly on the second derivative of the I-V curve.Since we want to include both scenarios (one fixed point and three fixed points) in the neighborhood of the critical point, then we need to impose that: Entropy 2015, 17, 7859-7874 In addition to the two BT conditions, this third condition puts the critical point in the transition point between one to three fixed points.Then, in its neighborhood, one finds many of the local bifurcation scenarios previously observed in neuron models: Andronov-Hopf bifurcation with one fixed point, Andronov-Hopf bifurcation with three fixed points and a saddle-node with three fixed points.When the BT normal form is calculated in this point, we obtain surprising results regarding its universality to reproduce single neuron dynamics.
) function is shown in green.The plots in black frames show the eigenvalues of the linear matrix L(u * ) in the complex plane.We omitted the greater in absolute value eigenvalue (far to the left) in order to better visualize the two critical ones.Tags 1-3 indicate the corresponding values of the variable u * .The plot in a red dashed frame shows the eigenvalues of the linear matrix L(u * ) in the complex plane when the model undergoes a BT bifurcation.

The Cubic Bogdanov-Takens Normal Form
Using the results in Elphick, Tirapegui, Brachet, Coullet and Iooss in [2], it is possible to characterize and calculate explicitly all of the non-linear terms for the BT normal form.In the case of the BT bifurcation, the equations obtained can be satisfied in more than one way, leaving an inherent freedom to incorporate the normal from (see Section 3, Supplementary Materials).We have two extreme choices: Arnold's or the Takens choice [22].Since in Arnold's choice, we gain a Hamiltonian intuition, we make this choice here and use the Arnold form of the cubic BT normal form, which can be written as: We will call the function F(u, α) the force and the function λ(u, β) the friction using the obvious analogy to mechanics.More explicitly, F (u, α) = α + γ 2 u 2 + γ 3 u 3 , and λ (u, The unfolding coefficients are α, β and γ 2 .When they vanish, we are at the critical point, and when they take small values, we are in the neighborhood of the critical point, where we expect the normal form to be valid quantitatively.Outside this neighborhood, one expects a qualitative description.We have three unfolding parameters, since the codimension is three.We can define a potential V(u, α), such that the force F (u, α) is given by F (u, α) = − ∂V(u,α)  ∂u , where Then, we can define the "energy" of the BT normal form as E(u, α) ≡ 1 2 u2 + V (u, α) and conclude that d dt E(u, α) = − u2 λ(u, β).This last equation shows that the BT normal form does not conserve the energy.The reader should notice that the BT normal form can undergo an injection or dissipation of energy depending on the sign of λ(u, β), and this sign can of course change when u changes.
The most common scaling studied in the literature is quadratic, i.e., F(u) = α + γ 2 u 2 and λ(u, β) = β + λ 1 u.The complete bifurcation diagram is known, and it exhibits saddle-node, Andronov-Hopf and homoclinic bifurcations [22,23].However, the dynamics observed in this normal form is not consistent with the dynamics observed in single neurons.First, it is not bounded in the complete phase space, which means that for some regions in the phase space, u can go to infinity (see Figure 2 left plot).Second, it does not display many of the distinctive global bifurcations observed in neurons as the saddle-node homoclinic bifurcation or the saddle homoclinic bifurcation.

Injection of energy zone
Dissipation of energy zone Since we know analytically the Jordan basis, we are able to calculate explicitly the cubic BT normal form using the methods developed in [2]; that is, to calculate for any model of the form of Equation ( 1), the explicit dependence of the unfolding parameters (α, β, γ 2 ) and the coefficients (γ 3 , λ 1 , λ 2 ) on the parameters of the model.We find that the critical variables of the normal form can be chosen to be (u, u), where u is the dimensionless membrane potential and u its time derivative (see Equation ( 1), Section 1.1, Supplementary Materials).The calculations, which are standard after the results in [2], are presented for the quadratic coefficient γ 2 (see Section 3, Supplementary Materials).Very surprisingly we found that the quadratic term of the force is proportional to the second derivative of the I-V curve, i.e.: This relation is very important, since it shows that in order to have one or three fixed points (which is what is observed), then in the unfolding of the normal form, we have to impose an additional condition, namely that the second derivative of the I-V function with respect to u vanishes.
As is shown in Figure 2 (right plot), for some values of u, i.e., , there will be only the injection of energy, and it can be shown analytically that the system cannot have a stable fixed point in this region, no matter the form of the force.Therefore, all of the orbits will escape from this region, but the orbits will not go to infinity, because as soon as they reach the regions , the system will dissipate energy quadratically with the variable u and the dynamics will be bounded to a certain region of phase space.Then this scaling will present a bounded dynamics for sensible force functions.Moreover, it has been shown that the BT normal form with a quadratic friction and cubic force presents a rich set of global bifurcations, some of which are characteristic of single neuron dynamics [24].The cubic BT normal form has the following form: This equation corresponds to the normal form of a particular BT bifurcation point of codimension three described by Equations ( 16)- (18).We expect that this normal form displays at least some of the fundamental aspects of single neuron dynamics, since we construct this critical point imposing constrains widely observed among neurons.However, it is not clear if this equation is able to reproduce the diversity in the observed dynamics of neurons.We show that it does in the next section.

Numerical Results and Comparison with the Morris-Lecar CB Model
Simulations of the cubic BT normal form show that the most typical global bifurcations mentioned in the literature for neuron models (i.e., saddle-node homoclinic bifurcation, the saddle-homoclinic bifurcation and the big homoclinic bifurcation) [4,7] are generically displayed by this equation.

Before the bifurcation
In the bifurcation After the bifurcation The top plot in Figure 3 shows the saddle-node homoclinic scenario.Before the bifurcation, we have three fixed points: a stable node, a saddle and an unstable focus.As is shown numerically, if we follow the orbits in green, these two fixed points have a heteroclinic connection.The unstable focus is inside the heteroclinic connection.When we move the parameter α, these fixed points collide, and the heteroclinic connection becomes a limit cycle; this is the point where the saddle-node homoclinic bifurcation occurs.The unstable focus remains inside the limit cycle.This scenario is a distinctive bifurcation of Class I neurons.
The middle plot in Figure 3 shows the big homoclinic scenario.As in the last scenarios, there are three fixed points: a stable node, a saddle and an unstable focus.When we decrease the parameter β, the unstable manifold of the saddle fixed point returns, getting closer to the stable manifold.When the unstable manifold and the stable manifold connect by the same orbit, then the homoclinic orbit formed traps the other ends of the unstable and stable manifolds of the saddle, then the big homoclinic bifurcation occurs.This homoclinic connection becomes a limit cycle with the three fixed points inside.The limit cycle is separated of the unstable focus by a heteroclinic separatrix, as is shown in the third graphic of the figure in black.The separatrix connects the stable manifold of the saddle point with the unstable focus, forming a closed figure.This connection is generic (codimension zero) because it connects a stable manifold (dimension one) with the unstable manifold (dimension two) of the focus.This separatrix may suffer a saddle-homoclinic bifurcation, and we finally obtain a big limit cycle, a stable focus and a limit cycle separatrix in between.Then, moving the parameters, this limit cycle could collide and disappear.Together, these scenarios are distinctive for Class II neurons.
The bottom plot in Figure 3 shows the saddle homoclinic scenario.There are three fixed points: a stable node, a saddle and an unstable focus.When the parameter β increases, the unstable manifold of the saddle fixed point returns, getting closer to the stable manifold.When the unstable manifold and the stable manifold connect by same orbit, the saddle homoclinic bifurcation occurs, and then, the homoclinic connection becomes a limit cycle.The unstable focus stays inside the limit cycle.The limit cycle (marked with 1 in the figure) and the stable node (marked with 2 in the figure) are separated by the stable manifold of the saddle point.This bistable behavior has been reported before in neural models [4,7].In general, this bistable scenario is not necessarily the outcome of the saddle homoclinic bifurcation, but in neuron models, this is often the case.This scenario is also distinctive for Class II neurons.
We have compared these results with the dynamics observed in the Morris-Lecar model [13].This CB model has been widely studied in computational neuroscience, and it is known that it displays much of the diversity observed in single neuron dynamics [4].This is a minimal CB model, in the sense that it has the minimal number of variables for presenting oscillations (two variables).The nondimensional form of the Morris-Lecar model reads: where m ∞ (u), n ∞ (u) and τ(u) are hyperbolic functions (see Equations ( 34)-(36), Section 4, Supplementary Materials).We transformed this model to a second order differential equation form (for details, see Section 4.2, Supplementary Materials) that reads: where G(u) is the gating function of the slow gating variable n (Equation (55), Section 4.2, Supplementary Materials).Surprisingly, the transformed equation has almost the same form of the cubic BT normal form: a cubic-like force and a quadratic-like friction (see Figure 4), but also a term proportional to u2 .In Figure 4, we showed the simultaneous simulation of the full Morris-Lecar model, the Morris-Lecar model without the term proportional to u2 and a reduced version, where we also removed the u term (to test the effect of the singularity in the friction term at u = u 1 ).The force term in all of the equations is the same, and interestingly, it is proportional to the injected current minus the I-V curve, which will vanish when it is evaluated at a fixed point.Figure 4 shows the clear cubic-like shape of the force and also the quadratic shape, in the relevant region for the dynamics, of the friction.The quadratic shape of the friction is present in both cases, with and without the − I− f (u) u−u 1 u term, with slight differences between the two friction curves.We found the same global bifurcation scenarios as in the cubic BT normal form (i.e., saddle-node homoclinic, big homoclinic and saddle homoclinic bifurcation) in the three versions of the model.The dynamics does not change qualitatively, but quantitatively, with differences in the specific value of I when a global bifurcation

Conclusions
Using as a starting point a generic CB model, we looked for a critical point that encloses the diversity of local bifurcations observed in single neuron dynamics.In order to do that, we first constrained the general model to be in a critical point that possesses both local saddle-node and Andronov-Hopf bifurcations in its neighborhood.By using theoretical considerations, we concluded that this critical point is the BT bifurcation.In order to include the experimentally-observed curvatures of the I-V curve within the neighborhood of the critical point, we further constrained the model to be in the transition between one to three fixed points.This led us naturally to conclude that the BT normal form at this particular point must be cubic.Then, the dynamics of any CB model that is in some neighborhood of this particular BT bifurcation should be qualitatively captured by this equation.Our simulations showed that, in fact, it robustly displays much of the features observed in single neuron dynamics.Furthermore, we numerically found that the Morris-Lecar model, one of the most widely-used CB models, is in practice equivalent to the cubic BT normal form for realistic parameters values.These results might be explained by the fact that we construct the critical point imposing general features that are robustly observed among neurons.Additionally, by so doing, we obtained a critical point that correspond to a realistic parameter range for CB models.Hence, moderate displacements from this point in the space of parameters by different CB models, which represent different types of neurons, will still lay in the region where this normal form is qualitatively valid.This is consistent with the fact that this model cannot describe more unusual scenarios when, for example, a third, very long time scale conductance is introduced, and CB models display a dynamics, named bursting.This type of dynamics is generically well described by models of at least three equations [25].We believe that it is necessary to impose a further constraint on the critical point to arrive at a three-variable normal form that successfully reproduces this type of behavior.Recently, we gave some steps toward that direction [1].

Figure 1 .
Figure 1.(A) The I-V curve (i.e., function I ∞ (u * )) for the Hodgkin and Huxley model for three different values of the potassium conductance g K .A red dot indicates a fixed point in the dynamics.(B) The Bogdanov-Takens (BT) bifurcation in the Hodgkin and Huxley model.The 1+ τ 0 (u * ) ∑ N j=1 G j (u * ) function is shown in purple.The 1 − ∑ N j=1 G j (u * )τ j (u *) function is shown in green.The plots in black frames show the eigenvalues of the linear matrix L(u * ) in the complex plane.We omitted the greater in absolute value eigenvalue (far to the left) in order to better visualize the two critical ones.Tags 1-3 indicate the corresponding values of the variable u * .The plot in a red dashed frame shows the eigenvalues of the linear matrix L(u * ) in the complex plane when the model undergoes a BT bifurcation.

Figure 2 .
Figure 2.The left graphic shows in green a linear friction curve.The right graphic shows in red a quadratic friction curve with 0 < λ 2 .The arrow indicates the points where the friction changes the sign in both cases.In pale red, we indicate the region where the system undergoes an injection of energy 0 < d dt E(u, α) and in pale blue the region where the system undergoes a dissipation of energy d dt E(u, α) < 0.

Figure 4 .
Figure 4. Simultaneously simulated full Morris-Lecar model (green), Morris-Lecar without the term proportional to u2 (yellow) and the reduced Morris-Lecar (purple) for three different values of parameters.The force term ( I− f (u) τ(u) ) is the same for the three equations, and we plot it in the phase plane in blue.For the full Morris-Lecar model and Morris-Lecar without the term in u2 , the friction curve is displayed in red, for the reduced model in magenta.The top, middle and bottom figures show a saddle-node homoclinic, big homoclinic and saddle homoclinic bifurcation, respectively, before, in and after the bifurcation.The arrows show the direction of increasing of I.The same shapes and color codes in Figure 3 are used.
By doing the following nonsingular translation m = m ∞ (u; σ) + x, any generic CB model equation (i.e., Equation (1)) can be transformed into the equivalent equation: