Stimulation-induced ectopicity and propagation windows in model damaged axons

Neural tissue injuries render voltage-gated Na+ channels (Nav) leaky, thereby altering excitability, disrupting propagation and causing neuropathic pain related ectopic activity. In both recombinant systems and native excitable membranes, membrane damage causes the kinetically-coupled activation and inactivation processes of Nav channels to undergo hyperpolarizing shifts. This damage-intensity dependent change, called coupled left-shift (CLS), yields a persistent or “subthreshold” Nav window conductance. Nodes of Ranvier simulations involving various degrees of mild CLS showed that, as the system’s channel/pump fluxes attempt to re-establish ion homeostasis, the CLS elicits hyperexcitability, subthreshold oscillations and neuropathic type action potential (AP) bursts. CLS-induced intermittent propagation failure was studied in simulations of stimulated axons, but pump contributions were ignored, leaving open an important question: does mild-injury (small CLS values, pumps functioning well) render propagation-competent but still quiescent axons vulnerable to further impairments as the system attempts to cope with its normal excitatory inputs? We probe this incipient diffuse axonal injury scenario using a 10-node myelinated axon model. Fully restabilized nodes with mild damage can, we show, become ectopic signal generators (“ectopic nodes”) because incoming APs stress Na+/K+ gradients, thereby altering spike thresholds. Comparable changes could contribute to acquired sodium channelopathies as diverse as epileptic phenomena and to the neuropathic amplification of normally benign sensory inputs. Input spike patterns, we found, propagate with good fidelity through an ectopically firing site only when their frequencies exceed the ectopic frequency. This “propagation window” is a robust phenomenon, occurring despite Gaussian noise, large jitter and the presence of several consecutive ectopic nodes. Electronic supplementary material The online version of this article (doi:10.1007/s10827-014-0521-9) contains supplementary material, which is available to authorized users.

[S1] where i = 1, …, N is the node number. The last term on the RHS is the inter-compartment current and uses V 0 = V 1 and V N+1 = V N boundary conditions. The six first terms are transmembrane currents: the two currents (Na + and K + ) controlled by voltage-gated channels (VGC), the currents mediated by the population of Na/K pumps, and three ohmic leak currents (sodium, potassium and unspecific). We neglect any possible internodal sodium (Zeng et al 2008), although our results could extend to unmyelinated axon with sodium clusters (Zeng et al 2009).
VGC currents. The voltage-gated sodium and potassium currents were modeled using the equations of the Hodgkin-Huxley (HH) model (1952), with the usual kinetics for the Na-and K-gating: 3 Na, Na Na, where m i and h i model the activation and inactivation (availability) of the Nav channels while n i is the Kv channel gating variable. However, in this work, the HH model is augmented with time-varying concentration-dependent Nernst potentials E Na and E K given by where the indexes o and i denote outer and inner concentrations respectively. Time-dependence of the concentrations is described below. In the HH model, the three gating variables evolve in time according to , , The forward () and backward () rate functions are considered to be functions of the voltage V i only that can be found in (Hodgkin & Huxley 1952) where the indexes i = node number were dropped, here and in all following equations, to ease notation and to avoid confusion with i = inner. Previous work (Boucher et al 2012) argued that HH kinetics approximates results obtained with more evolved models (eg, Taddese & Bean 2002) in which the kinetic coupling of activation and inactivation (availability) in Nav channels is represented using a manystate allosteric Markov model. The HH model, however, presents the advantage of using much less computational time.
Pump currents and concentration time dependence. Even though the activity of Na/K pumps shows a small dependence on voltage, in normal conditions it depends mostly on the concentration of the substrate (Läuger 1991). Therefore, we model the pump current at each node using the Michaelis-Menten enzyme kinetics model, as was also done in (Kager et al 2000): where the exponents show up because the pump exchanges three inner Na + ions for each two outer K + ions.
The point of introducing pump currents is to investigate time-dependent ionic concentration changes. Na + and K + concentrations at each node evolve according to where  Na,pump = 3 pump ,  K,pump = -2 pump ,  X is the VGC current given by Eq. S2 or S3, F = N A e = 96472 C/mol is the Faraday constant, r = 20 cm 2 /μL is the surface to volume ratio of the node of Ranvier and the inner and outer volumes are considered equal for simplicity. Note that 20 is a rather high value for r, chosen to shorten computation time without affecting phenomena qualitatively : a cylinder of length L and radius R = 5 μm would have r = 4 cm 2 /μL and a real axon, where ions can effectively diffuse from a volume longer than the nodal length L, would have an even smaller r value. Using a more realistic value, phenomena reported in Figs 1, S3 and S8 would occur over a time scale about one order of magnitude larger.
Concentration changes given by Eq. S8, in turn, affect the Nernst potential that appear in Eq. S2 and S3. Since these concentration changes are orders of magnitude slower than HH kinetics, they are neglected in some simulations. We did so by setting r = 0 in Eq. S8, not by returning to a standard HH model with  pump = 0 and I X,leak = 0. This ensures that constant concentrations are a limiting behavior of our model, not a feature of an entirely different model; it also promotes comparability of the results.
Leak currents. The sixth term in Eq. S1 is the regular nonspecific leak current of the HH model, given by leak leak However, when pumps are added in the model, specific leak currents I X,leak are also required: indeed, the total Na + and K + currents, found between square brackets in Eq. S8, need to both be zero when the node is at rest. Since the pump current is never zero and since there is a window current involving I Na and I K even at rest potential, specific leak currents are needed to balance total Na + and total K + currents separately (Boucher et al 2012). They are given by All equations were solved with RK4 method with variable timestep, implemented on a PC using code in C language. The variable timestep was adjusted as suggested in (Press et al 1992): the 5 th order term was monitored and kept just below a 410 -12 target. This term was then added to the 4 th order (RK4) result. The raw code is made available as Online Resource 3.
Damaged node(s) of Ranvier. In our computational model, irreversible CLS of activation and inactivation of the Nav channels is readily modeled by shifting the voltage dependence of  m ,  h ,  m and  h by an amount LS (i.e. replacing V by V + LS in Eq. S6, in the damaged node(s) only).
Individual channels are not modeled explicitly; rather, we assume all channels in the damaged (i = 6  Q) node of Ranvier are left-shifted by the same LS value. In some figures, nodes Q±1 were also damaged, with the same LS value as for node Q. We do not consider situations where a distribution of LS values occurs at a same node, even though we expect our results to be easily extended to those situations.
Parameter selection. Unless stated otherwise, all parameters in the above equations take the values listed in Online Resource 1. Of those, the classical parameters (C, g Na , g K ) were chosen identical to the ones in (Ochab-Marcinek et al 2009). The three parameters g Na,leak , g K,leak and  max needed to be chosen to cancel the total Na + current and the total K + currents in the resting state. These two constraints for three unknowns left one degree of freedom. The last constraint that we considered is qualitative: the node should respond readily to a stimulation current which is a small fraction of the total Na + current observed during an AP. This consideration led us to select a reasonably small I max . Finally, three remaining parameters (, N and the stimulation amplitude V stim , to be defined) will now be discussed in greater detail.
Simulations with a damaged axon will be compared with identical simulations on an intact (control) axon. This intact axon should behave normally: each time it is stimulated, it should fire an AP and each AP fired should be propagated faithfully in all nodes. Also, behavior should be independent of N, the total number of nodes. The last three parameters, , N and V stim , were chosen to ensure this expected behavior ( Fig. S1).
Stimulation at the initial node of the axon. To allow for the study of both periodic and aperiodic spike trains, node 1 has to be stimulated with timed delta functions, each causing an instantaneous V stim discontinuity in V 1 . To select an appropriate value for V stim , we fed periodic stimulation of various heights V stim and frequencies f stim , using r = 0 and  = 0, and kept track of the ratio of firing frequency to f stim . Fig. S1a shows (in green) the range where this ratio is 1:1 and shows the value selected, V stim = 25 mV. When r > 0 is restored, the system behaves normally for a prolonged stimulation period (Fig. S1b). When propagation is restored ( > 0), the limit of the 1:1 zone slightly shifts towards lower frequencies (see "stim failure" limit on Fig. S1c).
Outside the green range on Fig. S1a, too large V stim caused period doubling while too large stimulation frequencies caused partial response, with or without phase locking (other than 1:1). Interestingly, partial responses involved a Farey sequence of phase locked zones, surrounded by zones of aperiodic behavior, much like what Feingold et al (1988) obtained in an externally driven van der Pol oscillator (see also Aihara et al 1986;Kaplan et al 1996). Also, the green zone boundaries were slightly dependent on E Na and E K , as shown by the curves in Fig. S1a. The tested values of E Na and E K are typical of the depletion that occurs during a prolonged (300 ms) stimulation (red and green curves on Fig.  S1b).
Propagation along the axon and length of the axon. In Eq. S1, the saltatory propagation along the computational axon was modeled by considering each internode as a single ohmic conductance  linking two consecutive nodes of Ranvier. This is equivalent to assuming that the myelinated internodes have negligible capacitance. Consequently, their RC charging time constant is so small compared to the node of Ranvier's that the nodes can be assumed to directly charge one another. This approach was previously used in ( -20 20 60 0 300 600 900 1200 1500 To select an appropriate value for , we fed periodic stimulation (V stim = 25 mV) of various f stim , using r = 0, to axons having a range of  values, and kept track of the ratio f N /f 1 of the firing frequencies of node N and node 1. Fig. S1c shows (in green) the range where this ratio is 1:1 and the value selected,  = 0.3 mS/cm 2 . (We separately tested that  was still low enough to prevent subthreshold oscillations to be propagated.) For short axons, the limit of the green zone on Fig. S1c was found to be N-dependent, but not above N  10, the value we have thus chosen. At the right of the "stim failure" line, f 1 < f stim (zone of "stim failure") while at its left, despite the fact that f 1 = f stim , there is a zone where f 10 = f 1 (green zone) and another where f 10 < f 1 (zone of "propagation failure"). Therefore, at  = 0.3 mS/cm 2 , propagation failure can be distinguished from stimulation failure. As for stim failure, propagation failure involves a Farey sequence of phase locking zones. When r > 0 is restored, the system behaves normally for a prolonged stimulation period (Fig.  S1d). It does so for all firing frequencies  85.1 Hz.

Additional Results
Onset of ectopic activity. Fig. S2 shows the limit between quiescent and ectopic behaviors as a function of E Na and E K . Diminishing E K reduces the K + driving force and increases the rest potential, bringing it closer to firing threshold. Diminishing E Na partially inverts that effect, but occurs too slowly to prevent the onset of ectopicity. Note that once ectopicity is triggered, further decline of E K increases f Q (see also Fig. S8).  Fig. S3b is the same obtained with nodes 5-7 damaged, while the middle and lower panels are, respectively, the values of E K and E Na at the time of the next stimulation following the onset of ectopic activity. Note how ectopic activity is triggered when a limit value of E K is reached, whatever the value E Na reached at the same time (this is mostly obvious on low f stim curves).
-71.45  Fig. S1. Parameter tuning of the intact (control) system. (a) In green : range of stimulation strengths (V stim ) and stimulation (stim) frequencies for which the isolated first node fires each time it is stimulated. The limit of the green zone (bounded by the curves) is recomputed for three fixed E Na and E K values (r = 0), with  = 0. Value selected and used throughout this article : V stim = 25 mV. (b) Typical response of isolated node ( = 0) to prolonged periodic stimulation (denoted by black bar), once r = 20 cm 2 /μL is restored (V stim = 25 mV). Note that this r value is chosen large to promote rapid depletion of the Nernst potentials. (c) In green : range of internode conductances () and stim frequencies for which APs fired in node 1 are faithfully propagated to all N nodes. Limit of green zone is recomputed for various axon lengths (N), with r = 0. Values selected and used throughout this article: N = 10 nodes and  = 0.30 mS/cm 2 . Note that the control system behaves normally for all stim frequencies <85.1 Hz. (d) Typical response to stimulation of the 10-node intact axon, once r = 20 cm 2 /μL is restored (stim as in (b),  = 0.30 mS/cm 2 ). The top node (1) is the one being stimulated. Fig. S2. Ectopic frequency as a function of E Na and E K , drawn for LS = 3 mV, r = 0, t LS = t  = 0 and no stimulation. White zone = not ectopic. Note that diminishing E K brings the system into ectopicity, while diminishing E Na does the opposite.
Evaluation of the settling time in Fig. 3a. To draw Fig. 3a, the following procedure was followed: 1) Delays between the n th output spike in control and damaged axons were computed. 2) We tested if those delays converged towards a constant (potentially non zero) value (a small fluctuation was tolerated). If not, the settling time is infinite (longer than the simulation duration) and procedure stops. Otherwise, 1:1 phase locking is considered detected.
3) If phase locking is detected, we found the instant t lock when consecutive delays started differing by less than 0.5 ms (we ensured that this difference didn't re-increase later by conducting the search with decreasing t). Once this instant was located between two consecutive spikes, we linearly interpolated the precise t lock . 4) The settling time was then equal to t lock  t first , where t first is the instant when the first AP reaches node 10 in the control system. 5) Finally, the result of step 4 was found to be dependent on the phase of the ectopic node at the instant when stimulation begins (remember: since node 6 is ectopic, it already fires at frequency f Q before stimulation begins). To eliminate the influence of this factor, steps 1-4 were repeated 5 times, using equidistant phases obtained by changing t LS . The five results were then averaged to yield a point of Fig. 3a.
VP metric (output infidelity). The Victor-Purpura (VP) metric (1996) allows to compare two spike trains (S a and S b on Fig. S4) by evaluating the minimal series of elementary steps required to transform one of these trains into the other.
Two elementary steps are possible, a "cost" being associated to each: adding/removing a spike (cost = 1) or shifting a spike by Δt (cost = qΔt, where q is an arbitrary parameter). All possible ways to convert S a into S b using these steps are computed ( Fig. S4 shows a single one of them) and the minimal cost is considered a measure of the difference between S a and S b , their "VP distance".  Using the VP metric to evaluate output infidelity takes into account both settling time (during which a number of spikes are added) and phase shifts, thus producing a reasonable measure of output infidelity. In Fig. S5a, we show results for a range of q values.
As shown in (Victor & Purpura 1996), this "VP distance" matches all the mathematical conditions to qualify as a metric. Therefore, if one imagines an abstract space where each possible spike train is a point, the cost obtained can be seen as a distance between these two points.
Effect of parameters on output infidelity. Fig. 4 used q = 0.2 ms 1 and 500 ms of stimulation / comparison time. Fig.  S5 shows that one of our main results, namely the existence of a "propagation window", is independent of these two choices: except for very large q values, a propagation window is always observed (Fig. S5a). Also, increasing the duration of the stimulation (during which the output trains are compared), simply reduces the relative importance of the settling time, thus deepening the minimum band corresponding to the propagation window (Fig. S5b).
Note that the q = 0 curve on Fig. S5a shows zero output infidelity in the propagation window, meaning identical spike count. The non zero infidelity seen for identical f stim on Fig. 4a (where q > 0) is therefore due only to timing differences, which are in turn due to the settling time and to the phase shift that we described in Fig. 3.
Effect of local conductance. While Fig. 5 shows the effect of local increase or decrease in radius when LS = 4 mV, Fig.  S6 shows the same for other values of LS.

Generation of jitter and noise.
To produce each of the jittered inputs, we first generated a deterministic periodic spike train. Then, for each of the spikes in this input train, we used a random number generator to select a "displacement" from a Gaussian distribution of fixed standard deviation "stdev" (see legend on Fig. 6a). A total of 50 spike trains were thus obtained and saved, each with random displacements selected from the same Gaussian distribution. Each of those 50 input trains was then fed both to the control and damaged axons, producing each time two outputs that were compared using the usual VP metric to measure output infidelity. (This way, compared outputs are always affected by the same random effects.) Finally, the 50 infidelity values were averaged to produce one point on Fig. 6a. where W n is the n th number from that node's list, A is the noise amplitude (see legend on Fig. 6b) and h is the integration time step. Note that when noise was present, we turned off the adaptative time step and used h = 0.001 ms. As To help visualize the relative importance of jitter with respect to quasiperiod or the importance of noise relative to voltage excursions, Fig. S7 shows an example time series of each.

Output infidelity with free running ionic concentrations.
When r = 20 cm 2 /μL, pumps cannot speed up enough to prevent E Na and E K from declining. At first, the Nernst potentials are nearly unchanged. Therefore, mildly damaged axons, for which the "settling time" is the shortest, have no trouble phase locking 1:1. However, their ion gradients eventually deplete, f Q increases, the 1:1 phase locking condition is lost. Therefore, for the total stimulus duration (500 ms), the propagation window is shallower. Fig. S8 shows this, both with a time series and with an output infidelity diagram. in Fig. 6a, the same lists were used both for the intact and damaged system. The outputs were compared using the VP metric and the 50 output infidelities were averaged to produce a point on Fig. 6b.