From Squid to Mammals with the HH Model through the Nav Channels’ Half-Activation-Voltage Parameter

The model family analyzed in this work stems from the classical Hodgkin-Huxley model (HHM). for a single-compartment (space-clamp) and continuous variation of the voltage-gated sodium channels (Na v) half-activation-voltage parameter ΔV 1/2, which controls the window of sodium-influx currents. Unlike the baseline HHM, its parametric extension exhibits a richer multitude of dynamic regimes, such as multiple fixed points (FP’s), bi- and multi-stability (coexistence of FP’s and/or periodic orbits). Such diversity correlates with a number of functional properties of excitable neural tissue, such as the capacity or not to evoke an action potential (AP) from the resting state, by applying a minimal absolute rheobase current amplitude. The utility of the HHM rooted in the giant squid for the descriptions of the mammalian nervous system is of topical interest. We conclude that the model’s fundamental principles are still valid (up to using appropriate parameter values) for warmer-blooded species, without a pressing need for a substantial revision of the mathematical formulation. We demonstrate clearly that the continuous variation of the ΔV 1/2 parameter comes close to being equivalent with recent HHM ‘optimizations’. The neural dynamics phenomena described here are nontrivial. The model family analyzed in this work contains the classical HHM as a special case. The validity and applicability of the HHM to mammalian neurons can be achieved by picking the appropriate ΔV 1/2 parameter in a significantly broad range of values. For such large variations, in contrast to the classical HHM, the h and n gates’ dynamics may be uncoupled - i.e. the n gates may no longer be considered in mere linear correspondence to the h gates. ΔV 1/2 variation leads to a multitude of dynamic regimes—e.g. models with either 1 fixed point (FP) or with 3 FP’s. These may also coexist with stable and/or unstable periodic orbits. Hence, depending on the initial conditions, the system may behave as either purely excitable or as an oscillator. ΔV 1/2 variation leads to significant changes in the metabolic efficiency of an action potential (AP). Lower ΔV 1/2 values yield a larger range of AP response frequencies, and hence provide for more flexible neural coding. Such lower values also contribute to faster AP conduction velocities along neural fibers of otherwise comparable-diameter. The 3 FP case brings about an absolute rheobase current. In comparison in the classical HHM the rheobase current is only relative - i.e. excitability is lost after a finite amount of elapsed stimulation time. Lower ΔV 1/2 values translate in lower threshold currents from the resting state.


Introduction
The Hodgkin-Huxley (HH) model recently turned 60 [1,2]. It has been the only widely used prototype for microscopic single-compartment models [3]. Other popular models (e.g. the leaky integrate'n fire) do not contend for either detailed and type-specific descriptions of single neurons, or even less of their parts.
More recently, related to the detailed microscopic reconstructions, multi-compartmental neuronal models are built of electrically coupled HH-type compartments. The model has eloquently proven its quantitative capacity to convey an explanatory power to relatively simple neuron models, which account for more and more experimental findings [3]. Yet it has also been criticized, re-parameterized and subjected to multiple attempts to find it incomplete or inefficient [4,5]. Importantly, the HH model (Table 1 introduces all the commonly used abbreviations) simplicity was blamed, while critics sometimes omitted essential properties themselves [6]. The model's very applicability to mammalian neurons-with respect to metabolic requirements and flexibility of encoding, has been closely reexamined [2,7,8].
An unifying parametric framework is proposed in this work to systematically examine the impact of Na v subtypes' distributions on neuronal dynamics, and thence on excitability and refractoriness.
We systematically explore the influence of parameter variation on the fundamental biophysical properties of voltage-gated sodium (Na v ) channels, within an HH-type modeling framework. The Na v channels produce large and fast membrane currents essential in the generation and propagation of action potentials in excitable tissues. They contain a transmembrane alpha subunit which forms the channel pore. The latter subunit's genetic expression proved sufficient for the expression of whole functional channels of a given subtype [9]. Hence Na + channel nomenclature follows closely that of their alpha subunit (in this work we use the Na v1.X notation). With the help of gene engineering and selective expression of specific Na v1.X channels, significant experimental evidence has been accumulated on their (in)activation and localization properties [10][11][12].
Relatively little is known about the ways in which the different Na + channel subtypes are distributed and expressed toward functional axons in either a developmental or mature stage [13,14]. An exceptional wealth of evidence comes from epilepsy research [15][16][17][18][19][20], where Na + channel mutations have been associated with either gain-of-function or loss-of-function effects [19]-i.e. increased or decreased neuronal excitability in either excitatory or inhibitory populations (e.g. GABA interneurons or Purkinje cells). The Na v1.1 subtype, which is hypothesized to undergo such mutations, is also involved in an important developmental aspect. Namely, it gradually replaces the low threshold Na v1.3 subtype-which is only expressed during early development or excitable-tissue injury [21]. Interestingly, both the Na v1.1 and Na v1.3 subtypes are encoded on chromosome 2q24. It may also be tempting to speculate that the easily excitable Na v1.3 subtype is desirable during large neural network connectivity formation, but would lead to dynamic stability issues in an adult highly active and interconnected brain.
It is known that action-potentials (AP) are primarily initiated close to the axon hillock [22]. The cerebral cortex is densely packed with an estimated 50,000 neurons/mm 3 and at least 100 times as many neural processes [23]. A recent study [24] used HH models of two cortex-specific Na v1. 2 and Na v1. 6 subtypes. Na v1.2 sustained AP propagation, while Na v1.6 activation at lower membrane voltage (V) values contributed to AP initiation. The same Na v channel subtypes were then used in a multi-compartmental HH model of a neuron with a single dendrite, soma and axon morphology along a straight axis [25], which robustly interpreted and predicted the clinical effects of intra-cortical micro-stimulation (ICMS), following up on previous work on the subject [26]. It was found that AP initiation in the axon initial segment (AIS) was more likely and required lower stimulation current. Moreover, this was attributed to the higher density of the Na v1.6 subtype.
Benefitting from a detailed exploration of the literature on experimentally observed Na v channel properties in the central nervous system (CNS), we address the parameter V 1/2 controlling the membrane voltage V at which the Na + conductance attains its half-maximal value. This parameter has a direct impact on a number of fundamental properties. Some of these are straightforward to demonstrate from first biophysical principles. Others required subsequent bifurcation and phase-plane analysis, which are only enabled by models with just a few parameters. Nonetheless, this reductionist approach provides useful generalizations about key aspects under investigation-such as the HH model limits of metabolic efficiency or encoding. We do not claim that the model presented here is a definitive one, or that the different Na v subtypes are only due to simple shifts in V 1/2 . Nature has access to many complex molecular structures and there may be restrictions we are unaware of. Hence, overly simplistic models may fail to explain some observations. The proper model type depends on its goals and phenomenological scale [3]. Nonetheless, simple models can generate useful predictions and are easy to incorporate as building blocks in novel more elaborate paradigms. Thus, we use our parametric framework-within a physiological V 1/2 range-to address metabolic energy savings and other desirable properties such as fast AP transmission or large frequency-encoding capacity. We provide a comprehensive excitability-dynamics description of the Na v subtypes similar to the ones expressed in the mammalian brain (Na v1.X with X = 1, 2, or 6). Such analysis over a range of continuous V 1/2 variations seems in order, especially given that the V 1/2 values reported by different authors vary significantly [2, 10-12, 15, 17-21, 27-29]. assuming that one (and only) given Na v1.X subtype is expressed in the host cell, on which whole-cell voltage clamping techniques are applied. Very recent experimental and computational work on the molecular and functional effects of brain trauma points in quite a similar direction. Namely that Na v channels undergo irreversible hyperpolarizing V 1/2 shifts [30,31]. This may for example lead to ectopic excitability and propagation in damaged axons [32].
The presentation is laid out as follows: The next section (Methods) presents the model details and key assumptions. In a structured form the Results section defines a fundamental ensemble of related key biophysical properties of the parameterized single-compartment HH model, such as metabolic efficiency and excitability. The theoretical analysis presents the bifurcation structure as a function of the ΔV 1/2 parameter, and as a function of the constant (and hence known as bias) current I bias , applied intra-cellularly, which is another parameter receiving special attention in the relevant literature The Discussion provides a succinct summary of the most important findings, as well as a generalization of some results to multi-compartmental models.
To facilitate reading and provide for a wider audience, the most technical sections on automaticity and on the codimension-2 bifurcation structure in the ΔV 1/2 × I bias parameter plane are presented separately as Supplementary Results.

Single-node model
In this work, we explore the nonlinear dynamics for the following 4D system of ordinary differential equations (ODE): where I s stands for the stimulation current injected into the modeled compartment and is the total ionic current. The channel-gate state variables of the HH model m, h and n are described by: Matlab was used to numerically solve the model Eq (1) and AUTO [33] -for fixed point (FP) or periodic orbit (PO) continuation and bifurcation analysis. We use the single-compartment Hodgkin-Huxley model as described in [34]. The same model is also used in [24,25] as a building block. There we noticed an interesting pattern in the description of the Na v1.2 and Na v1.6 channel subtypes. For both the m and h gating variables, all the related V 1/2 parameters differed by exactly 13mV.
The Na v1.6 and Na v1.2 channel types can be seen as two instances of a generalized Na + ionic current model, in which the V 1/2 parameter of both the m and h gate variables is controlled by the parameter ΔV 1/2 with values of 0 mV and 13 mV, respectively (Fig 1 and Tables 2, 3 and 4; see also the gate-state dynamics maths in Box 1.

Model extensions
A Mixture model. provides a comparative perspective for the effects of parameter variation, and specifically concerning the effects of ΔV 1/2 to those of channel density by types and numbers (i.e. g Na variation). In this model extension, a fraction P 1 of all the Na + ion channels are of the Na v1.6 type (i.e. ΔV 1/2 = 0 mV), while the remaining Na + channels-i.e. a fraction 1 − P, are of a second different type of channel with a chosen ΔV 1/2 6 ¼ 0. Thus: I Na;S ½P; DV 1=2 ðtÞ ¼ ð1 À PÞI Na ½DV 1=2 ðtÞ þ PI Na ½0ðtÞ ð2Þ Such mixture models may provide predictions of the electrophysiological properties for neural processes with a particular experimentally observed distribution of voltage-gated sodium channel subtypes (e.g. as in [10,12,14,27,35]). From a compartment to a cell. The single-node model can be generalized to a multiplecompartment model of an entire axon. In models of higher dimensionality, the spatial distribution (relative proportions) of ionic current types will correspond to better or worse conditions for AP initiation and propagation.
To estimate the effect of ΔV 1/2 on AP propagation velocity (through the Na v currents), homogeneous cables of 100 identical compartments-each of length 25 μm (yielding a total cable length of 2.5 mm) were simulated (see also [36] for a detailed description of this type of cable models). The cable was subject to the no-flux boundary condition on one end, and the single compartment on the opposite end was stimulated (T STIM = 100 μs or 1 ms). For each value of the ΔV 1/2 parameter, the resting threshold was identified. An AP was reliably evoked and propagated to the very end of the cable. A propagating AP was signalled by small deviations of propagation velocity (range < 30% of mean value) in the section occupying the middle 50% of the total cable length.
To avoid overestimating the required stimulation current thresholds, the total duration of simulation was chosen as T STIM + 1.5 ms. However, this also provided for a little extra variability in the identified I THR , due to phenomena resulting from the cable's finite length and no-flux boundary conditions. Namely, lower (near-threshold) stimulation current latently reaches and  depolarizes compartments down the cable. As a result, when the AP is finally evoked it propagates faster than in a significantly "supra-threshold" case. Hence, to avoid excess variability in the estimated propagation velocity, a safety factor was applied: I STIM ¼ 1:5 ÂÎ THR .

Results
Our results demonstrate that the continuous variation of the ΔV 1/2 parameter comes close to being equivalent with recent HHM 'optimizations'. The dynamics of the parametric family exhibit a rich multitude of properties, which we have grouped in 4 main categories, namely: excitability, efficiency, frequency-encoding range and conduction velocity. The former three are discussed in a single-compartment (space-clamp) context. Conduction velocity and wholeneuron excitability are presented within a multi-compartment model extension.

Effects of ΔV 1/2 variation on neuronal excitability
Can the parameterized HH models' exhibit wide variation in their excitability? It is well known that the HHM's excitability to brief stimulation pulses is fully dependent on the fast activation of the m gates versus a slower inactivation of the h and activation of the n gates (Fig 1A and 1B). This can be demonstrated by a simplified model in which m is assumed to reach instantaneously m 1 (V) while h and n remain at their resting values h rest and n rest . Upon very brief I s duration this is a fair approximation of the system's behavior (see Ch.3 in [37]). Stimulating from rest [V rest , h rest , n rest ], the system dynamics is reduced to: C m dV dt ¼ I s À I ion ðV; m 1 ðVÞ; h rest ; n rest Þ ð3Þ Notes:: *1: The temperature-dependence factor. Please see the first note in Box 1. *2: Gate-state dynamics parameters of the K + or Na þ v1:X channels. Actual channel opening/closing rates are given by Boltzmann-like functions-for details see the second note in Box 1. *3: The Na þ v1:X channels' inactivating h gates have an asymptotic state, which is independent of the h gates' opening/closing rates-see the second and third notes in Box 1. For all ΔV 1/2 , the above approximate ionic current-denoted I ion,0 (V)-a function of just V, has 3 zero-crossings: V rest , V thr (see Fig 1D) and V up (not shown). The LOCAL minimum −I rh is reached for a V value located between V rest and V thr . While V rest and V up remain virtually unchanged for different values of ΔV 1/2 , V thr moves toward V rest and I rh diminishes as ΔV 1/2 is decreased. In this 1D approximation, one can see that the membrane will depolarize further toward V up if at the end of the stimulation V > V thr . The latter is possible only if I s > I rh . Hence lower ΔV 1/2 increases the excitability by reducing both V thr − V rest and I rh .
Strength-duration (SD) curves for the ΔV 1/2 -family, computed through simulations of the full HH dynamics, and for the whole range of For the short (e.g. T STIM = 100 μs) the increase of I THR with ΔV 1/2 is close-to-linear, while for the long (e.g. T STIM = 2000 μs) the increase is steeper than linear.
Toward approximate analytic SD curves, the model is simplified by assuming that all gate variables follow their asymptotic values, i.e.: This model is valid for low-amplitude stimulation of very long duration. I ion,1 (V) may have one or three zero crossings ( Fig 1C). For the one-zero-crossing case, the very existence of a threshold is a transient phenomenon, i.e. a finite V thr no longer exists, when the h and n gates Membrane capacitance, in μF/cm 2 : Maximum (*2) conductances, in mS/cm 2 : g Na Na + conductance 300 Currents, in μA/cm 2 : Notes: *1: Typical values are for the Na v1.6 model, [25]; see also Table 3 *2: These are dependent on (grow with) temperature, the values listed are for T = 23°C *3: Membrane voltage is either at its resting value V rest ; is depolarized (grows due to stimulation and/or activated sodium Na + ion channels); is repolarized (decays back to V rest , due to the potassium K + ion channels) *4: Ionic currents depend on both the membrane voltage and the dynamic state of the ion channels' gates.
See Table 3. , where T 0 = 23°C. Ã 2: Opening/closing rate of gates y in the K + or Na þ v1:X channels have Boltzmann-like templates (Á y subscripts dropped to simplify notation): where w = V m − V 1/2 , and the rates change in function of w: When both the opening and closing rates share the same V 1/2 parameter (e.g. the m and n gates, Table 2), the corresponding V-dependent asymptotic gate state and time constants are: y 0 1 ðwÞ is reciprocally dependent on k, and its maximum is obtained at w Ã = −kln(a/b) (hence with a = b, w Ã = 0). Hence dy 1 (w Ã )/dw is determined by the V 1/2 parameter, and y 1 (w Ã ) = 1/2-regardless of the actual a, b, k values.
Also for V m % V 1/2 , τ(V) = τ(w) attains its maximum (τ(w Ã )%τ(0) = 1/k(a + b)), since: , which vanishes when w = 0. Ã 3: Inactivating h gates of the Na þ v1:X channels have the asymptotic state respectively close and open. With three zero crossings −I ion,1 (V) also has a local minimum-a situation qualitatively similar to the I ion,0 (V) case above. Hence, just as with the simplified model of Eq (3), further depolarization of V toward V up,1 may be expected when V is beyond V thr,1 at the end of the stimulation. In this case I rh is precisely the value of the rheobase current.
Through I ion,1 (V) of Eq (4), a lower bound of the threshold current-the absolute rheobase -can be estimated.
Hence, SD[V shift ] may be estimated for relatively short durations (T STIM ! 0) or long durations (T STIM ! 1). In the latter case, this is only possible for models like Na v1.6 (ΔV 1/2 = 0), but not for models like Na v1.2 . This is one of the key properties pointed at in this work.
To evoke an AP, according to Eq (1) the stimulation current I s needs to overcome the dominating opposing ionic currents I ion . In search for patterns of lowest amplitude or highest energy efficiency, I s may just slightly outweigh I ion and still successfully trigger AP's [36]. However, this is highly dependent also on stimulation duration and Na V class.
Importantly, for the Na V1.6 subtype, but not for the Na V1.2 subtype-as demonstrated through the simplified model Eq (4)-a threshold value of the membrane voltage exists for any stimulation duration. Hence an AP can be evoked by applying almost as little current as the absolute rheobase current of a sufficient duration (Fig 2A).
Refractoriness (absolute or relative) is the other face of excitability. In a standard procedure to assess it, one applies a second stimulus S2 at different times along a test action potential AP S1 , induced by an initial supra-threshold stimulus S1. The minimal current I trh,S2 , needed for another active response, is then found for each S1 − S2 coupling time, or else the system is diagnosed as absolutely refractory-if current amplitude lies beyond some acceptable limit. I trh,S2 is Parameter-Dependent Na v Ion-Channel and Neural Dynamics a function of t = t S2 − t S1 coupling interval (assume t S1 = 0 for simplicity), as well as of the S2 duration T s,S2 . Whether it also depends on the parameters of the S1 stimulus is to be investigated. Alternatively, refractoriness can be studied by a one-dimensional approximation similar to Eq (3).
in which h 0 = h(V 0 ) and n 0 = n(V 0 ) are the values of the gate variables at each point V 0 along the AP S1 . As for stimulations from the resting state, a necessary (but not sufficient) condition for S2 to produce an AP is that Eq (5) has 3 FP's. When V thr does not exist, the system is absolutely refractory. Otherwise, a minimum stimulus current is needed.
We consider the case where AP S1 is evoked by a stimulus applied from the resting state. Both V thr (t)-calculated from Eq (5), and I thr,S2 (t) depend on the AP S1 , which in turn varies with T s,S1 (and the corresponding I s,S1 ). However, both functions become invariant once the time axis is aligned to t up -the time of maximum dV/dt derivative during the upstroke. Different I s,S1 and T s,S1 may modify the latency of AP S1 , but almost none of its time course beyond the upstroke. This invariance is maintained as long as T s,S1 does not become too long. Conversely, when T s,S1 ! 1, I s,S1 becomes what is known as bias current. It may induce automatic firing, which in turn precludes the use of the S1/S2 protocol (see also the Supplementary Results subsection on automaticity).
The invariance of the AP S1 time course means that a single I thr,S2 (t − t up ) curve (excitability is a function on the recovery time t − t up ) can be used to capture the system behavior for each T s,S2 and ΔV 1/2 . Rabinovitch et al. [38] referred to this invariant trajectory as hidden structure (HS) and proposed a method that extends to excitable systems the concept of the phase transition curve (PTC), widely used to analyze the forcing of non-linear oscillators [39,40].
Panel B of Fig 2 shows the variation of I thr,S2 as a function of t − t up and ΔV 1/2 for T STIM = 0.1 ms. Decreasing ΔV 1/2 slightly shortens the absolute refractory period. For any given t − t up , it also reduces the current needed to get a response. For any t − t up , there is a quasi-linear increase of the I thr,S2 as a function of ΔV 1/2 . For t − t up > % 2.5 ms, where excitability regained in all model Na V subtype cases, a higher ΔV 1/2 may mean up to a two-fold increase of the threshold current.

The metabolic efficiency of the HHM[ΔV 1/2 ] family
Can the parameterized HH models exhibit higher metabolic efficiency?
The Na + ions that enter the cell during an action potential (AP) will have to be returned back to the extra-cellular space. This is done by the Na + /K + -ATPase enzyme pump, while simultaneously bringing the repolarizing-phase K + ions back to maintain the resting potential. It binds 3 intracellular Na + and 2 extracellular K + ions. One ATP molecule is hydrolyzed, leading to phosphorylation and a conformational change in the pump, which exposes the respective ions to alternate sides of the cell's membrane. Thus the Na + /K + pump can be responsible for up to 70% of the neurons' energy expenditure.
Hence, the metabolic cost for firing a single AP is proportional to the Na + ionic charge that enters the cell during each AP, and also depends on the overlap of opposing Na + /K + ionic currents [41].
Such AP properties as the I Na /I K currents interplay and metabolic efficiency as a function of ΔV 1/2 are illustrated in Fig 3. 8ΔV 1/2 the AP is obtained via 1.5× the resting I THR for duration T STIM = 100 μs. As observed, larger ΔV 1/2 values translate into lower excitability: higher threshold potential V THR and threshold stimulation currents I THR . The very large latencies for ΔV 1/2 = 0 and -5 mV are due to very low threshold values I THR , which remain low even after a 50% "safety-factor" increase. The thick color traces in Panel A of the figure show the −I Na current profiles, aligned to the AP upstroke (time zero). The Na + current has a local peak just before the AP upstroke. there is a complete overlap of the Na + and K + currents (shown by color and dashed-black traces respectively in Fig 3B). The close I Na /I K interplay here gives an electrophysiological and physical perspective onto the nonlinear dynamics concept of the "slow manifold" [37,42,43].
As Fohlmeister and colleagues have observed [7,35], in the classical HH model the large overlapping Na + and K + currents lead to considerable (and rather 'suboptimal') metabolic energy expenditure.
However, smaller ΔV 1/2 values reduce both the Na + and K + currents through the AP's repolarization stage. First, lower ΔV 1/2 implies lower V THR . Second, the AP upstroke amplitude V up is invariant to ΔV 1/2 and hence h 1 (V up ) is significantly more closed for smaller ΔV 1/2 values (see Fig 1A).
The resting FP V rest remains almost identical in the Na v1.2 and Na v1.6 models (-77.03 vs -77.01 mV), but nevertheless corresponds to quite different values of h 1 (V rest ) (0.96 vs 0.76). Hence, the Na v1.6 h gates will close more during the AP upstroke, and will recover later and to a lesser value during the post-upstroke re-polarization. Furthermore, the Na v1.6 model has two additional depolarized FP (diamond and square in Fig 1C-denoted V thr,1 and V up,1 ). Compare the equilibria of the gate state h on Fig 4 for the typical models of Na v1.6 (ΔV 1/2 = 0), and Na v1.2 (ΔV 1/2 = 13). The latter figure contains the big picture of our bifurcation analysis (presented next in more detail), namely that the ΔV 1/2 family contains as a superset a wide variety of models with different properties, in which models behaving exactly like the classical HHM archetype are encountered over only a sub-interval of the ΔV 1/2 parameter full range of variation.
Efficiency is estimated by computing the overall Na + ions' charge Q Na (in micro-Coulomb's per cm 2 ), which is transferred during a single AP ( [41,44], T STIM = 100 μs). The integration is carried over the whole time span (30 ms) of the simulation (from AP triggering, incl. latency to the return to rest). The shape of the AP and hence Q Na are mostly invariant as long as I STIM ) I THR . A lower ΔV 1/2 clearly associates with a lower estimate of metabolic cost (see Fig 3, Panel C).
Finally, is this effect trivial? Can it be easily deduced by looking at the HHM Eq (1) alone? Here we argue that this is not the case: Without a significant effect on AP amplitude or duration (Fig 3A and 3B), the overall shape of the AP changes just slightly and a significant gain in metabolic energy expenditure is achieved by an earlier closing of the voltage-gated Na + channels. Since V THR is lower in this case, this earlier Na + inactivation is matched to a significantly lesser activation of the voltagegated K + channels. Parameter-Dependent Na v Ion-Channel and Neural Dynamics Bifurcation structure as a function of the ΔV 1/2 parameter The resting potential (V rest ) of the system is a fixed point (FP) of the system set by: It is located at the intersection of the monotonically increasing outward currents I K,1 (V) + I leak (V) with the bell-shaped ΔV 1/2 -dependent inward current I Na,1 (Fig 1C). The effect of ΔV 1/2 variation is a linear shift of the window-current I Na,1 [ΔV 1/2 ](V) toward lower or higher V values. This may create additional FP's. The number, location and stability of the FP's is bound to play a significant role in the dynamic properties of the system upon stimulation and re-polarization. The existence of the two unstable fixed points also changes the organization of the trajectories in the phase plane. As mentioned in the presentation of the bifurcation diagram below, V thr,1 has a single unstable direction. Perturbations along this direction (either positive or negative) produce two heteroclinic trajectories ending at V rest , corresponding to a passive and an active depolarization (Fig 5, the red and green trajectories). The latter is the limiting case of an action potential that would be obtained by a stimulus slightly beyond the rheobase current. Another set of trajectories connect the most depolarized fixed point (unstable focus) V up,1 to the stable resting state. All these invariant trajectories and their associated stable manifold do not exist for models like that of the Na v1.2 subtype. They add constraints onto the system dynamics since they cannot be crossed by any other trajectory. As we shall see below, the two additional FP's also affect the nature of the automatic regimes upon sustained stimulation. Fig 4A shows the bifurcation structure diagram (BD) of the of the 4D space-clamp model as a function of the ΔV 1/2 parameter. In addition to the resting point (the lower branch on the BD, see the "stable FP" arrow), a second depolarized FP (corresponding to V thr in panel A) appears through a saddle node bifurcation at the limit point LP1: ΔV 1/2 = 2.94 mV. From LP2: ΔV 1/2 = -9.5 mV, up to LP1, the S-shaped BD FP curve has three coexisting branches. Beyond LP2, a single depolarized FP remains. FP stability is also accounted for in Fig 4A. The middle branch (LP1 to LP2) is half-stable (see Table 5 and [43]) due to a single positive real eigenvalue. The lower FP branch is stable with four real negative eigenvalues, except for a tiny interval close to LP2 (there is a second subcritical Hopf bifurcation HB2 located 0.005 mV beyond LP2, not shown). The uppermost branch is unstable from LP1 down to the Hopf bifurcation point (HB: ΔV 1/2 = -11.05 mV), where it recovers stability. From LP1 to HB the top branch has either two positive real (ΔV 1/2 2 [1.7 mV, LP1]), or two positive-real complex eigenvalues (ΔV 1/2 2 [HB, 1.7 mV]).
The HB is subcritical (see Table 5) and produces a branch of unstable PO's connecting through a saddle-node of cycles (SNC) bifurcation to a high-amplitude stable PO at ΔV 1/2 = -35.5 mV. This stable PO branch disappears slightly beyond LP2 (around ΔV 1/2 = -9.16 mV) through a second SNC bifurcation (see Table 5, as explained in the Supplementary Resultssubsection on automaticity). In this work, the analysis is mostly restricted to ΔV 1/2 > LP2.
There are two distinct model sub-classes by the parameter value range: M1 ΔV 1/2 > LP1. The corresponding dynamic system has a single stable resting (hyperpolarized) FP, e.g. the Na v1.2 channel subtype.

Term Description
FP stability A given FP is stable if all the eigenvalues of the associated Jacobian have non-negative real parts; an FP is called half-stable (e.g. in [43]) when the above condition is not met for some of the eigenvalues -in that case the associated eigen-directions diverge from the FP

Codimension 1 bifurcations
Saddle-node, or fold, or limit point (LP) Two fixed points appear or disappear coalescing onto a single half-stable FP-having stable, unstable and neutral (zero-real-part) eigenvalues, which in its turn disappears Transcritical bifurcation A given FP may change stability with parameters variation standard mechanism, related to the SN Hopf bifurcation (HB) Supercritical Subcritical a given PO's amplitude decreases until the PO is reduced to a point and disappears; the PO's period approaches a finite limit (given by the purely imaginary eigenvalues) for the critical parameter value the HB-related PO appears at the HB and its amplitude increases gradually, coexisting with the former FP which has become unstable. within a parameter range the HB-related PO's-a stable and an unstable of lower-amplitude, coexist with a stable FP, the parameter-dependent system trajectory jumps to a distant attractor, and hysteresis-like phenomena occur as the parameter is either increased or decreased. offers an alternative way to create or suppress the folding and hence the appearance of the middle and topmost FP branches. Likewise, a variation of the relative proportion of different channel subtypes in a mixture model is equivalent to nominal conductivity variations. Hence, an estimate of the g Na value yielding a transition from one to 3 FP's is useful to understand the dynamics of the mixed model. From Fig 1C, and with ΔV 1/2 = 0 mV, g Na has to be decreased by a third to suppress the BD from folding, but with ΔV 1/2 = 13 mV, it has to be increased fourfold to create the fold.
Except for a small interval of ΔV 1/2 2 [HB2, -9.16 mV] where a stable FP and a stable PO coexist, the two model sub-classes necessarily return to their resting state after a sub-theshold stimulus. However, the supra-threshold dynamics (leading to an AP) will be dramatically influenced by the existence or not of the 2 additional FP's (see the Supplementary Results' subsection on automaticity).

Lower ΔV 1/2 brings about a large frequency-encoding range
Can the parameterized classical HH model exhibit the low-frequency automatic firing needed to account for certain regimes observed experimentally in mammalian species?
Before we answer this question, let us introduce another preliminary bifurcation analysis result. It is known [37], that in the HH model periodic firing may be obtained by injecting a constant I bias -known as bias current within a certain range (see also Fig A in S1 File). In the literature one can frequently encounter bifurcation analysis of the HHM with respect to I bias .
Panels A and C in Fig 6 result from exactly such bifurcation analysis for the typical models of Na v1.6 (ΔV 1/2 = 0), and Na v1.2 (ΔV 1/2 = 13). Here a higher ΔV 1/2 translates into incapacity for low-frequency automatic firing and a narrow overall frequency-encoding range (see Panel B in Fig 4), consistently with earlier observations [7].
Panel B in Fig 6 illustrates the period length of automatic firing as a function of I bias . In the Na v1.6 case (see also Panel C), stable automatic firing (blue dots in Fig 6B) is lost through a cycle saddle node, which accounts for the finitely increased length of the related cycle period, which is a desirable property of the parametric model family. On the other hand, the unstable cycle (green dots in Fig 6B) disappears into an homoclinic bifurcation, when I bias tends to the forementioned absolute rheobase value. As is well known, homoclinic bifurcations may be associated with infinite increase of the related period. Fig 7 illustrates the limits of automatic firing and of periodic stimulation maintaining the 1:1 response pattern (also known as pacing) as a function of ΔV 1/2 . In Panel A of the figure, the periodically-applied stimulation current I STIM required for an 1:1 response-at the given minimal pacing period τ 11 , depends on pacing frequency. Clearly, I STIM ! I THR (the resting threshold, see also Fig 2), as τ 11 ! 1. Notice that, with near-threshold current values, higher ΔV 1/2 allow periodic stimulation at higher frequency.
However, the situation is reversed for automatic regimes. The period of automatic firing is a monotonously decreasing function of I bias with a minimum well below 5 ms. This means that some (less excitable) models may not have automatic firing at a frequencies below 200Hz. Importantly, this is a monotonicity property specific to the model at hand, contrasting it to other models (e.g. the Fitzhugh-Nagumo model)-(see Fig 7B-itself, a simplified presentation of Fig D.B in S1 File. The full detail of periodic firing variability and bifurcation structure as a function of ΔV 1/2 and I bias are provided as Supplementary Results).
How about the maximum period that can be achieved? Applying the same monotonicity property for the period, it would be longest for the lowest I bias for which automatic firing is still possible. Lower ΔV 1/2 provide for such regimes with progressively lower firing frequency (longer inter-AP periods). As seen in Panel B of Fig 7, a very reasonable ΔV 1/2 parameter variation translates into a five-fold variation of period lengths. Thus, the model family is capable of handling automatic firing at a frequency less than 20Hz.
Hence lower-ΔV 1/2 models bring about a significantly larger frequency-coding range, which in addition requires less metabolic resources.

On Hogkin's maximum-velocity hypothesis
Can the parameterized HH models conduct AP's faster?
The effect of ΔV 1/2 on AP propagation velocity (through the Na v currents), was studied in a multiple-compartment model as described in the Methods. To reliably evoke propagating AP's for all values of the ΔV 1/2 parameter, the resting stimulation-current thresholds I THR were identified.  Parameter-Dependent Na v Ion-Channel and Neural Dynamics (ΔV 1/2 = 2 mV, please see Methods for a detailed interpretation of the obtained voltage-distribution patterns in either case).
According to Hogkin's maximum-velocity hypothesis, the ion channels' ion-types and density evolved through natural selection to maximize the AP conduction velocity. However, it has been argued that gating current limits the increase of Na + channel density as a means to gain in conduction velocity. Moreover, the computed optimal density is more than twice higher than in the real squid, and with capacitance reduced by neglecting gating current, optimal density is even higher [41].
A summary of resting threshold I THR and mid-axon mean AP conduction velocity as a function of ΔV 1/2 is presented in Fig 9. Notice the clear linear trend of increasing I THR with ΔV 1/2 for both T STIM cases. Despite the large estimation variability (as illustrated and discussed in Fig  8) a significant trend of deccelerating AP propagation with ΔV 1/2 can be acknowledged.
Hence, evolution may play another card to solve many a problem at once-avoid the overhead of having too many channels, while simultaneously gaining in AP transmission speed, and even in efficiency (as just shown above). All of it by a simple shift of the sodium current window along the membrane potential dimension.

Mixture model
From the preceding analysis, the single-channel model has 3 FP's for ΔV 1/2 < LP1 (Fig 4A). This section analyzes mixed models with sodium current. For ease of reference Eq (2) is Parameter-Dependent Na v Ion-Channel and Neural Dynamics  Fig 8. A generalization perspective of ΔV 1/2 variation effects-from a compartment to a cell. A When I STIM is lower, it latently reaches and depolarizes compartments down the cable. This is clearly visible from the voltage traces, which are depolarized up to halfway down. As a result, when the AP is finally evoked it propagates faster. B "supra-threshold" case with a safety factor applied I STIM ¼ 1:5 ÂÎ THR . doi:10.1371/journal.pone.0143570.g008 recalled from the Methods section: The B.D. on Fig 10A shows that, for ΔV 1/2 = 13 mV, the 2 upper FP branches appear at saddlenode bifurcation (see Table 5) for P % 0.648, which is close to the minimum relative conductivity (g Na /g K = 0.663) needed for the single-channel model with ΔV 1/2 = 0 to enter the M3 class. In the voltage range of peak I Na,1 [ΔV 1/2 = 0](V), the contribution of I Na,1 [ΔV 1/2 ](V) is so low that it does not change much the total I Na,1 profile. Hence, the upper FP's rely almost exclusively on the of the M3 current contribution. In codimension-2 (see Fig 10C), as ΔV 1/2 decreases, the contribution of I Na,1 [ΔV 1/2 ](V) becomes larger, thereby reducing the proportion of M3 current needed to preserve the 3 FP's. In Fig 10C, the required proportion P of Na v1.6 channels (recall the Methods section) decreases abruptly as ΔV 1/2 < 5 mV. Clearly, the latter is due to the intuitive fact that such ΔV 1/2 values mean that the second channel type is itself tending to the M3 class-a tendency which is completely realized with ΔV 1/2 < LP1.
The relative resting threshold κ(P) = I THR (P, ΔV 1/2 )/I THR (1,ΔV 1/2 ) of the mixture model is a function of the P parameter (T STIM = 0.1 ms). There is an approximately twofold difference between the thresholds of the M3 and M1 single-channel models (Fig 10B, see also Fig 2), but more than 70% of this difference vanishes for P as low as 0.2. This could be understood by considering the 1D resting version of the mixed model-see Eq (3) with the m gates of the two channel subtypes at their saturation values, and all other gates at their resting values. As in Fig 10. Mixture models. A: Red trace: Existence/creation P parameter-range for the 2 additional fixed points (saddle-node in the middle branch, and unstable center/focus in the top branch of the B.D.'s) With ΔV 1/2 = 13 mV, P ! 0.648, which is very close to the minimum relative Na/K conductivity ratio (0.663) needed for 3 FP's in the single-Na v -subtype model with ΔV 1/2 = 0 mV. In the voltage range of peak M3 current, the contribution from the Na v1.2 channels is low and does not change much the total current profile. Hence, the upper branch of FP's is due the contribution of the Na v1.6 channels. B: Resting threshold of the mixed models as a function of the P parameter. The figure shows κ(P) = I THR (P, ΔV 1/2 )/I THR (1,ΔV 1/2 ) of the mixed model for T STIM = 100 μs. Notice that there is an approximately twofold difference between the thresholds of the M3 and M1 single-Na vsubtype models More than 60% of this difference vanishes for P as low as 0.1. The difference becomes minimal for ΔV 1/2 < 4 mV (see also Panel C) or P ! 0.1. C: As the ΔV 1/2 parameter is decreased (from 13 mV down to 0), the M1 contribution becomes larger, thereby reducing the proportion of M3 current needed for 3FP-dynamics. Notice that the proportion of Na v1.6 channels may decrease significantly when ΔV 1/2 < 4 mV. Importantly, this is also about where the transition from M1 to M3 occurs (ΔV 1/2 = LP 1 % -2.94 mV. doi:10.1371/journal.pone.0143570.g010 section on neuronal excitability, both I thr and V thr − V rest provide measure of excitability. Both I thr (Fig 10B) and V thr − V rest (data not shown) have decreased by more that 70% at P = 0.2.

A generalized model
The model family analyzed in this work stems from the classical Hodgkin-Huxley model (HHM). The latter has been criticized and its very applicability to mammalian neurons has been questioned. Various related HHM optimization directions have been suggested.
Here, we argue that: • The model family analyzed in this work contains the classical HHM as a special case • That the validity and applicability of the HHM to mammalian neurons can be achieved simply by picking the ΔV 1/2 parameter value in an appropriate range.
With respect to the first claim above, from Fig 1 the knowledgeable reader can rapidly see that the quantitative (and hence the qualitative) picture may also be very unlike that of the classical HHM. Namely the interplay between the asymptotic states y 1 (V) for the sodium and potassium channels is far more involved: In contrast to the classical HHM, the h and n gates' dynamics are uncoupled. I.e. the n gates can no longer be considered in mere linear correspondence to the h gates as is often assumed when working with the HHM and its simplifications such as the Fitzhugh-Nagumo model [45][46][47].
In the latter n = %1 − h lead to a model of reduced dimensionality. Such linear correspondence would however be inadequate over the whole range of ΔV 1/2 parameter variations. For the model considered here, regressing the n(t) variation during an AP (data not shown) by a polynomial in 1 − h requires at least order 2. I.e. there is nonlinear (parabolic) relationship between the n and h dynamic variables during an AP. Moreover, around the resting state, the K + channels are fully closed (n r est = 0), which is quite unlike the classical HHM (where n r est % 0.2). Hence, the resting FP is fully determined by the leak parameters (conductance and Nernst potential). This explains the invariance of V r est with ΔV 1/2 .
A similar perspective stems also from Fig 4, where the classical HHM corresponds to a narrow range of values for the ΔV 1/2 parameter out of the significantly broader variation range considered by this work.
As to the second claim, in the Results section we dedicated paper space and work effort to demonstrate that indeed essential model family properties, in function of the given ΔV 1/2 parameter value, covary in similar directions as the HHM optimizations published in the literature (see also the Summary and Conclusions below).
An general discussion point is also that many of the observed neural dynamics phenomena are nontrivial-i.e. cannot be explained singlehandedly by the model equations. They are the complex outcome from the interplay of the Na V and K V channels, their properties and the related multitude of key parameter choices.

Parametric-model relevance in the light of previous work
How does our parametric analysis of the HH model compare to previous work? Where does it bring new insight and where is it similar? A body of previous work examined bifurcation structure in the Hodgkin-Huxley model. One arguably popular parameter has been I bias . Rinzel and Miller examined the I bias range for which stable and unstable periodic solutions exist and their findings are qualitatively similar to our own results [48].
The two most relevant other parametric models, which we believe provide parallels and grounds for contrasts are the codimension-2 analysis in the parametric plane I bias × E K [49,50] and the two-dimensional I Na,p + I K HH-type model [37] with its either high-or low-threshold I K (i.e. two E K levels). The latter model is equivalent in many respects to the well-known I Na,p + I K Morris and Lecar model [51].
The particular interest of these two alternative parameter models is in their choice of "mirror" parameters: The potassium's Nernst potential E K is manipulated toward changing the I K sign-switching point. This may have similar effects to ΔV 1/2 variation, which affects the position of the sodium current window relative (closer or farther) to the inversion point of the opposing K + current. However, the two parameters are also qualitatively different as the ΔV 1/2 variation affects the onset (and offset) of sodium current through the dynamics of Na + channel opening and closing gates.
Both [49,50] and [37] contain a detailed analysis of the related nonlinear dynamics structure, as well as introduce the "menagerie" of codimension-1 and codimension-2 bifurcation types. This provides an excellent baseline toward interpreting the results of our work. Finally, other previous work (e.g. [52]) examined the variational effects of nominal g K conductance, relative to nominal g Na . As demonstrated in the Results' section on mixture models, this may also be qualitatively similar to a ΔV 1/2 variation.
The key difference brought about by ΔV 1/2 variation is that it affects the interplay of sodium and potassium currents not only directly (like compatible E K or g K variation) but also indirectly through the V-dependent temporal dynamics of the HH gates.

Summary and Conclusions
This work is about bridging theory to practice in an attempt to go all the way from the model mathematics (nonlinear dynamics), through simulations and to the experimental practice.
Here we summarize the gist of our key findings.

Summary of ΔV 1/2 -related dynamic properties
In this work we explore the effects of Na + ion-channel properties on neuronal dynamics. Starting from two well-studied channel types ( [24,25]) we observed that a key difference between the Na v1.2 and Na v1.6 types is their half-activation-voltage parameter V 1/2 . Hence, we introduced a generalization parameter. This resulted in a parametric family of models exhibiting a continuous variation of sodium current-influx properties.
The introduction of the ΔV 1/2 parameter provides for a comprehensive analysis of neuronal excitability, refractoriness and periodic-stimulation dynamics for a space-clamp model. As seen, this analysis can also be extended further to multiple-compartment models of an entire neuron. In models of higher dimensionality, the spatial distribution of ionic current types and their relative proportions will yield specific spatio-temporal patterns, corresponding to better or worse conditions for AP initiation and propagation. The approach used in this paper can be generalized further to produce these profiles of threshold stimulation-field and stimulationcurrent values, that lead to successful AP initiation and propagation [36].
The following key effects were determined and explored: Continuous ΔV 1/2 variation leads to a bifurcation diagram ( Fig 4A) containing a multitude of dynamic regimes. One of the most noteworthy differences is that for ΔV 1/2 2 [LP2,LP1] 3 FP's coexist. These may also coexist with stable periodic orbits. This means that, depending on the initial conditions, the system may behave as either purely excitable or as an oscillator.
Absolute vs relative rheobase current: In the range LP2 < ΔV 1/2 < LP1 (the M3 class), the simplified model of Eq (5) can be applied, and −I ion,1 (V) has a local minimum −I rh between V rest and V THR,1 (Fig 1C and 1D). Hence an AP can be evoked by applying as little current as the absolute rheobase current I rh for a sufficient amount of time T STIM (Fig 2).
The mere FP count corresponds to nontrivial differences in dynamic organization and its complexity. E.g. supernormality of I THR,S2 as a function of t − t up , esp. for the lowest ΔV 1/2 (data not shown).
A complete picture of refractoriness is given by: • the application of the S1-S2 simulation protocol (Fig 7, panel A and Fig 2, panel B). Fig 7  presents the two faces of periodic AP firing in the model. Importantly, and as noted in the Results' subsection, the M1 model sub-class is not only less excitable. Its lowest periodic firing frequency does not get much below 100 Hz. Hence, at stimulation durations of about 10 ms the models of this sub-class will already behave as oscillators. We return to this point in the summary on 'relative' refractoriness below.
• the continuation (Fig A in S1 File) of the bifurcation diagram of Fig 4A, with respect to bias current and periodic regimes (see the Supplementary Results section on the codimension-2 bifurcation structure in the ΔV 1/2 × I bias parameter plane). Fig A in S1 File demonstrates transitions between pure excitability and oscillation determined by the I bias and ΔV 1/2 values. The automatic-firing region of I bias variation becomes gradually larger as ΔV 1/2 is increased from 0 to 13 mV (Fig A in S1 File). However the corresponding inter-spike period shows very little variation (see Panel B in Fig 4).
Absolute values of I STIM as a function t − t up and ΔV 1/2 (Fig 7, panel A) imply that "1:1" responses of lower-ΔV 1/2 models can be maintained for either shorter pacing periods or lower stimulation currents. However, maintaining 8ΔV 1/2 a similar proportion of I STIM relative to the resting threshold I THR,0 [ΔV 1/2 ](T STIM ) (see also Fig 2), Fig 7B predicts that the higher-ΔV 1/2 models would be the first to behave as oscillators as T STIM is increased.
'Relative' refractoriness may be the consequence of more complex dynamic organization. Thus low ΔV 1/2 means higher excitability and the existence of an absolute rheobase (Figs 1C and 2). However, the higher number of FP's implies more features and-hence, constraints in phase space (Fig 5).
This is intersects into the main 3D surface, which shows the minimum 1:1 pacing I STIM currents as a function of ΔV 1/2 and the pacing interval τ 11 (the black-yellow line in Fig 7A).
For I STIM close to the I THR value (for a given model), the M1 sub-class is paced at a faster frequency. Such higher 'relative' refractoriness is also clearly present in Fig 7B where the longest automatic-firing period is substantially shorter for ΔV 1/2 = 13 mV, than it is for ΔV 1/2 = 0 mV.
Sometimes this may not be desired (see the Discussion points about a large frequencyencoding range, above). However, when I STIM is much above the I THR value, the trend can be completely reversed, given that the M3 model sub-class is both more excitable and has a very large frequency range (Fig 7B and Fig D.B in S1 File). The higher 'relative' refractoriness of the M3 sub-class may be in part due to the h gates' dynamics, since the position of the resting FP (V rest ) does not vary much with ΔV 1/2 . Hence, the h gates' may need a longer-lasting recovery phase of the membrane to regain excitability (Fig 3A and 3B).
Simplified models such as the ones relying on the concept of "hidden structure" (a generalization of the PTC concept for non-linear oscillators, data not shown) capture the essential parts of the complex dynamics-such as the invariance of the evoked AP and the latency dependence on I STIM compared to I THR .
On the need for a methodology to study and predict neuronal dynamics determined by ion-channel distributions Nonlinear dynamics mechanisms may explain some experimental observations in a predictive way. A methodology to study the effect of ion-channel distributions on neuronal dynamics, excitability and refractoriness holds promise for experimental neuroscience practice (e.g. [53]). Hence, such modeling and analysis are worth pursuing. Functional electric stimulation is at the threshold of a new realm of possibilities. The paradigm is rapidly shifting from classical work ( [54,55]), as rapid stimulation (typically through trains of * 300 Hz) to evoke single AP's in isolated neurons is largely suboptimal.
Recent experimental ICMS work (e.g. [53]) confirmed the modeling predictions (e.g. [22,25,26]) that AP's are evoked almost exclusively from axonal targets, even if this results antidromic AP propagation. Due to refractoriness, an outcome of 'ad hoc' high-frequency stimulation trains may be depression instead of activation. A given electrode(s) geometry and relative position with respect to the targeted excitable tissue will imply a given Na v channel density (as the latter depends on the underlying cell(s) morphology). From our analysis, it may be inferred that particular positions and stimulation protocols may be appropriate for a given desired outcome.
Together with critical aspects such as the long-term compatibility, functional visual prosthetics require effective and reliable (1:1) stimulation of select (precisely targeted) neural subpopulations. Such stimulation not only uses optimally little energy (from an optimal-control point of view) and has minimum undesirable side-effects (e.g. tissue damage or involuntary saccades, [56]). Selective stimulation with lowest possible currents also means higher possibility for conveying information through visual percepts and specific visual features: e.g. oriented lines conveying shape or colored stimuli versus bright but color-less and feature-less phosphenes ( [57]). Finally, a very large impact on the type of possibilities is determined by the size, type and number of the stimulation electrodes that are used, as well as by the precision of their guidance to their neuronal targets. We demonstrated here that with reasonably high I STIM values the limit periods of 1:1 pacing may be as short as * 5 ms. Hence stimulation train frequency can be as high as 200 Hz. The latter result is consistent with the experimentally observed fact that past such frequencies the subjective visual percept reaches a plateau, beyond which it starts getting weaker. This is suggestive that the neuronal targets of stimulation will respond with activation ratios lower than unity to such 'overdrive'.
Theory needs to build bridges to practice-from modeling (mathematics, nonlinear dynamics, simulations) to the applications such as functional electric stimulation. This case is further illustrated by Fig F in S1 File. An important feature here is that, slightly before the two HB's indeed coalesce, their type changes from subcritical to supercritical, which is accompanied by the disappearance of two extra branches of stable/unstable PO's (see Fig F in S1 File and SNC5 in Fig B.A in S1 File).
Dashed-black lines These illustrate the parameter-range for the saddle-nodes LP1 and LP2 through which are created the 2 additional FP branches (saddle in the middle branch and unstable center/focus in the top branch of the B.D.'s, see Fig 4A and Fig C in S1 File). Here, there are also 2 related codimension-2 bifurcation phenomena. First, the two LP's coalesce in what is known as Cusp (see Table 5). Second, for low-enough ΔV 1/2 LP1 coalesces with HB1 through a codimension-2 bifurcation known as Takens-Bogdanov (see Table 5).
Red and blue lines the ΔV 1/2 × I bias range of stable/unstable PO's determined by the saddlenodes for cycles SNC1 and SNC2 (the red locus) and SNC3 and SNC4 (the blue locus, see also Fig B in S1 File) Notice that outside (below and over) the HB area, there is at least bi-stability of a stable fixed point with stable PO (or multi-stability between the stable FP and two stable PO's, see also