Spontaneous Excitation Patterns Computed for Axons with Injury-like Impairments of Sodium Channels and Na/K Pumps

In injured neurons, “leaky” voltage-gated sodium channels (Nav) underlie dysfunctional excitability that ranges from spontaneous subthreshold oscillations (STO), to ectopic (sometimes paroxysmal) excitation, to depolarizing block. In recombinant systems, mechanical injury to Nav1.6-rich membranes causes cytoplasmic Na+-loading and “Nav-CLS”, i.e., coupled left-(hyperpolarizing)-shift of Nav activation and availability. Metabolic injury of hippocampal neurons (epileptic discharge) results in comparable impairment: left-shifted activation and availability and hence left-shifted INa-window. A recent computation study revealed that CLS-based INa-window left-shift dissipates ion gradients and impairs excitability. Here, via dynamical analyses, we focus on sustained excitability patterns in mildly damaged nodes, in particular with more realistic Gaussian-distributed Nav-CLS to mimic “smeared” injury intensity. Since our interest is axons that might survive injury, pumps (sine qua non for live axons) are included. In some simulations, pump efficacy and system volumes are varied. Impacts of current noise inputs are also characterized. The diverse modes of spontaneous rhythmic activity evident in these scenarios are studied using bifurcation analysis. For “mild CLS injury”, a prominent feature is slow pump/leak-mediated EIon oscillations. These slow oscillations yield dynamic firing thresholds that underlie complex voltage STO and bursting behaviors. Thus, Nav-CLS, a biophysically justified mode of injury, in parallel with functioning pumps, robustly engenders an emergent slow process that triggers a plethora of pathological excitability patterns. This minimalist “device” could have physiological analogs. At first nodes of Ranvier and at nociceptors, e.g., localized lipid-tuning that modulated Nav midpoints could produce Nav-CLS, as could co-expression of appropriately differing Nav isoforms.

In hippocampal neurons, a Nav-CLS based window current left-shift occurs following metabolic insult from prolonged epilepsy-like stimulation [15]. Hippocampal neurons exposed to the industrial compound, melamine, also show Nav-CLS and exhibit hyperexcitability [16]. In neurons and muscle fibers, leftshift of Nav availability (i.e., steady-state inactivation) [17] is reported for many damaging and/or toxic conditions (see Table 1 of [3]) and for membrane fluidizing conditions [18]; fast activation also becomes left-shifted in muscles (e.g. [19]) and neurons (e.g. [20]). As predicted when membrane injury causes Nav-CLS [11], sick excitable cell pathologies range from hyperexcitability through ectopic excitation to depolarizing block. Sucrose-gap voltage clamp of nodes of Ranvier subjected to ischemic and other insults could, we have suggested [3], be used to test to what extent the initial Nav-leak of injured Nav-rich native membranes is explained by Nav-CLS. Injured-bilayer I CLS-window could help explain [3] why, invariably, clinical Nav inhibitors protective against excitotoxicity are lipophilic [21]; moreover, if bilayer structure matters, then retuning therapeutic strategies to better target channels in damaged bilayer could pay dividends [3].
Here, our goal is to understand in dynamical terms the mechanisms of abnormal steady-state excitability in membrane modeled with CLS-type injuries of various intensities. In particular, we wish to determine if -in conjunction with axolemmal Na/K pumps -CLS-injury predicts forms of abnormal excitability reported for ''mildly'' injured excitable cells, i.e., cells that, with appropriate remediation, could be salvaged. In healthy unstimulated axons, quiescence is the norm and spontaneous firing is neuropathic and ectopic.
Several recent studies involving simulation and dynamical modeling addressed abnormal or paroxysmal discharge in injured neurons. Injured DRG neuron simulations [22] showed that subthreshold oscillations (STOs) resembling ones recorded in vivo and relevant to neuropathic pain can be generated via depolarizing stimuli to a membrane with three distinctive g Na s. Importantly, though, for STOs, the slower two g Na s required first order, rather than the usual third order activation kinetics. We address this later. Similarly, for two Nav-type and two Kv-type conductances, Choi and Waxman [23] showed responses to current injection (mimicking sensory input) that included STOs with the spikes; there too, for STOs to occur, one of the g Na s required firstorder activation kinetics.
STOs typically occur in conjunction with bursts of spikes, appearing in healthy neurons, (e.g., entorhinal cortex [24]) as Figure 1. Schematic of a mechanically-injured node of Ranvier depicted with a mix of intact-looking and severely-blebbed axolemma (as labelled) such as seen in transmission electromicrographs of stretch-injured optic nerve nodes [2]. In pipette aspiration bleb injury, the cortical actomyosin-spectrin skeleton progressively detaches [6]. Our model considers a node as one equipotential compartment in which actual spatial arrangements of pumps and channels are irrelevant. However, the fraction of Nav channels in the injured portion of the membrane, along with the severity of their gating abnormality, are model parameters. doi:10.1371/journal.pcbi.1002664.g001

Author Summary
Nerve cells damaged by trauma, stroke, epilepsy, inflammatory conditions etc, have chronically leaky sodium channels that eventually kill. The usual job of sodium channels is to make brief voltage signals -action potentials-for long distance propagation. After sodium channels open to generate action potentials, sodium pumps work harder to re-establish the intracellular/extracellular sodium imbalance that is, literally, the neuron's battery for firing action potentials. Wherever tissue damage renders membranes overly fluid, we hypothesize, sodium channels become chronically leaky. Our experimental findings justify this. In fluidized membranes, sodium channel voltage sensors respond too easily, letting channels spend too much time open. Channels leak, pumps respond. By mathematical modeling, we show that in damaged channel-rich membranes the continual pump/leak counterplay would trigger the kinds of bizarre intermittent action potential bursts typical of injured neurons. Arising ectopically from injury regions, such neuropathic firing is unrelated to events in the external world. Drugs that can silence these deleterious electrical barrages without blocking healthy action potentials are needed. If fluidized membranes house the problematic leaky sodium channels, then drug side effects could be diminished by using drugs that accumulate most avidly into fluidized membranes, and that bind their targets with highest affinity there.
small V m oscillations for a discrete range of conditions. Dynamically-speaking, STOs emerge where there is bistability between tonic firing and quiescence [25]. For allodynia-related neuropathic pain, a dynamical study of firing patterns showed E Anion changes having a key stability role [26]. Firing abnormalities of demyelinated axon were mimicked in a dynamical modeling study with fixed Nernst potentials and g Na /g leak ratios varied [27]. Regarding paroxysmal spiking, the same group [28] explored initiation and termination of afterdischarges and spontaneous activity; there, a non-inactivating g Na induced bistability and to mimic [ion] homeostasis following transient Na + -loading, an exponential decay provided a slow relaxation toward initial [Na + ] intracellular . Kä ger et al [29] investigated epilepsy-like spontaneous bursting patterns and elevated [K + ] extracellular , and Barreto, Cressman and colleagues [30,31] performed a dynamical study for oscillating [ion]s connected with epileptic seizures.
Our Nav-CLS-based model of membrane injury uses a Hodgkin-Huxley axon without or with a Na/K pump as developed previously in [11]. That model focused on transitory responses to abrupt CLS ''injury'' (left-shift intensity and the fraction of total g Na affected were varied) and revealed how CLSinjury affects V m (t) and ion homeostasis. Once the injured system (no pumps) stabilized, depolarizing current was applied, and based on spontaneous and stimulated behaviors, Nav-CLS injury activity was classified under various regimes ( Figures 5 and 9 in [11]). CLS-injury caused hyperexcitability, ectopic excitation or depolarizing block, and interfered with saltatory propagation of normal action potential (AP) traffic. Varying extracellular volume can significantly alter the quantitative measures of the bursting, consistent with a major role for pump activity. Unexpectedly, it was discovered that Nav-CLS in conjunction with pumps could generate STOs. Since injured neuron neuropathic pain discharge is characterized by diverse STO phenomena, we undertook the dynamical analysis of Nav-CLS injury presented here.
We model acute injury: even in egregiously blebbed native membranes, Nav channels retain their ability to respond to voltage [32,33]. Keeping g Na-max constant corresponds to an absence of de novo Nav channel expression or membrane remodelling. Immunochemistry confirms an abundance of axolemmal Nav channels at damaged nerve-end neuromas, sites considered to be loci of ectopic excitation [34]. Our CLS-injury model is consistent with present views of voltage-gated channel structure/function [35][36][37]. Nav-CLS injury in situ would likely be of non-uniform intensity, a point we address computationally in several ways. Except for some instances of modified pump intensity, ''impairment'' refers exclusively to Nav-CLS.
Na/K pump activity is sensitive to ATP levels, injury and bilayer lipid characteristics [38][39][40], but how mild membrane injury would affect pump rates is unclear. So here we vary the maximal pump rate to explore its effect on tonic firing and bursting (in this paper, both firing patterns are spontaneous).
The diverse spontaneous steady-state rhythmic activity patterns that became evident in our injury scenarios are studied using bifurcation analyses. For mild CLS injury, a prominent and robust feature was the emergent slow process of pump/leak-mediated E Ion oscillations [11]. As we will illustrate, these slowly oscillating E Ion values cause time-dependent changes in the firing thresholds, which in turn will lead to complex spontaneous STO and bursting behaviors; the complexity of the observed patterns is further exacerbated by current noise.
In summary, we study, via dynamical analysis, a simple, robust, and biophysically explicable model of mild axonal injury in which pumps remain functional, i.e., we examine a collection of mild Nav-CLS injury scenarios. The major finding is that interactions between the window current based Nav-leak and Na/K pump activity engender slow oscillations of the Na + and K + driving forces that spontaneously trigger a plethora of known neuropathological excitability patterns. The following results are based on simulations of a Nav-CLS model of injured node described in the Methods section and in Table 1 which is a list of parameters.

Membrane potential oscillations in CLS injured nodes for fixed E ion
Simulated nodes with Nav-CLS injury were shown [11] to exhibit activity regimes that include tonic and burst firing, plus a quiescent steady state. In this section we explore the node with simulated more nuanced versions of Nav-CLS injury and find that rhythmic firing combined with small-amplitude V m oscillations are commonplace. We start with cases of constant [Ion] (pump turned off), approximating axonal situations over brief times. We choose cases that provide snapshots of dynamics at selected fixed reversal potentials that could occur at different stages of injury. For fixed reversal potentials (E Ion ) we take three sets of values corresponding to different transmembrane Na + and K + gradients: (i) E Na = 50 mV, E K = 277 mV; (ii) E Na = 42 mV, E K = 277 mV when initial injury has elicited a change in I Na thence [Na + ]s; and (iii) E Na = 42 mV, E K = 271 mV. In this latter case, I Kpump , because of its dependence on [Na + ]s, changes, which results in E K changes.
Nav-CLS intensities induced by membrane damage would vary within and between axons, as suggested in Figure 8Aiii in [8]. To reflect this, three sets of LS i and f i are applied. To incorporate the idea of ''smeared-out LS'', the fraction f i of the population are given a Gaussian distribution of LS with mean6SD of 1.360.4 mV, 861 mV and 1561 mV. The numerical results are plotted in Figure 2. For the first two E Ion pairs (i) and (ii) (left and middle columns), V i dynamics in the injured node vary from quiescence to repetitive spiking, with the firing rate increasing and the spike amplitude decreasing as mean(LS) increases. However, for E Ion pair (iii) (right column), this node shows changes from STF at high amplitude (again, with increasing firing rate and decreasing  Figure 2. Nav-CLS induced spontaneous activities of injured node: steady state with fixed E Ion , tonic spiking, tonic subthreshold oscillations (STOs). By setting I maxpump , g Naleak , and g Kleak to zero and eliminating Eqs (11) and (12), the E Ion are artificially maintained at fixed values. resting potential (RP)), unstable fixed point, stable limit cycle (SLC, i.e., tonically firing APs), unstable limit cycle (ULC) and Hopf bifurcation point (HB). When E Na alone changes (J to K), the bifurcation structure shows only slight changes in the amplitude of SLC and the locations of subcritical HB on both ends. When E K changes (K to L), the HBs on both ends shift to the left (i.e. towards relatively smaller E Na ) and the previously subcritical HB (right side) becomes supercritical. Across the 3 columns, note that if the system did have pumps, interactions between the E Ion and I pump would continually and slowly change E Ion thereby repeatedly evoking these activity patterns. doi:10.1371/journal.pcbi.1002664.g002 amplitude) to a very low amplitude oscillation ( Figure 2I) that shares features of subthreshold oscillations (STOs).
The question of what happens to the threshold for STF for these changes in E Ion is raised by these results. The bottom row, panels J-K-L, displays the bifurcation diagram of the voltage excursions as a function of left-shift; the E Ion values for each panel are those used in the three panels above it. For these diagrams, the left shift is the bifurcation parameter; all channels (f = 1) assume this shift (there is no spread, i.e. SD (standard deviation) = 0). Each diagram is constructed by choosing a value of LS, running a simulation, letting transient behavior die out, then measuring various characteristics of the steady-state voltage solution. In particular, here we measure the value of a fixed point (resting potential, RP) or the minimum and maximum values of a V m oscillation. This procedure is repeated over a range of LS values, and the characteristic values are plotted against LS.
As LS increases, we see in panels J, K, and L a Hopf bifurcation (HB) i.e. a transition between a stable resting potential (RP), or in the diagram, stable fixed point, and tonic firing. Further, there is a narrow range of LS where both RP and tonic firing coexist; this is known as a ''subcritical'' HB. Subcritical HBs are seen on the left part of each diagram, and on the right part of diagrams J and K. In panel L however, the right part displays no such coexistence, and the HB is ''supercritical''. This means that as LS increases further, there is an abrupt (although continuous) transition from tonic firing to RP. Upon decreasing E Na (panels J to K), the threshold moves to higher LS values, and the spike height is reduced. However, when E K is reduced (panel L), the threshold moves to a lower LS, which explains the onset of firing in panel C. The decrease of the repolarizing current also continuously raises the minima and lowers the maxima of the voltage excursions. This underlies the trend to low amplitude spiking behavior (cf. panels CRFRI). Thus, in axons where injury has caused Nav-CLS and is fostering gradual rundown (depletion) of the ion gradients, a range of abnormal excitability patterns could be expected to appear then disappear.
Note that the mean LS in Figure 2I corresponds to the limit cycle regime ( Figure 2L) -although it is very close to the fixed point (quiescent) regime which is to the right of the supercritical Hopf bifurcation in Figure 2L. This is why the spike amplitude is so small -but technically they are still spikes, not STOs; this would be the case even if the Hopf bifurcation were subcritical (as in the right part of Figures 2J and K). There is no injected current here, and the system can generate tonic firing of low amplitude spikes, so we cannot speak of depolarization block. However, depolarization block can be produced by injecting a current of sufficient magnitude (not shown).
Since the E Ion are controlled by intra-and extracellular [Ion]s (Eq (13)) which in turn depend on pump currents via Eq (10), these results imply: a) I pump will substantially affect excitability at injured nodes [41], and b) in the presence of a small rundown of the K + gradient, Nav-CLS is sufficient to produce STO-like voltage excursions. We will see in the two sections on bursting below (Bursting in injured nodes and Bifurcation analysis of bursting) that these excursions and STOs share a common origin.

V m oscillations in a CLS injured node with the pumps on
Whereas the transition from quiescence to tonic firing when E K moves closer to zero ( Figure 2B to 2C) is expected, since the move causes a small depolarization, the appearance of STO-like excursions ( Figure 2H to 2I) is less intuitively obvious. To better understand how pump current affects the emergence of small oscillations, we simplify the above multiple Gaussian-distributed g Na populations by forming two subpopulations: LS i = [0,15]mV and f i = [0.5,0.5], signifying that half of Nav channels have no shift, and the rest have a 15 mV CLS. This choice for g Na populations is taken from the tonic firing regime ( Figure 9B in [11]) that occupies the largest part (,40%) portion of the LS/f plane depicting five spontaneous activity regimes of injured nodes. We turn on the Na/K pump and mimic injury to pump function by varying the parameter I maxpump in Eq (10). I maxpump is a functional maximum, a rate that, as per Eq 13 would only be attained for [K + ] o = 0 and [Na + ] i = 0 (unrealistic in living systems). It differs among axons [42] and will vary with the cellular surface area, Na/ K-ATPase density, bilayer lipid composition [39,42] and fluidity [43], and temperature [44] among other factor [45]. Importantly, I maxpump can be reversibly decreased experimentally by specific inhibitors (e.g. strophanthidin, ouabain). Figure 3 demonstrates that when I maxpump falls from 95 to 30 mA/cm 2 , repetitive AP activity (dotted line) is replaced by tonic low-amplitude V m oscillations (solid line). While the V m -oscillation amplitude increases as I maxpump increases ( Figure 3D), the firing rate decreases. This suggests that, if I pump were to decrease, repetitive APs could change to STOs. Inhibiting the Na/K pump can further lead to a complete cessation of firing. This is in general agreement with experimental results [42]. It is noteworthy that because [Na + ] i and [K + ] o (Eq (10)) limit maximal pump current, the upper values of I pump are considerably less than I maxpump . In Figure 3B, e.g., with I maxpump at 95 mA/cm 2 , I pump fluctuates around 22.760.5 mA/cm 2 . The section Effects of varying maximal pump current discusses the impact of I maxpump on the bursting regime.

Bursting in injured nodes
After repetitive firing activity, the electrogenic Na/K pump mediates the restoration of the ion gradients on a relatively slower time scale. Figures 2E and 2F suggest that pump current, by driving E Ion changes, could also induce patterns of bursting in which volleys of high-amplitude spikes are separated by periods of low-amplitude STOs.
In this section we pursue this idea by simulating nodes with three g Na populations (LS i = [26.5,2,0]mV and f i = [0.2,0.08,0.72]) and I maxpump set at 95 mA/cm 2 . In Figure 4A a node injured this way produces bursts of spikes separated by STOs. The silent and active phases of such bursts are associated with pump current dynamics and the left-shifted window conductance (m 3 h). In the silent phase where V m is near the RP (260 mV), the Nav-CLS channels are a ''leaky'' Nav component ( Figure 4B and 4C) flowing into the injured node, but without depolarizing it above firing threshold, because outward I Napump slightly exceeds the inward sodium current (the sum of I Na and I Naleak ). The result: V m is destabilized and STOs develop. I Napump keeps decreasing (Figure 4Aii), becoming less than inward sodium current at about the time marked (pink star) in Figure 4Ai. Therefore Na + accumulates inside and K + outside the node, resulting in E Na and E K changes, based on Eq (13).
When changes in E Ion further reduce the firing threshold (system dynamics will be explained in the following section) the injured node would cross the diminished threshold and produce APs. With the consequent [Na + ] i increase, I Napump starts increasing. It eventually prevents the node from crossing the threshold and APs cease (green dot in Figure 4Aii) and V m falls back to the STO state around the RP. Repeated switching between clusters of spikes and STO intervals constitutes a form of bursting behavior that, given the involvement of the pump, should be sensitive to the surface to volume ratio, as per the numerical results of Kä ger et al [29]. With intracellular and extracellular nodal volumes tripled (Vol i = Vol o = 3610 29 m 3 ) and the surface area fixed (6610 28 m 3 ) ( Figure 4D), bursting period (BP -from the beginning of one burst to the beginning of the next) and burst duration (BD) (Figure 4Ai) are prolonged. This confirms that, for the Nav-CLS model of injury, both the burst parameters are indeed sensitive to the surface-to-volume ratio.
Interestingly, the solutions in Figures 4A and 4D do not exactly repeat. The number of spikes during the active phase of the burst varies slightly among bursts. Additionally, missed spikes at the end of the burst phase relate to the existence of period-doubled solutions (PD) as we will see in Figure 5C. Even though below we characterize them by an approximate bursting ''period'', BP, these solutions are likely chaotic.

Bifurcation analysis of bursting
As seen above, bursting is followed by quiescence and STO's are observed at the onset and offset of bursting. We now expose the dynamical structures underlying these behaviors. During firing, the ion gradients are slowly depleted due to the window current in the injured node. The pumps cannot compensate for the loss, and E Na and E K move towards zero. During quiescence the pumps replenish the concentration gradients, and the absolute values of E Na and E K increase. Spontaneous rhythmic bursts in Figure 4A are shown three-dimensionally in Figure 5A with dynamic (i.e. time-varying) E Ion . A bifurcation analysis based on a decomposition of the full dynamics into slowly and rapidly changing variables (e.g. E Ion and V m respectively, Figure 5B and 5C) helps identify the burst mechanism. This is known as a slowfast decomposition ( [24] and references therein). We first focus on the initial phase of the burst (e.g. pink star, Figure 4A). For now, E K and I pump are regarded as fixed, i.e. as constants taking their (arbitrarily chosen) values at the pink star (283.17 mV,14.96 mA/ cm 2 ). We pick a value of E Na , run a simulation, wait until transients die out, and measure the minimum and maximum values of V m . This procedure is repeated for the range of E Na values on the abscissa of Figure 5B. On the ordinate we plot the corresponding minima and maxima of V m . A striking feature is that in certain regions of the diagram there are more than two ordinate values; this is a hallmark of bistability as discussed further below. Since E Na = 40.2 mV at the pink star, the diagram demonstrates that the CLS injury places the system in a superthreshold oscillation (i.e. spiking) regime. The spiking trajectory at the pink star in Figure 4A and 5A corresponds to the stable limit cycle (i.e., stable periodic oscillation) in Figure 5B.
At the onset of firing (pink star in Figure 5A), the pumps are still driving E Na upwards and past (i.e. to the right of) the Hopf bifurcation point (see Figure 5B). In that region, the fixed point, around which the STO oscillates, has become unstable; the trajectory moves towards the only possible stable state, a stable limit cycle corresponding to repetitive firing (APs).
During firing E Na is decreasing and consequently this limit cycle shifts to the left and approaches the bistability regime. During this process, both amplitude and frequency of the limit cycle decrease. In the bistability region of Figure 5B, the system has two stable states: a steady fixed point or resting potential (i.e. RP, marked by solid line) and a superthreshold oscillation (i.e. APs, or tonic firing). Between these states is an unstable limit cycle. The system phase space is then split into two sets of points. In the first set, typically near the RP, phase space solutions converge to the RP. In the other, typically further from the RP, solutions converge to the tonically firing solution. In both these sets, initial conditions close to the unstable limit cycle yield solutions that diverge away from it in an oscillatory manner: in the first set the oscillation grows into the tonic solution, and below, oscillations damp out towards the RP. Further, a perturbation from the RP must ''jump over'' this unstable limit cycle in order to lead to tonic firing, and vice-versa.
Near the offset of firing (green dot), with E Na still decreasing, the dynamics move toward a region of period doubling (PD). Sometimes the trajectory falls off the limit cycle on the upper branch before reaching this region such as in Figure 5C. Other times during the same solution the trajectory enters briefly into the period doubling regime. This manifests itself as a longer interspike interval at the end of the firing phase (see Figure 4A, second burst). For parameter choices corresponding to bursting the trajectories are not periodic, i.e. they appear to be chaotic unless the transients exceed our simulation times. The precise dynamics in this region warrant more extensive investigation beyond the scope of this paper. Note also that the AP amplitude decreases together with E Na . This is expected as the ion concentration gradients drive the action potential. This suggests that the upper branch of the PD regime is not reached by the trajectories near firing offset.
Let us inspect Figure 5 in more detail. Figure 5A shows that both E Na and E K are dynamic, changing slowly during a bursting period. Figures 5B and 5C capture the behavior of the system near the two transitions (pink star and green dot) in Figure 5A. The pink bistability region in Figure 5B -corresponding to the pink star -sits between a limit point (LP: where a ''saddle-node of limit cycles'' bifurcation occurs), and a Hopf bifurcation point (HB: where a transition between a stable fixed point and a stable limit cycle occurs). The green region in Figure 5C lies between (at left) an area (PD) of complex dynamics where many period-doubling bifurcations occur and (at right) a HB. The unstable limit cycle solution (marked by open dots) between the stable fixed point and stable limit cycle is the aforementioned threshold boundary. The arrow in Figure 5B shows the direction in which E Na changes during the firing phase (E K changes as well). Likewise, the blue arrow in Figure 5C shows how E Na changes during the STO phase. During the bursting period, the bifurcation diagram in Figure 5B slowly morphs into the Figure 5C diagram and back. These two diagrams, snapshots of the system's possible steady state behaviors at the two ends of the bursting cycle, help us understand how the whole solution plays out during an entire bursting cycle.
Interestingly, from Figure 5B one might think spiking would continue until the point LP is reached, where E Na is almost nil. However, the pump's dependence on [Ion]s limits this decay, such that a PD region appears and truncates the spiking phase at a relatively large value of E Na (note the difference in E Na scales in Figure 5B and 5C).
The two-parameter phase diagram in Figure 5D provides a view of the burst dynamics from the E Ion perspective. Here the pink star and green dot have the same meaning as in Figure 5A-C (see caption). The black loop connecting pink star to green dot is the projection of the 3D trajectory in Figure 5A onto the E Na 2E K plane, and displays the E Na and E K trajectories during bursting. The solid and dashed curves are, respectively, sets of all the LP and HB points for the range of E Ion of interest (exemplars of these points were seen in Figure 5B,C). The area bounded by them is the bistability regime. This involves computation of many bifurcation diagrams along the entire bursting trajectory, choosing from the trajectory the two extremal pump values (pink for the lowest and green for the highest). One can imagine that, when the system state moves from pink star to green dot, the periodically varying I pump pushes the bistability regime between the pink-curvebounded and green-curve-bounded areas. Thus, injured node dynamics switch between different modes, resulting in bursting behavior.
Nav-CLS phase diagram: 3 g Na subpopulations Various spontaneous activities of an injured node with three g Na subpopulations (LS i = [26.5,2,0]mV) comprising different fractions are presented in Figure 6A. The fractions (f 1 and f 2 ) corresponding to the Nav-CLS populations form the xand yaxes; the fraction of intact channels is simply f 3 = 12f 1 2f 2 . When the percentage of channels with deeper injury (i.e. LS i = 26.5 mV) increases, the state of the injured node goes from RPRburstingRbistability (APs or RP, depending on the initial condition)R tonic firingRbistabil-ityRRP. For two random points from the burst and bistability regimes, V m is plotted in Figures 6B and 6C, respectively. Overall, STOs can arise from various combinations involving decay to or growth away from RP, but the particular behavior seen depends sensitively on the percentage and degree of injury.

Nav-CLS and susceptibility to noise-induced bursts
Neuronal noise arises naturally from a wide variety of sources including channel noise, synaptic noise, electrogenic ion pumps, and thermal noise. In this subsection we investigate the influence of current noise on spontaneous electric activity of a node with Nav-CLS injury by including a simplified Gaussian white noise sj on the right side of Eq 1 as an additive term. This additive noise j has zero mean, and its autocorrelation function is ,j(t)j(s). = d(t2s) where d(?) is the Dirac delta function; s is the noise strength (in mA/cm 2 ). For a node, the white noise can be taken as describing conductance and thermal noise. By definition it has equal power at all time scales, in particular beyond the fastest time scales present in our deterministic (i.e. noiseless) node model. Furthermore, because the central limit theorem applies to the sum of different independent fluctuation sources acting on V m , its amplitude distribution is assumed to be Gaussian (see for instance, [46]). For a deterministic bursting activity with burst period BP = 10.2 s illustrated in Figure 6B, additive noise dramatically decreases BP ( Figure 7A). It also makes BP irregular because the crossing of the threshold by the V m trajectory is now a probabilistic process, especially in the vicinity of the threshold [47].
When the node operates in a bistable regime, as in Figure 6C, this additive noise causes stochastic switching between the quiescent and tonic firing modes ( Figure 7B); noise-induced bursts appear. In fact, the injured node produces bursting in both the burst and bistability regimes of Figure 6A. With higher noise intensity, the mean BPs of noise-induced bursts is shorter. Such bursting is also referred to as noise-induced mixed-mode oscillations [48]. This phenomenon has been found in experimental settings, for example the bistable squid giant axon [49], electroreceptors of paddlefish [50] and entorhinal cortex neurons [51], as well as in excitable neuron models such as the Hodgkin-Huxley model [49], the FitzHugh-Nagumo model [48], and an excitable spine model [52].
The expanded sections in the right column of Figure 7 demonstrate that, upon introduction of noise, the silent phase of the burst behaviour in Figure 6B and the steady state in Figure 6C turn into STOs. This is observed in injured primary sensory neurons [53], and constitutes coherence resonance (CR) (or autonomous stochastic resonance) wherein noise alone excites subthreshold oscillatory responses in excitable systems (see [47,54] and references therein). Noise continually modifies the firing pattern because it randomly pushes the system into and out of different regions of phase space, -similar behavior has been described for thermoreceptors [46]. In terms of Figure 5C, when the system is on the lower STO branch, the noise induces the behavior seen at the nearby Hopf bifurcation. In summary, for our Nav-CLS system, noise can elicit a random mixture of the behaviors seen in the phase diagram in Figure 6.

Nav-CLS variance and diverse patterns of ectopic excitation
Next, we explored the spontaneous activity of an injured node with different Nav-CLS variances. Gaussian distributions of Nav-CLS were imposed with different means and standard deviations (SDs). As illustrated in Figure 8A Figure 8K). The corresponding mean V m excursions decrease with increased mean(LS) (Figure 8L). With a smear of Nav-CLS that encompasses a wider range of LS values (SD = 10 mV: Figure 8J; see Figure 8A)), the injured node fires tonically, as in the lower SD case for the same mean(LS) ( Figure 8G). Evidently the larger LS components of the SD = 10 mV distribution do not sufficiently dominate to yield the type of quiescence associated with a (smaller SD) 30 mV LS ( Figure 8I).

Effects of varying maximal pump current
In Figure 3D we saw that decreased I maxpump suppresses spontaneous tonic firing with CLS injury. In Figure 9 we look at the effect of varying the pump activity on bursting patterns in a mildly injured node (LS = 2.560.5 mV). I maxpump = 95 mA/cm 2 (the value used in Figure 9C) is our standard or healthy cell reference. I maxpump depends on the density of pump expression in the membrane, which will differ among axons and also for different regions of a given axon. It also varies with a multitude of factors, including injury conditions that affect the turnover rate at individual pump molecules. To explore these possibilities in a mildly injured node, we varied I maxpump over a wide range below and above the standard value in Figure 9C. Figure 9B shows that a lower I maxpump is associated with longer interburst intervals (IBI: see Figure 4) than for the standard case ( Figure 9C). Figure 9A, an extreme example in which excitability is lost, has an infinite IBI. Increasing I maxpump , by contrast, initiates burst of spikes and shortens both burst duration (BD) and burst period (BP) ( Figure 9D,E), reflecting the hyperpolarization produced by the electrogenic Na/K pump. The differences between Figure 9B and 9C are consistent with experimental observations in rat spinal interneurons where, during the course of pump blockade, the burst rate decreases (larger BP) [42]. It is noteworthy that for the uninjured system, varying maximal pump current had little impact on the steady-state V m (see Table 2). What vary substantially are the steady state E Ion . This is consistent with the mechanistic interpretation that I pump oscillations responding to chronic Nav-CLS-based Nav-leak underlie the array of bursting patterns in the many simulations above. The basic phenomena (ectopic APs and STOs) are extremely robust properties of nodal membranes with mild Nav-CLS and a pump, and appear in various manifestations without the need to invoke novel conductances or Nav channels with peculiar kinetics.

Nav-CLS, pseudo-first order activation and subthreshold persistent current
Unlike nodes, somata express multiple Nav isoforms, and in recent modeling work of injured somata, Kovalsky et al [21] showed that a system with three distinctive g Na sub-types (but no Kv and no Na/K pumps) yields STOs, provided the two slow g Na sub-types are given 1st rather than 3rd order activation. This provided the system with the ''several delayed components'' needed to mimic the repertoire of STO-triggered phenomena of live neurons. To simulate Nav1.7 and Nav1.8 interacting in a fashion that elicits STOs, Choi and Waxman [23] also depicted one of the g Na sub-types with 1st order activation. Interestingly, in our system, left-shifting a fraction of the m 3 h-based I window acts like introducing a lower order activation process (see Figure 10). What physical change in injured axons that could yield such pseudo-first order m ' (V) curves for total g Na , even though individual channels continue to exhibit their normal m 3 activation? Heterogeneity in the extent of bilayer damage could cause this by fluidizing a fraction of the nodal bilayer. This would yield a subthreshold persistent I Na . Differences in bilayer packing in different parts of the membrane [3,36] also arise developmentally in some circumstances. In a given region of plasma membrane, a membrane remodeling process that places some Nav channels in (say) cholesterol-poor [55] and others in cholesterol-rich subdomains could, in principle, create or ''tune'' physiologically useful subthreshold persistent currents.

Discussion
Here we demonstrated that different patterns of spontaneous activity can be expected from mildly damaged nodes of Ranvier in which some or all the Hodgkin-Huxley type Nav channels have undergone a coupled left-shift in their activation and inactivation kinetics. This is called Nav-CLS [11] and as a form of Nav-leak, could result in the ectopic activity that emanates from injured nodes. As detailed in the Introduction, there is growing evidence that activation and inactivation are left-shifted in various chemical and mechanical injuries in different cell types and Nav isoforms [10,[12][13][14][15][16][17][18][19]. In more focused studies, it has been demonstrated that these left-shifts are coupled [8][9][10], and this CLS in Nav1.6 (the nodal isoform) is the basis of our study [11]. In diverse sick excitable cells, it has been hypothesized that Nav-CLS occurs wherever channel-embedding bilayer becomes structurally degraded and fluidized as a result of mechanical or chemical insults [3,4].
The following are the main findings of this study which is, essentially, a dynamical analysis of a standard node of Ranvier model with Na/K pumps, under conditions of coupled-left-shift (CLS) injury to some or all of the Nav channels and injury to the Na/K pump in terms of lowered output. The analysis uses simulations and bifurcation theory.
1. With E Ion values fixed (i.e. pump activity irrelevant) at standard or mildly depleted values (10-15%), LS injury levels up to 20 mV (including for ''smeared'', i.e., Gaussian distributions of CLS) elicit either quiescence or tonic firing. With increasing LS, a transition from quiescence to tonic firing occurs via a subcritical Hopf bifurcation, then at yet higher LS, tonic firing is replaced by quiescence via another subcritical Hopf bifurcation. If E Na and E K both approach zero, that bifurcation is supercritical; in its vicinity, tonic APs of arbitrarily small amplitude occur. 2. Maximal pump current is an important determinant of nodal excitability. When increased above a certain value, it yields an abrupt transition from quiescence to a relatively high rate of tonic firing; this is a supercritical Hopf bifurcation. 3. There exist CLS injury conditions where the system is bistable, exhibiting either quiescence or tonic firing. In this bistability regime, behavior depends on initial conditions in a complicated manner not explored here. Stimuli, like spikes arriving from an upstream node, could thus toggle the system these two behaviors. Figure 5. Burst dynamics explained with three-dimensional V m trajectories and bifurcation diagrams (for bursts in Figure 4A with I maxpump = 95 mA/cm 2 ). Note: x-axis scales are different is each plot. (A) The first burst of spikes in Figure 4A plotted as a function of time are now plotted in 3-D as a function of E Na and E K (blue arrow: direction of V m trajectory). (B) Bifurcation diagram for fixed E K and I pump as per Figure 4A pink star; E Na is the slowly varying parameter. Lines and points are labeled (see also abbreviation list) and have the same meanings in (C). For E Na at its Figure 4A pink star value (40.2 mV), the only stable solution is a periodic orbit (pink oval), with each cycle corresponding to a spike (in a burst). The oval-loop symbolizes one cycle of V m oscillation at fixed E Na . During a burst the E Na decline shifts orbits leftward into the bistability regime (pink area) where there exist two stable solutions: a limit cycle and a stable fixed point. (C) Bifurcation diagram for E K and I pump fixed as per Figure 4A green dot.
For E Na at its Figure 4A green dot value (37.27 mV), the system (large green dot) is within the bistability regime (green area). V m , attracted by this fixed point, has STOs (drawn as the green loop) until E Na increases and superthreshold-oscillations (spikes) return. The PD region corresponds to period-doubling bifurcations. (D) Two-parameter phase diagram for E K and E Na . Pink solid and dashed curves represent LP (saddle-node bifurcation) and HB, respectively, when I pump is fixed as in B. The green solid, dashed, and dash-dotted curves represent LP, HB and PD, respectively, when I pump is fixed as in C. With varying I pump the bistability regime shifts from the pink to the green area (the zone with both colors is the overlap of these two areas). The gray area between two green dash-dotted curves is a zone with PD bifurcations. The black loop shows E K and E Na orbits during a burst. doi:10.1371/journal.pcbi.1002664.g005

Spontaneous bursting activities occur for small and for large LS injuries, while intermediate intensity injuries elicit tonic firing.
In the case of bursting, time-varying E Ion values due to pump activity provide a slow almost periodic element in the system; E Ion values approach zero during the burst and grow during the interburst interval. Slow-fast decomposition of the full dynamics reveals that these slow E Ion bring the fast spiking dynamics into and out of a limit cycle (tonic firing) regime. 5. STOs can appear with single injured Nav population or for complex CLS injuries. STOs that bracket bursts of firing are associated with the complex eigenvalues governing dynamics in the vicinity of the lower branch of fixed point equilibria. Burst onset is associated with a subcritical Hopf bifurcation while offset is organized around a saddle-node bifurcation of limit cycles. Period-doubled solutions exist near this latter bifurcation, and the limit cycle can collide with these solutions, producing abrupt changes in the interval between successive spikes. Bursting solutions tend not to repeat over many bursts, and thus appear to be chaotic (not analyzed in detail). 6. Including small levels of Gaussian white noise in simulations causes burst characteristics to fluctuate and yields more longlived STOs, consistent with the reported propensity of injured axons to exhibit STO phenomena. 7. Inhomogeneous CLS is sufficient to explain a rich variety of axon excitability behaviors. Combinations of CLS-injured Nav channels, all exhibiting realistic m 3 activation behavior, yield flattened g Na (V) curves similar to those of other models that, of necessity, include at least one g Na with m 1 activation.

Nav-CLS and pathophysiological Nav leak
Regions of membrane with mild Nav-CLS have a ''subthreshold'' persistent g Na (a Nav-leak conductance) due to g window operating abnormally close to the cell's normal resting V m ; when only a small fraction of Nav channels are affected, the cell's g window plot has a shoulder as in Figure 4C; such plots strongly resemble coupling-induced subthreshold persistent conductance of pacemaker neurons [56]. Like that physiological pacemaker, our injury model includes no slow mode (or non-inactivating) g Na . Calculations [3] show, in any case, that at voltages near typical AP thresholds, left-shifted I window would dominate over equivalently left-shifted non-inactivating I Na . In much of this study, we mimicked ''smeared'' extents of bilayer damage via Gaussian distributions of Nav-CLS, reflecting the notion that nodal axolemmal injury is unlikely to exhibit a unique intensity. Smear, though expected in situ [2,8], is not actually required for neuropathic activity. Use of assorted smears confirmed that CLS-fuzziness in no way abrogates dysfunctionality and it demonstrated the robustness of Nav-CLS as a model for leak. Nav-CLS injury thus provides a minimalist explanation for multiple spontaneous modes of firing in which APs, bursts and subthreshold oscillations (STOs) occur in various patternings. In severe cases, it also predicts depolarizing block [11]. We used bistability diagrams and various bifurcation analyses to help clarify the mechanistic underpinnings of the Nav-CLS-induced dynamics.
A mildly injured excitable cell cannot survive for hours or days (the time course, e.g., of diffuse axonal injury) without operational Na/K pumps. In electrical models of axon injury, pumps do not only maintain ion gradients, they also contribute to slowly varying excitability patterns of the damaged-but-viable axons [11,29,30]. Pursuing this further here, we found the following: mild Nav-CLS induces STOs through a gradual reduction of firing threshold mediated by slow E Ion dynamics whose characteristic times depend on pump-mediated fluctuations of extra-and intracellular ion levels. Threshold voltages for firing are not fixed but depend on E Ion . E Ion fluctuations induced by Nav-CLS produced, in a bistable fashion, both STOs and bursts of full-sized APs. With mild Nav-CLS injuries, adding Gaussian white noise to the current-balance  Figure 6C, but here, noise (s = 0.4 mA/cm 2 ) causes V m to randomly visit quiescent or repetitive firing states, hence producing noisy bursts. doi:10.1371/journal.pcbi.1002664.g007 equation induced erratic bursting, or, where injury was already causing sustained bursting, noise influenced the rates and durations of the ectopic signals. These effects of noise are consistent with those seen in a recent computational study of neuropathic firing in a hyper-excitable system with more ion channel species than used here, but without pumps [57]. In the following we discuss and provide more context for the above results. Problematic rhythmicity patterns associated with mild axonal injury Trauma-induced [58] and ischemia-induced [59] Nav leak triggers Ca 2+ -excitotoxicity, the proximate cause, in diffuse axonal injury and oxygen glucose deprivation, of axon demise. A mechanistic understanding of early stages of Nav channel dysfunction, when salvage is possible [60], is needed. Although pathological TTX-sensitive Na + leaks are attributed to increased ''persistent I Na '' (e.g. [61,62]) or slow-inactivating I Na [63] it is not known how membrane damage elicits chronic Nav-leaks. A molecular-level understanding is lacking as to what constitutes leaky native Nav channels. A common Nav-leak culprit in diverse sick excitable cell conditions is, we propose, Nav-CLS [3]. For recombinant Nav1.6, this phenomenon occurs as an immediate consequence of mechanical membrane damage [8,11]. Any fraction of I window (an attribute of fast g Na ) that has become left-shifted acts as a leak [3,8] that could trigger positive feedback excitation.
In previous modeling [11], we established that mild axonal Nav-CLS would effectively dissipate [Na + ] gradients, bringing on excitotoxicity. Total nodal g Na was modeled by an intact fraction plus a second fraction with a discrete level of CLS. Various ''Nav-CLS injuries'' were imposed, V m (t) was plotted as new steadystates were attained, then excitability was probed with current injection. Mild Nav-CLS caused spontaneous APs (tonic or damped trains). With Na/K pump activity included, it caused AP bursts whose specifics varied with compartmental volumes. From temporal patterns of [Na + ] gradient rundown, it was evident that Nav-CLS type injuries could cause axon demise by overwhelming an axon's capacity for ion homeostasis. Injured nodes, moreover, generated pernicious ectopic signals and impaired transmission fidelity during saltatory propagation. Pathological activity patterns, plotted as a function of the extent (in mV) of Nav-CLS, and of the fraction (0-1) of affected nodal Nav channels, fell into several regimes.
What can we infer about responses to stimulating inputs when there is Nav-CLS injury? Effects of stimulation, a multi-faceted issue, was addressed for nodes with fixed reversal potentials in a preliminary study [11] involving two Nav populations, one intact and one with a specified CLS injury. Spontaneous firing occurs over a range of damage, then beyond a critical threshold level (LS c ) the system is quiescent. Stimulating with constant current yields a larger maximal firing rate but, due to depolarizing block, quiescence occurs with less damage (LS,LS c ). Effects of stimulation are currently under study for bursting nodes (model with pumps as in Figure 5). Because bursting is a globally attracting  Table 2. Steady-states for intact nodes (V m = 259.9 mV) as I maxpump varies (units as in Table 1).  behaviour, a single spike as a form of input may transiently change a bursting cycle but not its steady state. For the regime we considered, bifurcation analysis reveals that a single spike cannot ''trigger'' bursting behavior.

Nav-CLS and neuropathic pain
Though ectopic signals in vivo are generated at both axons and somata [64], peripheral neuropathic (as distinct from nociceptive) pain mechanisms are mostly studied via ectopic activity in damaged somata such as those of primary mechanoreceptor afferents [62,65,66]. A closely-related issue of relevance here is that of low threshold afferents, dubbed ''algoneurons'' [67]. These are nonpain sensory neurons that can serve double functions, signalling pain only when peripherally hypersensitized due to trauma, inflammation and the like. Spontaneous activity patterns described here are like those of both damaged mechano-afferents neurons and hypersensitized algoneurons. This is a post hoc assessment; we monitored behaviors across parameter space, but did not adjust parameters to mimic particular neuron activities. Clearly, it would be interesting to determine if ''hypersensitization'' in algoneurons corresponds to Nav-CLS in damaged peripheral endings.

Mechanisms for generating STOs
Emerging from that ''injury''-modified Hodgkin-Huxley axon (two g Na sub-populations and the Na/K pump) was a critical feature of the electrical dysfunction of injured excitable membranes, abnormal STO phenomena [11]. We lacked, however, a dynamical understanding of the phenomenon and its robustness; given the prevalence of STOs in neuropathic firing (plus physiological correlates) this represented a significant gap. Typically, STOs and injury have been modeled with non-inactivating and/ or slow-gating g Na . Our model -Nav-CLS -includes only a fastgating g Na to which a simple biophysically-justified (and biophysically-explicable) modification has been applied at varying intensities. Here, using smeared-intensity versions of Nav-CLS to model injury, we generated diverse STO phenomena.
The ''spared but impaired axons'' of injured tissue are targets for Nav antagonists (see [21]) and such mildly-injured axons were our focus here. We examined sustained spontaneous (ectopic) behaviors, extending our analysis first to three g Na sub-populations (i.e., for otherwise identical g Na mechanisms, the Nav population was given three extents of Nav-CLS), then to various Gaussiandistributed (''smeared'') Nav-CLS populations. Realistically untidy, but conceptually simple, all these Nav-CLS injury variants exhibited the key components of the injury repertoire. STOs arise with ( Figure 8D) or without smear. Along with the basic spontaneous activities (APs and bursts), injured nodes spontaneously produced sustained STOs ( Figure 3) and STOs as parts of bursts consisting of clusters of spikes and STOs (Figure 4).
Pump/CLS-leak: a slow ''component'' of an oscillatory system Na/K pump turnover rates are sensitive to substrate concentrations, membrane potential [68] and lipid packing [39,40]. Quiescent (''at rest'') axons maintain a steady basal pump activity that counteracts all system leaks. In our model I pump is always substantially smaller than I Na and I K (,22-fold in Figure 4A) and pump rates oscillate significantly more slowly than APs (e.g., Figures 8C,D, ,2 orders of magnitude slower). In our simulations of injury, small amplitude slowly oscillating I pump plays an essential role in mediating spontaneous V m excursions (Figures 3, 4 and 9). Larger magnitudes for I maxpump (see Figure 3D and Figure 9) accelerate the efflux/influx of Na/K ions. This was associated with dynamical changes, from RP to STOs to bursts or tonic APs. Figure 3D in particular is a bifurcation diagram; it characterizes the transition from quiescence to tonic firing as I maxpump increases. This shows clearly the increase in the amplitude of voltage excursions V m and the decrease in firing rate. Figure 9 shows the effect of I maxpump in the bursting regime. As I maxpump increases, we go from quiescence to bursting. All characteristics of the burst summarized in Figure 4 shorten with this increase.
In the injured nodes, therefore, a reduced I maxpump (due, say, to ATP depletion and/or to a more fluidized pump-embedding bilayer) would modify ectopic excitation patterns. In line with this, in endogenously bursting snail neurons, inhibiting the pump shortens burst periods and increases spike frequency [69], though unlike injured nodes, pacemaker neurons have multiple slow processes that confer resilience on specific oscillatory patterns. Dynamic transitions of I pump are seen in a number of rhythmically active cells including sinoatrial node cells [44] and trigeminal motoneurons [70].
Insofar as Na/K pumps show greater ATPase activity in membranes with greater lipid packing order [39,40], membrane injury that caused Nav-CLS should simultaneously impair pump activity. To better understand ectopia arising from injured nodes or axon initial segments, it would be helpful to know if Na/K pump inhibitors (e.g. ouabain, but also generic membrane fluidizers [43,71]) reduce ectopic burst firing in mildly damaged axons. In this regard, we note that peripherally-applied ethanol (shown recently to inhibit neuronal Na/K-ATPase [72]) induces a ''feeling no pain'' condition when injected locally to relieve peripheral trigeminal neuralgia [73].
Bifurcation analysis: Nav-CLS/pump STOs help trigger spikes by decreasing firing threshold Although neuropathic STOs trigger APs in many neurons [58,74], STOs are not principally a neuropathic phenomenon. In the CNS [75], STOs are a filtering mechanism by which neurons decode oscillatory inputs from other neurons in a frequency dependent manner [76]. Disparate underlying mechanisms for STOs have been noted, including Ca channel I window [77] and HCN-channel based I h [75], and, of particular interest here, a fast g Na -based [56] subthreshold persistent I Na [78]). Our mild Nav-CLS-based I window resembles that mechanism, but generates STOs and AP bursts only in conjunction with the electrogenic pump. STOs arising spontaneously at nodes are, by definition, pathologic, and given the dominance of nodal conductances by g Na (specifically, Nav1.6 channels) would likely be simpler in nature than the STOs elicited (by current stimulation) in multi-Nav isoform sensory neuron somata [21] or in double-Nav nociceptive terminals [23].
STOs simulated here involve a small-amplitude rhythmic mechanism that develops with mild Nav-CLS I window (Figures 3  and 4) operating near what would normally be the RP level. Bifurcation analysis of this system enables us to propose that slow dynamics of the firing threshold are responsible for the timing of spikes. In a system where the E Ion vary, the firing threshold varies too, gradually decreasing in the interval between successive STOs (Figure 9). This is broadly consistent with experimental findings of a decreased spike rheobase in injured dorsal root ganglion neurons [79] and, more general terms, seems consistent with the injury-induced hypersensitization of algoneurons [67]. Models involving fixed E Ion cannot exhibit spontaneous oscillations of firing threshold, but in vivo, and especially in fine processes (e.g., Nav-rich nociceptive endings), fluctuating [Ion]s and pump currents are expected once channels activate.
Leak/pump counterplay and the oscillating threshold: implications for healthy membranes Analyses of the electrical excursion diversity seen in excitable cells gradually going through their death throes have, historically, been fruitful in providing sketches of the behavioral diversity of healthy excitable cells. For example, the Hodgkin [80] classification of various spiking behaviors in isolated crayfish axons and the Morris-Lecar [81] analysis of giant barnacle myocyte voltage oscillations provide the springboard for analytical descriptions of oscillatory types in mammalian interneurons [82]. Here we simulated sick cell behaviors in a context where large fluxes into small volumes made it crucial to include the contributions of the electrogenic Na/K pump. From this emerged an STO mechanism that triggers bursts of spikes in a manner not, to our knowledge, previously recognized, i.e., the oscillating threshold. Since threshold occurs at V m values where fluxes unbalance, threshold varies with E Ion variations. An ultrasimple biological exemplar of this point is a g Ca -based biphasic AP (Figure 3c of [83]) that shows different thresholds for its depolarizing and hyperpolarizing excursions and that is modeled simply by letting [Ca 2+ ] intracellular accumulate with g Ca activation (Figure 3d and Equation (8) of [81]).
Here we explicitly considered injury-induced Nav-CLS in a setting where Nav1.6 would be the exclusive Nav isoform. But Nav-CLS may have resonances for healthy excitable membranes with two or more Nav isoforms with different voltage midpoints, such as seen in axon initial segments [33] and nociceptive nerve endings [23]. Moreover, even with a single isoform, nonhomogeneous physiological modulation (e.g., from inhomogeneous lipid packing density) could yield CLS-like situations that might affect tuning (e.g,. in the cochlea [83,84]). As in vivo lipid imaging becomes refined, it may be possible to detect lipidpacking changes before and after injurious stimuli at nodes or axon initial segments. This is important given that the activities of voltage-gated channels and the Na/K pump are both modulated by lipid packing [37,39]. The kinetic phenotypes of the membrane proteins simulated here are, of course, constrained by amino acid sequence, but they also vary according to the bilayer mechanical milieu. Our work here suggests how, in the face of bilayer injury, the sensitivity of Nav channels to bilayer mechanics can be perilous. By the same token, it may spell opportunity, with developmentally controlled CLS-type kinetic modulations spawning behavioral amplifications beyond what would initially be expected from a limited set of membrane proteins.

Node of Ranvier with variable amounts of membrane damage
Axonal voltage excursions are modeled at an individual node of Ranvier with the Hodgkin-Huxley model [85] using the values listed in Ochab-Marcinek et al [86]. We start with the basic equation for membrane potential (V m ): where C is the nodal membrane capacitance. The total I Na for all Nav channels and the potassium current I K are given as: I K~gK (V m {E K ), where g K~gK n 4 : ð3Þ g Na and g K are the maximal conductances of the Nav and Kv channels, respectively; E Na and E K are the sodium and potassium reversal potentials, respectively; and n 4 gives the steady-state open probability of potassium channels. By the HH formulation, gating variables m, h, and n evolve according to: The forward (a m and a h ) and backward rate functions (b m and b h ) describe the first-order transitions between the deactivated and activated states of activation (m) and inactivation (h) processes in Eqs (4)(5). They are determined by the membrane potential V m : The rate functions a n and b n of potassium gating variable n in Eq (6) also depend on V m : Experimental findings for recombinant I Na from nodal type Nav1.6 channels [8] show that, with mechanical membrane damage (''injury''), transient I Na at a given V m irreversibly accelerates because the activation and steady state inactivation (availability) processes undergo irreversible hyperpolarizing (''left'')-shifts. Consistent with this, the steady-state product of activation and availability (m 3 h(V m ) tR' ), also called the window conductance, also left-shifts. This injury response we term ''coupled left-shift'' (CLS). This CLS can be easily modeled within the above HH model by shifting the voltage dependences of the activation (m) and inactivation (h) functions by LS, i.e replacing V m by (V m +LS) in Eqs (7). In other words whatever used to be observed at a given V m in any function of m and h is now observed at (V m 2LS). This CLS is the basis of our model for axon damage. In patch clamp experiments the activation/inactivation time course data are explained with a single LS value for each degree of injury [8]. However, axon injury is spatially non-homogeneous [2] and is therefore likely to result in spatially inhomogeneous CLS. To account for this possibility we introduce fractions f i of Nav channels in damaged membrane, with left-shift intensity LS i (in millivolts). A Gaussian distribution of LS i seems a realistic model of axon injury, hence we explore various Gaussians. But discrete populations with a small number of LS values are more amenable to a thorough analysis, such as the bifurcation analysis done on a three population system exhibiting bursting dynamics (see Figure 5).
Specifically for a population of left-shifted Nav 1.6 channels {f i , LS i } i = 1,N , for each f i of intensity LS i , independent activation variables m i (i = 1,2,…N) and inactivation variables h i (i = 1,2,…N) are introduced. The total I Na for all Nav channels is now given as: Three classes of Nav-CLS are considered here: (a) two g Na populations: one intact, one left-shifted by a specified fixed amount (as in [11]); (b) three populations: one intact, two with discrete LS values; and (c) assorted versions of LS of all or specified fractions of the total Nav population, but with f i (LS) following Gaussian distributions to mimic the expectation that nodal injury would produce ''smeared'' (as opposed to discrete) extents of ''Nav-CLS''. At the node of Ranvier, g Na dominates over g K , but previously [11] we briefly discussed the effects of also applying a LS to the gating of the Kv channels.
The following features are also included in the model. Leakage currents [29], I Naleak and I Kleak , account for the background permeabilities needed to bring the basal system to a quiescent steady-state V m with an active pump present, along with the usual unspecific leakage current I leak . These leak currents take the form: where conductances g Naleak , g Kleak and g leak are constants. Parameter values are summarized in Table 1.
With an active Na/K pump, the model system has dynamic intracellular and extracellular Na + and K + concentrations, denoted by [Na + ] i , [Na + ] o , [K + ] i , and [K + ] o . The Na/K pump balances the accumulation and depletion of ions in the intra-and extracellular space by exchanging 3 internal Na + for 2 external K + ions. Thus [K + ] o and [Na + ] i determine the pump amplitude [29,87]: where I max pump is the maximal current generated by the pump, K M K and K M Na are Michaelis-Menten kinetic constants, and the Na + and K + currents flowing through the pump are I Napump = 3I pump and I Kpump = 22I pump . Ion fluxes across the node of Ranvier membrane cause the ion concentration changes: where F is the Faraday constant, A is the surface area of the nodal membrane, and Vol i (or Vol o ) is the intracellular (or extracellular) volume of the node of Ranvier under study (we take Vol i~V ol o for simplicity). Ion concentration dynamics alter the reversal potentials E Ion in Eqs (2-3) by the Nernst equations: Model equations are solved using a 4-th order Runge-Kutta scheme with 10 23 ms as time step in MATLAB and bifurcation diagrams were plotted using XPPAUT.