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 (Nav) half-activation-voltage parameter ΔV1/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 ΔV1/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 ΔV1/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. ΔV1/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. ΔV1/2 variation leads to significant changes in the metabolic efficiency of an action potential (AP). Lower ΔV1/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 ΔV1/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].

Table 1

Commonly used abbreviations.

Symbol	Description
0D	zero-dimensional, i.e. single-compartment or space clamp models; whose spatial extents are confined to a point
1D	cable-like, multi-compartment spatial structure; homo-morphic to line
2D etc.	two- or more dimensional, refers to the number of states that describe the excitable system’s dynamics
AP	Action potential
B.D.	bifurcation diagram
BVDP	the Bonhoeffer-Van der Pol oscillator-dynamics model; also known as the Fitzhugh-Nagumo model
ES	Electrical stimulation
FP	Fixed point of system dynamics → vanishing derivative(s)
HH or HHM	Hodgkin and Huxley’s [model of excitable membranes]
ML or MLM	Morris and Lecar’s [model of the barnacle giant muscle fiber]
ODE	Ordinary Differential equation; see also PDE
PDE	Differential equation involving partial derivatives; see also ODE
PO	Periodic orbit (or limit cycle)—Closed (starting and ending at the same point in phase space) dynamic trajectory The period of the PO may be finite or →∞. In the latter case it may be a hetero- or homo-clinic (starting and ending at either two distinct FP’s, or the same single half-stable FP, respectively)
PTC	phase transition curve
RGC	retinal ganglion cells
RHS	right-hand side
SD	strength-duration [curve]
S.T.	such that
W.R.T.	with respect to

An unifying parametric framework is proposed in this work to systematically examine the impact of Na 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 ) channels, within an HH-type modeling framework. The Na 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 notation). With the help of gene engineering and selective expression of specific Na channels, significant experimental evidence has been accumulated on their (in)activation and localization properties [10-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-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 subtype, which is hypothesized to undergo such mutations, is also involved in an important developmental aspect. Namely, it gradually replaces the low threshold Na subtype—which is only expressed during early development or excitable-tissue injury [21]. Interestingly, both the Na and Na subtypes are encoded on chromosome 2q24. It may also be tempting to speculate that the easily excitable Na 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 and Na subtypes. Na sustained AP propagation, while Na activation at lower membrane voltage (V) values contributed to AP initiation. The same Na 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 subtype.

Benefitting from a detailed exploration of the literature on experimentally observed Na 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 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 subtypes similar to the ones expressed in the mammalian brain (Na 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 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 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 , 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 parameter plane are presented separately as Supplementary Results.

Methods

Single-node model

In this work, we explore the nonlinear dynamics for the following 4D system of ordinary differential equations (ODE): where I 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 and Na channel subtypes. For both the m and h gating variables, all the related V 1/2 parameters differed by exactly 13mV.

From whole-cell expression of the Na subtype, V 1/2 = −29.2±1.8 mV and V 1/2 = −29.4±1.6 mV (as reported by [10] and [11] respectively). For the Na subtype, V 1/2 = −18.32 mV [29]. We also noticed that the cited modeling and experimentally observed data preserved a quite similar relative difference in V 1/2—between Na and Na . Furthermore, [29] demonstrated the remarkable effects of the pentapeptide QYNAD. At 200 μM concentration, the Na subtype V 1/2 was increased to −9.15 mV—an almost 9 mV shift toward decreased neuronal excitability. Since QYNAD (Gln-Tyr-Asn-Ala-Asp) is endogenously present in human cerebrospinal fluid (CSF) [29], such chemical interactions may also occur naturally and thence have significant effects on neural dynamics.

The Na and Na 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.

Fig 1

Gate variables as functions of the membrane voltage (V).

A: Asymptotic state y ∞(V). B: Time constants τ (V) The subscripts 1 and 2 of the m and h variables refer respectively to the Na and Na ion channels. Panels C and D illustrate the ionic currents interplay in the model as a function of membrane voltage V, computed for different example values of ΔV 1/2. The color traces use a color-code compatible to panels A, B for the two key Na channel currents; the latter plots are inverted to facilitate the interpretation of equilibria. C: I [ΔV 1/2](V), computed for the asymptotic state y ∞(V) of all gate variables The thin black trace shows the total re-polarization current I + I Fixed points (FP’s) locations are illustrated through the interplay of I (V) components D: I [ΔV 1/2](V)—computed for m = m ∞[ΔV 1/2](V) similarly to I (V), but for the reduced (1D) model at the resting state of the slower h and n gates, i.e. h = h = h ∞[ΔV 1/2](V ) and n = n = n ∞(V ), where V = the resting potential. The circle markers in both Panels C and D indicate the resting-state FP V . The diamond markers indicate the second AP-threshold FP V , which for I (V) (panel C) exists only for the M3-class models (see text and Fig 4). The gray trace is for ΔV 1/2 = -8.7 mV →LP2 (see Fig 4).

Table 2

The V 1/2 parameter values [mV], by ion-channel type—from [24, 25].

Rate parameter	Na v1.6		Na v1.2		K
	m	h	m	h	n
α	-41	-48	-28	-35	25
β	-41	-73	-28	-60	25
h ∞	-	-70	-	-57	-

Table 3

Gate-state dynamics parameters (see also Box 1—key notes on gate-state dynamics).

Notation	Variable description	Value
Temperature dependence:
Q10	Q 10 constant (*1)	2.3
K +: n-gate (*2)
a n	n-gate max opening rate	0.02
b n	n-gate min closing rate	0.002
k n	n-gate membrane-voltage-change constant k	9
Nav1.X+: m-gate (*2)
a m	m-gate max opening rate	0.182
b m	m-gate min closing rate	0.124
k m	m-gate voltage-change constant k	6
Nav1.X+: h-gate (*2,*3)
a h	h-gate max opening rate	0.024
b h	h-gate min closing rate	0.0091
k h	h-gate constant k (*3) used toward τ h only	5
k h∞	constant k (*3) for the asymptotic gate-state h ∞ only	6.2

Notes::

*1: The temperature-dependence factor. Please see the first note in Box 1.

*2: Gate-state dynamics parameters of the K + or channels. Actual channel opening/closing rates are given by Boltzmann-like functions—for details see the second note in Box 1.

*3: The 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.

Table 4

Definition and notation for the key HH variables.

Notation	Variable description and units	Typical value (*1)
Potentials, in mV:
V m or V	Membrane voltage	(*3)
V rest	Membrane resting voltage	-77
E K	K + Nernst potential	-90
E Na	Na + Nernst potential	60.0
E Leak	Leak reversal potential	-70
Membrane capacitance, in μF/cm 2:
C m or C	Membrane capacity
c	Membrane capacitance	1
Maximum (*2) conductances, in mS/cm 2:
g K	K + conductance	150
g Na	Na + conductance	300
g Leak	Leak conductance	0.033
Currents, in μA/cm 2:
I K	K + Ionic Current (*4)	g K × n × (V m − E K)
I Na	Na + Ionic Current	g Na × m 3 h × (V m − E Na)
I Leak	Leak Current	g Leak × (V m − E Leak)

Notes:

*1: Typical values are for the Na 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 ; is depolarized (grows due to stimulation and/or activated sodium Na + ion channels); is repolarized (decays back to V , 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.

Fig 4

State-space dynamics and model regimes as a function of ΔV 1/2.

These bifurcation diagrams illustrate the FP’s and PO’s of the V and h model states as a function of ΔV 1/2. Thin black line: stable FP’s; thick cyan: the unstable-FP branches (with one or two (real or complex) positive-real eigenvalues. LP1 and LP2: lower and upper limits of the middle branch. HB: Hopf bifurcation. See the text for more details. Thin and thick dashed lines: minimum and maximum model-state value for the stable and unstable cycles respectively. Unstable cycles are born from FP’s at a subcritical HB and metamorphose to stable at a cycle saddle-node (fold).

Box 1: Gate-state dynamics mathematical details

*1: Temperature dependence assumed linear with slope , where T 0 = 23°C.

*2: Opening/closing rate of gates y in the K + or channels have Boltzmann-like templates (⋅ subscripts dropped to simplify notation): where w = V − 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: α′(w), β′(w), and y ∞(w) are all sigmoidal functions of w. Regardless of a, b, k values, rate-change dependence on w is most important at w = 0, where α′(w) = a/2 and β′(w) = −b/2. Compare to α′(w) ≈ β′(w) ≈ 0 for w ≫ 0.

is reciprocally dependent on k, and its maximum is obtained at w* = −kln(a/b) (hence with a = b, w* = 0). Hence dy ∞(w*)/dw is determined by the V 1/2 parameter, and y ∞(w*) = 1/2—regardless of the actual a, b, k values. Also for V ≈ V 1/2, τ(V) = τ(w) attains its maximum (τ(w*) ≈ τ(0) = 1/k(a + b)), since: E.g. for a = b = k = 1, the numerator of τ′(w) reduces to u 2+2uw − 1 (u = e −), which vanishes when w = 0.

*3: Inactivating h gates of the channels have the asymptotic state

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 variation). In this model extension, a fraction P ≤ 1 of all the Na + ion channels are of the Na 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 ≠ 0. Thus: 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 multiple-compartment 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 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 = 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 + 1.5 ms. However, this also provided for a little extra variability in the identified I , 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: .

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 whole-neuron 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 ∞(V) while h and n remain at their resting values h and n . Upon very brief I duration this is a fair approximation of the system’s behavior (see Ch.3 in [37]). Stimulating from rest [V , h , n ], the system dynamics is reduced to: For all ΔV 1/2, the above approximate ionic current—denoted I (V)—a function of just V, has 3 zero-crossings: V , V (see Fig 1D) and V (not shown). The LOCAL minimum −I is reached for a V value located between V and V . While V and V remain virtually unchanged for different values of ΔV 1/2, V moves toward V and I diminishes as ΔV 1/2 is decreased. In this 1D approximation, one can see that the membrane will depolarize further toward V if at the end of the stimulation V > V . The latter is possible only if I > I . Hence lower ΔV 1/2 increases the excitability by reducing both V − V and I .

Strength-duration (SD) curves for the ΔV 1/2-family, computed through simulations of the full HH dynamics, and for the whole range of ΔV 1/2 ∈ [−8, 40] mV, are shown in Panel A of Fig 2. For the short (e.g. T = 100 μs) the increase of I with ΔV 1/2 is close-to-linear, while for the long (e.g. T = 2000 μs) the increase is steeper than linear.

Fig 2

Excitability as a function of ΔV 1/2.

A Strength/duration (SD) curves; B I (t − t ) T = T = 100 μs. ∀ΔV 1/2 the test AP is obtained by 1.5× the resting I (see Panel A).

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 (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 no longer exists, when the h and n gates respectively close and open. With three zero crossings −I (V) also has a local minimum—a situation qualitatively similar to the I (V) case above. Hence, just as with the simplified model of Eq (3), further depolarization of V toward V may be expected when V is beyond V at the end of the stimulation. In this case I is precisely the value of the rheobase current.

Through I (V) of Eq (4), a lower bound of the threshold current—the absolute rheobase—can be estimated.

Hence, SD[V ] may be estimated for relatively short durations (T → 0) or long durations (T → ∞). In the latter case, this is only possible for models like Na (ΔV 1/2 = 0), but not for models like Na . This is one of the key properties pointed at in this work.

To evoke an AP, according to Eq (1) the stimulation current I needs to overcome the dominating opposing ionic currents I . In search for patterns of lowest amplitude or highest energy efficiency, I may just slightly outweigh I and still successfully trigger AP’s [36]. However, this is highly dependent also on stimulation duration and Na class.

Importantly, for the Na subtype, but not for the Na 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 , induced by an initial supra-threshold stimulus S1. The minimal current I , 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 is a function of t = t − t coupling interval (assume t = 0 for simplicity), as well as of the S2 duration T . 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 .

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 does not exist, the system is absolutely refractory. Otherwise, a minimum stimulus current is needed.

We consider the case where AP is evoked by a stimulus applied from the resting state. Both V (t)—calculated from Eq (5), and I (t) depend on the AP , which in turn varies with T (and the corresponding I ). However, both functions become invariant once the time axis is aligned to t —the time of maximum dV/dt derivative during the upstroke. Different I and T may modify the latency of AP , but almost none of its time course beyond the upstroke. This invariance is maintained as long as T does not become too long. Conversely, when T → ∞, I 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 time course means that a single I (t − t ) curve (excitability is a function on the recovery time t − t ) can be used to capture the system behavior for each T 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 as a function of t − t and ΔV 1/2 for T = 0.1 ms. Decreasing ΔV 1/2 slightly shortens the absolute refractory period. For any given t − t , it also reduces the current needed to get a response. For any t − t , there is a quasi-linear increase of the I as a function of ΔV 1/2. For t − t > ≈ 2.5 ms, where excitability regained in all model Na 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 /I currents interplay and metabolic efficiency as a function of ΔV 1/2 are illustrated in Fig 3. ∀ΔV 1/2 the AP is obtained via 1.5× the resting I for duration T = 100 μs. As observed, larger ΔV 1/2 values translate into lower excitability: higher threshold potential V and threshold stimulation currents I . The very large latencies for ΔV 1/2 = 0 and -5 mV are due to very low threshold values I , which remain low even after a 50% “safety-factor” increase. The thick color traces in Panel A of the figure show the −I 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 /I interplay here gives an electrophysiological and physical perspective onto the nonlinear dynamics concept of the “slow manifold” [37, 42, 43].

Fig 3

AP properties and Metabolic efficiency as a function of ΔV 1/2.

A Membrane voltage during the AP; B I / I interplay; C Overall Metabolic efficiency Thick color: −I , aligned to the AP upstroke (time zero); Dashed black: the respective I Q —charge transferred into the modeled compartment during a single AP. AP shape and Q nearly invariant with I ≫ I .

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 re-polarization stage. First, lower ΔV 1/2 implies lower V . Second, the AP upstroke amplitude V is invariant to ΔV 1/2 and hence h ∞(V ) is significantly more closed for smaller ΔV 1/2 values (see Fig 1A).

The resting FP V remains almost identical in the Na and Na models (-77.03 vs -77.01 mV), but nevertheless corresponds to quite different values of h ∞(V ) (0.96 vs 0.76). Hence, the Na 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 model has two additional depolarized FP (diamond and square in Fig 1C—denoted V and V ).

Compare the equilibria of the gate state h on Fig 4 for the typical models of Na (ΔV 1/2 = 0), and Na (Δ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 (in micro-Coulomb’s per cm 2), which is transferred during a single AP ([41, 44], T = 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 are mostly invariant as long as I ≫ I . 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 is lower in this case, this earlier Na + inactivation is matched to a significantly lesser activation of the voltage-gated K + channels.

Bifurcation structure as a function of the ΔV 1/2 parameter

The resting potential (V ) 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 (V) + I (V) with the bell-shaped ΔV 1/2-dependent inward current I (Fig 1C). The effect of ΔV 1/2 variation is a linear shift of the window-current I [Δ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 has a single unstable direction. Perturbations along this direction (either positive or negative) produce two heteroclinic trajectories ending at V , 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 to the stable resting state. All these invariant trajectories and their associated stable manifold do not exist for models like that of the Na 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 5

Complex dynamic organization for ΔV 1/2 = 0 mV.

Two heteroclinic trajectories, starting at V (red diamond) and ending at V (green circle): correspond to passive—green and active depolarization—red. Trajectories from V (blue square)—thin dotted blue connect into V . All trajectories avoid and “slide” on the slow manifold (SM) (wireframe mesh), constructed by solving I (V, m ∞(V), h, n) = 0 ∀(n, h). Its ‘cubic-like’ shape defines the excitation threshold.

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 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 ∈ [1.7 mV, LP1]), or two positive-real complex eigenvalues (ΔV 1/2 ∈ [HB, 1.7 mV]).

Table 5

Bifurcation-terms Glossary [43, 49].

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.
Saddle loop or homoclinic	a given PO’s amplitude increases until it captures a saddle point and disappears; the PO’s period a →∞ as the parameter approaches the critical value
Double cycle or saddle-node of cycles (SNC)	Two POs appear or disappear coalescing onto a single half-stable PO, which in its turn disappears
Codimension 2 bifurcations
Cusp (C)	Three fixed points appear or disappear this is a sort of LP generalized by the presence of the second parameter
Takens-Bogdanov (TB)	An LP and HB coalesce into an FP yielding an ODE Jacobian with two zero eigenvalues; the HB-related PO becomes a homoclinic.
Generalized Hopf (GHB)	Two HB coalesce, as the second parameter reaches a critical value

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 Results—subsection 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:

ΔV 1/2 > LP1. The corresponding dynamic system has a single stable resting (hyperpolarized) FP, e.g. the Na channel subtype.

ΔV 1/2 ∈ [LP2, LP1] yields 1 stable resting FP, and two unstable FP’s, e.g. the Na channel subtype. As can be inferred by Fig 1C, increase or decrease of the nominal conductivity g 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 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 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 ∈ [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 —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 .

Panels A and C in Fig 6 result from exactly such bifurcation analysis for the typical models of Na (ΔV 1/2 = 0), and Na (Δ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].

Fig 6

I BD’s for the Na and Na ion channels (ΔV 1/2 = 13 and 0 mV, resp.).

Panel B in Fig 6 illustrates the period length of automatic firing as a function of I . In the Na 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 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 required for an 1:1 response—at the given minimal pacing period τ 11, depends on pacing frequency. Clearly, I → I (the resting threshold, see also Fig 2), as τ 11 → ∞. Notice that, with near-threshold current values, higher ΔV 1/2 allow periodic stimulation at higher frequency.

Fig 7

Limits of automatic firing and 1:1 pacing as a function of ΔV 1/2.

A Stimulation current (also function of frequency); B Min. I (→ Max. period length) for stable cycles I required for an 1:1 response at pacing period τ 11. T = 100 μs; I → I as τ 11 ≫ 0; Semi-transparent quasi-planar surface = 1.5× the resting I .

However, the situation is reversed for automatic regimes. The period of automatic firing is a monotonously decreasing function of I 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 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 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 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 were identified. Fig 8 illustrates this for two post hoc runs of the central step of the I -search algorithm (ΔV 1/2 = 2 mV, please see Methods for a detailed interpretation of the obtained voltage-distribution patterns in either case).

Fig 8

A generalization perspective of ΔV 1/2 variation effects—from a compartment to a cell.

A When I 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 .

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 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 with ΔV 1/2 for both T 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.

Fig 9

Resting threshold current I and mid-axon mean AP conduction velocity.

Notice the clear linear trend of I increasing with ΔV 1/2 for both T cases. A significant trend of accelerating AP propagation with ΔV 1/2 can be acknowledged. To prevent excessive results variability (see Fig 8 and text), conduction velocity was computed for I = 1.5 × I .

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 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 saddle-node bifurcation (see Table 5) for P ≈ 0.648, which is close to the minimum relative conductivity (g /g = 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 [ΔV 1/2 = 0](V), the contribution of I [ΔV 1/2](V) is so low that it does not change much the total I 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 [Δ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 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.

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 -subtype model with ΔV 1/2 = 0 mV. In the voltage range of peak M3 current, the contribution from the Na 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 channels. B: Resting threshold of the mixed models as a function of the P parameter. The figure shows κ(P) = I (P, ΔV 1/2)/I (1,ΔV 1/2) of the mixed model for T = 100 μs. Notice that there is an approximately twofold difference between the thresholds of the M3 and M1 single-Na -subtype 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 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.

The relative resting threshold κ(P) = I (P, ΔV 1/2)/I (1,ΔV 1/2) of the mixture model is a function of the P parameter (T = 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 section on neuronal excitability, both I and V − V provide measure of excitability. Both I (Fig 10B) and V − V (data not shown) have decreased by more that 70% at P = 0.2.

Discussion

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 ∞(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-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 est = 0), which is quite unlike the classical HHM (where n est ≈ 0.2). Hence, the resting FP is fully determined by the leak parameters (conductance and Nernst potential). This explains the invariance of V 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 and K 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 . Rinzel and Miller examined the I 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 × E [49, 50] and the two-dimensional I + I HH-type model [37] with its either high- or low-threshold I (i.e. two E levels). The latter model is equivalent in many respects to the well-known I + I 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 is manipulated toward changing the I 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 conductance, relative to nominal g . 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 or g 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 and Na 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 stimulation-current 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 ∈ [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 (V) has a local minimum −I between V and V (Fig 1C and 1D). Hence an AP can be evoked by applying as little current as the absolute rheobase current I for a sufficient amount of time T (Fig 2).

In general, lower ΔV 1/2 values translate into lower threshold current I [ΔV 1/2](T ) from the resting state (Fig 2).

The mere FP count corresponds to nontrivial differences in dynamic organization and its complexity. E.g. supernormality of I as a function of t − t , 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 parameter plane). Fig A in S1 File demonstrates transitions between pure excitability and oscillation determined by the I and ΔV 1/2 values.

The automatic-firing region of I 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 as a function t − t 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 ∀ΔV 1/2 a similar proportion of I relative to the resting threshold I [ΔV 1/2](T ) (see also Fig 2), Fig 7B predicts that the higher-ΔV 1/2 models would be the first to behave as oscillators as T 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).

The the semi-transparent surface in Fig 7A represents the (1D) dependence of resting threshold I on ΔV 1/2 for T = 100 μs, applying ∀ΔV 1/2 a 150% safety factor for robust AP-triggering—i.e. This is intersects into the main 3D surface, which shows the minimum 1:1 pacing I currents as a function of ΔV 1/2 and the pacing interval τ 11 (the black-yellow line in Fig 7A).

For I close to the I 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 frequency-encoding range, above).

However, when I is much above the I 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 ) 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 compared to I .

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 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 sub-populations. 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 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.

Appendix A: Supplementary Results

The effect of ΔV 1/2 on automaticity

Applying the S1/S2 protocol to study refractoriness implicitly assumes that the system does not reach an automatic regime—i.e. produce multiple AP’s during either of the S1 or S2 stimuli. For brief stimulation, the system reverts to the non-stimulated dynamics during re-polarization. However, as T increases, new attractors associated to long stimulation cannot be ignored—the system may now approach them. Hence, the outcome depends on the duration of the transient trajectories leading to these attractors. The latter are revealed by the bifurcation structure as a function of a constant bias current parameter (I ) Such current is known to change the attractors of the system (e.g. [37], Chapter 4).

Figs B and C in S1 File present examples of the rich bifurcation structure created by a constant stimulation current I for different values of ΔV 1/2. As ΔV 1/2 decreases, the FP curve changes from a monotonous increasing single-branch (Fig C in S1 File, top panel: ΔV 1/2 = 13 mV) to a three-branch organization (e.g ΔV 1/2 = 3 and -8 mV). The stability of the hyperpolarized FP is lost through a subcritical Hopf bifurcation (see Table 5). The most depolarized FP also becomes stable by a supercritical Hopf (HB) at a very high value of bias current. High amplitude stable PO’s also exist in an I range specific to each ΔV 1/2. These stable cycles are created by a saddle-node of cycles bifurcation (see Table 5), which produces also unstable cycles. The unstable PO’s may connect with the subcritical HB (of the lower branch), as in the case of ΔV 1/2 = 13 mV. The interaction of the unstable PO’s with the low-ΔV 1/2 (yielding a three-branch FP structure) is interesting. They may end through a homoclinic bifurcation (see Table 5), when hitting the middle FP branch, as in the case of ΔV 1/2 = 3 mV. In the latter case, the unstable cycles created at the HB1 also disappear through a similar homoclinic bifurcation, but at a different (higher) bias current. Hence, different situations are possible depending on the values of ΔV 1/2 and I : a unique stable FP; bi-stability between a stable FP and a stable PO, which in turn are separated by an unstable cycle and/or two unstable fixed points; a unique stable cycle with one or three unstable fixed points (see the zoomed-inset in Fig C in S1 File, the respective phase-plane trajectories, corresponding to the “menagerie” of possibilities, are shown in Figs H,I and J in S1 File).

How long does the transient take to reach a stable cycle from the resting state? Some clue is provided by the cycle period. As seen in Fig D.B in S1 File, the period is a monotonously decreasing function of I with a minimum well below 5 ms.

The codimension-2 bifurcation structure in the ΔV 1/2 × I parameter space

I stands for the magnitude a constant current injected into the modeled compartment and is known as bias current.

Already from the phase portrait of typical sub- and supra-threshold trajectories (Fig 5) one gets a pretty colorful idea about the dynamic-structure richness of the parametric HH-type model. Consistently to [49, 50], in this analysis we found a very rich a codimension-2 bifurcation structure.

Fig A in S1 File presents the full picture of the rather complex codimension-2 bifurcation structure in the (ΔV 1/2, I ) parameter-plane (see also Fig B in S1 File). The main features are:

Cyan trace shows the codimension-2 positions of the two Hopf bifurcations HB1 and HB2 as a function of the model parameters. Note that for the low ΔV 1/2 there is an only-apparent “collision” of the two (but see also Fig D.A in S1 File). While for very high ΔV 1/2 the two HB’s indeed coalesce, which in codimension-2 is known as Generalized Hopf (see Table 5). 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 range of stable/unstable PO’s determined by the saddle-nodes 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 Fig B.A in S1 File and the case of ΔV 1/2 = 3 mV in Fig C in S1 File). For very large I currents (beyond the blue trace) only the stable depolarized fixed point remains (apart when the island described next exists).

For very large ΔV 1/2 the stable/unstable PO’s between SNC1 and SNC2 detach from the FP’s locus to form an island (see Fig G in S1 File).

Supplementary Figures.

Fig A The Δ Bifurcation structure in the ΔV 1/2 × I parameter plane Cyan trace: the Hopf bifurcations (HB) as a function of the I and ΔV 1/2 model parameters (see also Fig D.A). Dashed-black lines: Existence/creation parameter-range of the 2 additional fixed points (saddle-node in the middle branch, and unstable center/focus in the top branch of the B.D.’s, see Fig 4A and 4C) Red and blue lines: the ΔV 1/2 × I range of stable/unstable PO’s (see also Fig G) Notice that outside (below and over) the HB area, there is bistability with a stable fixed point, while for very large I currents (beyond the blue trace) only the stable depolarized fixed point remains (see also Fig C). For very large ΔV 1/2 the stable/unstable PO’s between SNC1 and SNC2 detach from the FP’s locus to form an island (see also Fig F and G)

Fig B The Bifurcation Glossary (Codimension 1 and possibly 2) terms illustrated for two Δ black lines: stable FP’s. They lose (or recover) stability via Hopf Bifurcations red lines: unstable FP’s. In A, unstable FP’s appear (or disappear) via Saddle node (SN) bifurcations at LP1 and LP2. blue lines: Minimum and Maximum voltage of stable periodic orbits. These can appear (or disappear) through SNC bifurcations or supercritical HB’s (e.g. the low amplitude stable cycles in B). cyan lines: Minimum and Maximum voltage of unstable periodic orbits. These can appear (or disappear) through SNC bifurcations or subcritical HB (e.g. HB2 in A), or homoclinic bifurcation (e.g. the lower open ends in A).

Fig C BD’s for a set of Δ Inset: Zoom in of the ΔV 1/2 = 3 mV case for I ∈ [7, 14]μA/cm 2

Fig D Combined BD’s and periods of the high-amplitude cycle—in the Δ A: Combined BD’s B: max. cycle-periods as a function of ΔV 1/2 × I parameters (B) Note that the periods of the low-amplitude cycles (around SNC3 and SNC4) have not been continued.

Fig E BD’s for a set more of Δ (see also Fig C)

Fig F Generalized Hopf in codimension-2: Collision of two Hopf Bifurcations (HB) As the two Hopf’s become closer in parameter space, one unstable PO (UPO) is lost. Moreover, the HB type changes from sub-critical to super-critical (see also Fig G)

Fig G The creation of an A little before the 2 HB’s collide, the small-amplitude SPO detaches from the UPO. The latter forms an island together with the large-amplitude SPO. (see also Fig F)

Fig H Very complex dynamic organization for ΔV 1/2 = 3 mV and I = 9.5 μA/cm 2 A: Temporal evolution of membrane voltage V for the SPO (magenta, green trace on Panel C) and the UPO (black, thick cyan trace on Panel C). B: Gate states dynamics for the SPO (green trace on Panel C). Notice the very limited range of the n gate’s variation. C: The unstable FP invariant directions yield transient trajectories which converge mostly to the SPO. This is due to the UPO (thick cyan trace). The transient trajectories starting at the UPO’s single unstable invariant direction converge respectively to the resting FP and to the SPO (thinner dashed cyan traces). See also Fig H.A

Fig I Examples of dynamic objects (phase-space trajectories) corresponding to bifurcations in codimension-2. I variation for ΔV 1/2 = 3 mV results in a complex organization (see the B.D. in Fig H—middle row and Inset). The figure presents the rich diversity of dynamic objects (e.g. heteroclinic, homoclinic, stable and unstable periodic orbits) by simulating them for 3 values of the I parameter: 7.44, 7.49 and 9.068 μA/cm 2. Heteroclinic (dashed red), homoclinic (blue), stable (green) and unstable (cyan and black, Panel A only) periodic orbits

Fig J An unstable limit cycle as a phase-space separatrix A: Due to the existence of an UPO (cyan), all trajectories (black and red) initiated on the unstable FP’s (red triangle) invariant directions deviate toward the resting FP (blue circle) and not to the SPO (green). Incidently (and like in the B.D.’s on Fig F.C), the UPO’s own invariant direction yields 2 transient trajectories (data not shown) which converge respectively to the resting FP and to the SPO. B: An UPO is absent here unlike Panel A. Hence, the same initial conditions as in Panel A yield transient trajectories half of which (2 out of 4) converge to the SPO.

(PDF)

Click here for additional data file.