Internodal sodium channels ensure active processes under myelin manifesting in depolarizing afterpotentials

ARTICLE IN PRESS

Journal of Theoretical Biology 235 (2005) 451–462 www.elsevier.com/locate/yjtbi

Internodal sodium channels ensure active processes under myelin manifesting in depolarizing afterpotentials Alexander G. Dimitrov Centre of Biomedical Engineering, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Bl. 105, Sofia 1113, Bulgaria Received 28 July 2004; received in revised form 14 December 2004; accepted 28 January 2005 Available online 8 March 2005

Abstract The current opinion about processes in myelinated axon is that action potential saltatorially propagates between nodes of Ranvier and passively charges internodal axolemma thus causing depolarizing afterpotentials (DAP). Demyelination blocks the conduction that gives additional argument in favor of hypothesis that internode is not able to be activated by the existing internodal sodium channels. The results of our modeling study shows that, when periaxonal space is sufficiently narrow, saltatorial action potential is able to activate internodes. Low density of internodal sodium channels is sufficient to generate active internodal waves that slowly propagate from nodes towards corresponding midinternodes where they collide. The periaxonal width that stops internodal wave propagation (about 400 nm) is significantly larger than the highest value of the physiological range for this parameter (30 nm). Internodal activation is directly manifested as transmembrane internodal potential or as a full-sized action potential in periaxonal space where it can hardly be detected, and only as a small deflection in intracellular space. However, changes in the periaxonal potential cause transmyelin currents that lead to significant DAP. The shape and amplitude of DAP depends on myelin parameters and densities of internodal channels. Several technical parameters affect the results of calculations. Internodal spatial segmentation has to be sufficiently fine (at most 20 mm) for the model to be able to simulate internodal activation. We employ 338 internodal segments as compared with up to 21 used in previous models. Ionic accumulation together with related diffusive and electrical processes alter the calculated DAP amplitude. Inclusion of these processes in calculations demands such increase in the total number of segments that the numerical methods used up to now become unapplicable. To overcome the problem, an iterative implicit approach is proposed. It reduces a matrix of general type in multi-cable models to tridiagonal one and accelerates calculations considerably. r 2005 Elsevier Ltd. All rights reserved. Keywords: Myelinated axon; Multi-cable model; Internodal channels; Accumulation; Numerical method

1. Introduction Myelinated axon is known to conduct action potential (AP) saltatorially. Depolarizing afterpotential (DAP) was found in intracellular recordings at the end of AP (Barrett and Barrett, 1982; Blight and Someya, 1985). Barrett and Barrett (1982) suggested that the DAP has passive nature, i.e. it is mediated by charging the internodal axonal membrane capacitance through a resistive current pathway beneath or through the myelin Tel.: +359 2 979 36 10; fax: +359 2 723 787.

E-mail address: [email protected]. 0022-5193/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2005.01.024

sheath. To test this hypothesis, multi-cable models of myelinated axon were created (Blight, 1985; Halter and Clark, 1991). In these models, the myelin was separated from the axolemma by a cable corresponding to periaxonal space (Fig. 1), thus the state of the internodal axolemma could be incorporated into them. Depolarization of a small portion of the internodal axolemma near the nodes of Ranvier (Fig. 2a) was sufficient for the appearance of a small declining DAP in the modeled AP (see also Nygren and Halter, 1999; Stephanova, 2001; Stephanova and Bostock, 1995; Zhou and Chiu, 2001). Nevertheless, the hypothesis about the exclusively passive nature of DAP is not yet confirmed. Moreover,

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

452

40

Extracellular space (Cable 3, V3=0) Intracellular space (Cable 1, V1) NODE (1.5 µm)

…

MYSA (4 µm)

FLUT (75 µm)

STIN (1800 µm)

···

Axolemma

Transmembrane potential (mV)

Periaxonal space (Cable 2, V2)

0

-40 STIN

-80

Myelin wraps

(a)

0

1

2

3

4

5

2 3 Time (ms)

4

5

Volume for ion accumulation in MYSA 40

···

Fig. 1. Myelinated axon geometry and notation used in the model description. The NODE, MYSA (myelin sheath attachment), FLUT (fluted internodal), and STIN (stereotype internodal) regions in the axon are used as defined by Berthold and Rydmark (1983, 1995). Dashed lines separate regions of specific geometry. Dashed-dotted lines reflect the non-cylindrical shape of the FLUT region. Numbers in brackets are lengths of the corresponding regions.

using this hypothesis to elucidate the nature of DAP, David et al. (1995) concluded that it must be modified to take into account axolemmal potassium conductance that remained activated following the action potential. According to experimental results, however, DAP can initially increase and then decrease (David et al., 1995). McIntyre et al. (2002) found that the passive mechanism alone is insufficient to reproduce experimentally measured DAPs. An active contribution to DAP is also necessary and these authors introduced additional conductance in the nodal membrane, attributing it to persistent sodium channels. However, they mentioned some conflicting results between the density of persistent sodium channels and some characteristics of supernormal period in sensory nerve fibers. The inclusion of fast potassium conductance in the internodal region of the model resulted in decrease in amplitude and duration of DAP. This moved some characteristics of the fiber excitability farther away from the experimental ones (McIntyre et al., 2002, Fig. 9). On the other hand, studies of changes in post-tetanic excitability and ectopic discharges in myelinated axon (Baker, 2000; Bostock and Bergmans, 1994) were interpreted as caused by potassium accumulation at the internodes of myelinated axon, since similar axonal behavior was found in conditions of high peri-internodal potassium in rat (David et al., 1993) and lizard (David et al., 1992). In addition, a large number of voltage-gated channels was found in internodal axolemma of myelinated axon (Burke et al., 2001; Chiu and Schwarz, 1987; Roper and Schwarz 1989; Shrager, 1995; Vabnick et al., 1997;

Transmembrane potential (mV)

Internode

0

-40

-80 (b)

STIN

0

1

Fig. 2. Transmembrane potentials across axolemma in various regions (NODE, MYSA, FLUT, STIN) of a myelinated fiber simulated with different number of internodal segments and identical internodal 1 sodium channel density (30 of the nodal one). No diffusive processes were taken into account. (a) Twenty-one segments. The potentials propagate decrementally into the internode. Similar decremental 1 potentials were obtained for low sodium channel density (434 as in Halter and Clark (1991)) irrespective of the number of internodal segments, and for low number (21) of internodal segments irrespective of internodal sodium channel density. No internodal wave was present. (b) Three hundred and thirty-eight segments. Internodal action potential induced by the saltatorial AP, propagated non-decrementally into the STIN region.

Vabnick and Shrager, 1998). The internodal sodium channel density was found to be within 2–6% of the nodal one (Burke et al., 2001; Chiu and Schwarz, 1987; Shrager, 1995; Vabnick et al., 1997; Vabnick and Shrager, 1998). It was insufficient for transmission of the AP through demyelinated fiber. However, under certain conditions of focal demyelination, continuous conduction through such a region was detected (Bostock and Sears, 1978; Shrager, 1995). Hines and Shrager (1991) and Shrager (1995) related this conduction with formation of interrupted extracellular layer around demyelinated portion of axon at earliest stages of Schwann cell proliferation and

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

attachment. The proximal nodes were thus loaded less by the demyelinated internode. This was tested on a two-cable model of Xenopus-myelinated fibers (Hines and Shrager, 1991) one internode of which was either demyelinated or partially covered by 1–2 wraps of myelin. The authors, however, connected the nodes with normal internodes through a very high resistance. As a result, the internodal sodium channels were no longer significantly activated when the internode was covered by myelin along its entire length. The authors, therefore, concluded that in normal conditions, the conduction in internodes was passive. Thus, the only function associated with internodal sodium channels in literature, was continuous conduction in demyelinated axons. Different mechanism for activation of internodal sodium channels could be suggested. The formation of extracellular layer restricts the volume conductor around internodal axolemma. If the membrane characteristics of axons are identical, then the narrower periaxonal space should have a higher resistance. This would shorten the length constant along the internodal axolemma. Then, potential profile should be concentrated over a shorter length, therefore demanding less charge for depolarization of a smaller membrane area. This could result in higher depolarization of normal internodal membrane as compared to that in a demyelinated axon that is surrounded by a larger volume conductor. This increases the chance for the normal internodes to be activated. These facts suggest that AP could be able to activate internodal axolemma covered by myelin. To describe the shape of potential profile of a short spatial length in details, fine segmentation of the simulated fiber is required. This may be a reason why the existing multi-cable models fail to predict the internodal activation. If the internodal activation exists, it should inevitably lead to significant changes in ions’ concentrations due to the small volume of periaxonal space. Therefore, to be correct, one has to model diffusive processes as well. The aim of this paper is to test the hypotheses on existence of active processes in internodal axolemma of myelinated axons and on their possible relations to DAP using a proper mathematical modeling.

2. Method We have to choose whether to make our own model that, in our opinion, would best describe a myelinated axon, or to follow an existing model. The drawback of a new model is that the results obtained could be considered as model specific. That is why, we choose the multi-cable model of Halter and Clark (1991) that has served as a basis for most of the contemporary models of mammalian myelinated fiber. Geometry of

453

the axon and the cables used are presented in Fig. 1. Halter and Clark (1991) consider nodes and distributed internodes as corresponding capacitances in parallel with leakage, fast sodium, and fast and slow potassium conductances. Myelin is considered as resistance and capacitance in parallel. Besides, the model of Halter and Clark (1991) is accompanied by the best description of numerical method used. Unfortunately, for large systems that involve diffusive processes, this method is not effective (Nygren and Halter, 1999). Since not only we would like to calculate diffusive processes, but also significantly increase the number of internodal segments, we had to develop a new numerical method. This method is presented in Appendix A. Starting with the two-cable model of Halter and Clark (1991), we have adapted it for our purposes. We introduced accumulation of all major ions (Na+, K+, Cl
, Ca2+) in all cables. This enables us to calculate the conductivity of the media based on actual ionic concentrations and corresponding mobilities. As the origin of the leakage current is still unknown, we could not take into account the accumulation of its carriers. The initial conditions specified by Halter and Clark (1991) were not at equilibrium and we allowed potentials and concentrations to evolve to equilibrium state. Therefore, the conductivities used in our simulations were similar but not identical to those used by Halter and Clark (1991). We started all calculations from the equilibrium state. In the simulations performed without accumulation, all diffusive currents were assumed unchanging in time. The full system of equations for electrical and diffusive cables is presented in Appendix B. We have used two different segment divisions: traditional, with 21 internodal segments (Halter and Clark, 1991); and detailed, with 338 internodal segments. The NODE, MYSA, FLUT, STIN, FLUT, and MYSA parts of an axon (Berthold and Rydmark, 1983) (Fig. 1) were represented with 1,4,4,5,4,4 segments in the traditional case, and with 1,4,75,180,75,4 segments in the detailed one. The latter corresponded to 1 mm long segments for MYSA and FLUT and 10 mm long segments for STIN. Halter and Clark (1991) represented the FLUT region through 4 segments with a step change in diameter at the MYSA-FLUT transition. We made this transition 5 mm long and changed the diameter there gradually. Unless otherwise specified, we used geometrical and functional parameters of the myelinated axon as specified in Table 2. Our functional parameters differ from those in Halter and Clark (1991) by internodal sodium channel
1  density which was changed from 0.23% to 3.33% 30 of the nodal value. In our simulations we varied the following parameters: segment length; width of periaxonal space; internodal sodium, potassium and leakage channels’ densities; myelin conductance. The simulated axon was symmetrical in respect of the

ARTICLE IN PRESS 454

A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

node of stimulation (having 21 nodes was equivalent to a fiber with 41 nodes), and had sealed end boundary condition at the other end.

120 Saltatorial AP 80

3. Results 3.1. Active waves propagating along an internode

3.2. Range of parameters within which the internodal AP exists The internodal sodium channel density has to be at least 1.5% of the nodal density for initiation of the internodal waves. To ensure propagation of the waves in the STIN region, the sodium channel density there has to be at least 0.5% of the nodal one. It decreased up to 0.3% of the nodal value when the potassium and leakage channel density in STIN were reduced 10 times. The periaxonal width has to be less than 400 nm. To simulate the internodal waves, the spatial segments in STIN have to be at most 20 mm long. 3.3. Internodal AP is not yet detected. Why? The current, caused by the internodal transmembrane potential, passes through both a narrow periaxonal high-resistance space, and a considerably smallerresistance intracellular space. Thus, the main potential difference was produced in the periaxonal, but not in the intracellular space. As a result, internodal AP could be seen as a full-sized negative wave at the background of the saltatorial AP in the periaxonal space (Fig. 3), but only as a small deflection in the intracellular space (Fig. 5a, deflections marked by arrows). This means, that in a large intracellular space, where the detections are easier, the internodal AP is practically invisible. On the contrary, the internodal AP is pronounced in a very narrow periaxonal space, where the correct measurements are hard.

40 Periaxonal potential (mV)

1 When the internodal sodium channel density was 30 of the nodal density and the number of internodal segments was 338, the model demonstrated clear excitation of internodal axolemma (Fig 2b). In fact myelinated axon produced two types of APs: classical (saltatorial) AP, known to jump fast from node to node across myelin; and internodal APs induced by the saltatorial one at each node. As a result, two internodal waves originated at each node. They slowly propagated along the axolemma in opposite directions towards the midinternodes where the internodal waves from neighboring nodes collided. When we used 21 internodal segments, our simulations showed no signs of internodal axolemma activation irrespective of internodal sodium channel density (Fig. 2a).

level 1

0 NODE MYSA level 2

-40

Internodal AP

FLUT

-80

STIN

-120

0

1

2 3 Time (ms)

4

5

Fig. 3. Potentials in periaxonal space simulated at different points along the axon. Saltatorial and internodal action potentials are clearly visible. Diffusive processes were taken into account. They lead to substantial difference in the levels of potentials in the beginning and end of the internodal AP mainly because of changes in the potassium equilibrium potential.

3.4. Effects of internodal AP on saltatorial AP The development of internodal AP in the regions near nodes affected the falling phase of the saltatorial AP spike, making it prolonged at the foot (compare transmembrane potentials in Figs. 2a and b, and intracellular potentials represented in Fig. 4 for 21 and 338 segments). Nevertheless, the saltatorial spike duration defined at the 50% level of the spike amplitude was only slightly affected by the internodal APs (Fig. 4). Following the changes in periaxonal potential produced by the moving internodal waves, currents through myelin occurred. Up to the waves’ collision, they induced currents that were depolarizing for axon. This resulted in DAP whose magnitude and slope were sensitive to myelin resistance (Fig 5a, the upper three pairs of curves). When ion accumulation was not included in the model, the DAP was significantly attenuated (Fig. 5a, lowest pair of curves). Myelin capacity practically did not affect the DAP level. During the waves’ collision, the balance between the sodium and potassium currents along the axolemma was

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

40

-70

455

Saltatorial spike truncated

Intracellular potential (mV)

Intracellular potential (mV)

Decay during collision 21 segments

0

338 segments no accumulation 338 segments accumulation

-40

-72

of the internodal waves

-74

R my =0.5R R my =R

-76

R my =2R no accumulation (R my =R)

-80

0

2

4 Time (ms)

6

8

(a)

-78

0

-70

disturbed. Transmembrane currents became mostly potassium along the disappearing internodal waves, i.e. outcoming. This resulted in an abrupt decay of intracellular potential (Fig. 5a). The decay occurred at about 5.5 ms after initiation of the saltatorial spike under the parameters used. The magnitude of decay was sensitive to the myelin capacity (data not shown). After the collision finishing, the DAP decay in our model depended mainly on leakage-driven potassium clearance from the periaxonal space. The DAP and its decay were sensitive to density of the internodal channels. The decay was abrupt (Fig. 5a) when distribution of potassium channel density (Table 2) corresponded to that accepted by Halter and Clark (1991). However, Roper and Schwartz (1987) reported significantly lower values for potassium and leakage channel densities in the internodes. When we reduced the values for potassium and leakage channel densities 10 times in the STIN region and set the sodium channel density there to be 0.5% of the nodal value, the decay was gradual (Fig. 5b). The velocity of propagation of the internodal waves depended on the density of the internodal channels as well. This could be seen from the different position of the DAP decay in Figs. 5a and b.

4. Discussion 4.1. Active internodal waves The existence of narrow periaxonal space is an important prerequisite for the generation of active internodal AP. When the periaxonal width is 400 nm

10

15

20

Saltatorial spike truncated Decay during collision

Intracellular potential (mV)

Fig. 4. Effect of the presence of internodal wave (curves for 338 internodal segments) and diffusive processes on the intracellular potentials shown for the 11th node. Internodal sodium channel density 1 was 30 of the nodal one. The curve for the 21 segments corresponds to the case without internodal wave. The arrowhead marks the start of the internodal waves’ collision at the midinternode.

5

-72

-74

-76

-78 (b)

of the internodal waves

0

5

10

15

20

Time (ms)

Fig. 5. Manifestation of depolarizing afterpotentials in intracellular potentials. (a) Under various myelin resistances (Rmy) without (lowest pair of curves) and with (upper three pairs of curves) taking into account the diffusive processes. The density of sodium channels along 1 the entire internode was 30 (3.33%) of the nodal value. The distributions of potassium channel density and of leakage current corresponded to those given in Table 2. ‘R’ denotes the myelin resistance used by Halter and Clark (1991). (b) Under changed parameters in the STIN region: sodium channel density was equal to 1 200 (0.5%) of the nodal value; potassium and leakage channel densities were reduced 10 times. Curves are given in pairs corresponding to node (interrupted curves) and point of internode (solid lines) 600 mm away from the node. Arrows indicate a small deflection corresponding to the direct manifestation of the internodal AP spike in the internodal intracellular detections.

or more, our model predicts no internodal AP. This is consistent with the lack of AP propagation in demyelinated axons. In normal axons the actual periaxonal width varies from a few nm up to 30 nm (Berthold and Rydmark, 1995). In our model, the saltatorial AP induces passive charging of internodal axolemma only until initiation of internodal AP. Nevertheless, significant DAP is predicted. Thus, internodal AP provides an active mechanism for DAP creation, which is alternative to the passive one. Internodal AP can assure also internodal potassium conductance activation found by

ARTICLE IN PRESS 456

A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

David et al. (1995), as well as significant internodal potassium accumulation after post-tetanic activation of an axon (Baker, 2000; Bostock and Bergmans, 1994). An indirect support for an active internodal process is the presence of Na–K pump along the internodal axolemma in normal mammalian myelinated fibers (Mata et al., 1991). This suggests extensive movements of sodium and potassium ions across the axolemma that are normally associated with active processes. On the other hand, studying internodal sodium channels in demyelinated rabbit axons, Chiu and Schwarz (1987) concluded that it might be possible for an internodal action potential to be generated. It was also suggested (Hines and Shrager, 1991; Shrager, 1995) that the safety factor of demyelinated internodal membrane is very close to 1.

channels (18,000/mm2 instead of 2000/mm2) was used to compensate for the unrealistic dynamics of sodium channels. This means that the actual dynamics will require much lower value for both nodal and internodal sodium channel density for existence of saltatorial and internodal active waves. Therefore, internodal channels, whose densities are less than 25/mm2, are not to be ignored. The absence of pumps in the model, where all ions may diffuse and accumulate, led to equilibrium when potassium, leakage and resting potentials coincided. Therefore, no potassium current could hyperpolarize axon and induce afterhyperpolarization of the AP in our model.

4.2. Limitations of the model

Like in experimental recordings, in our simulations, DAPs can initially increase, be almost unchanged, or decrease (Fig. 5a) depending on parameter values. Decay of the DAP is often not as abrupt as shown in Fig. 5a, although a rather steep decay could be seen, for example, in experiments by David et al. (1995, their Fig. 9C, control, at the middle of the DAP after the first AP). Reduction of the sodium and potassium channel densities in STIN region slowed down the process of waves’ collision and reduced potassium component of current accompanying the process. This resulted in DAP (Fig. 5b) with a much more gradual decay. Highfrequency attenuation at the microelectrode tip as reported by David et al. (1995) could make the decay during and after collision still more gradual. Further study of the effects of different fiber parameters on internodal AP and DAP will be considered in another paper.

Like other multi-cable models of mammalian myelinated axon (Halter and Clark, 1991; McIntyre et al., 2002; Nygren and Halter, 1999; Stephanova, 2001; Stephanova and Bostock, 1995; Zhou and Chiu, 2001), our model has difficulties with the sodium channel dynamics. These channels comprise a diverse family of related proteins encoded by different genes, and have different physiological properties, functions, and highly dynamic activation (Waxman et al., 2000). Ten different a subunits have been identified as Nav1.x family of channels, which are similar functionally and differ in their activation kinetics (Salzer, 2003). The models usually include only one type of the fast voltage-gated sodium channels with certain dynamics. To ensure realistic properties, particularly propagation velocity of the saltatorial AP, all of the models modify the sodium channel dynamics. It is done either directly (McIntyre et al., 2002)—by shifting characteristics of sodium channel along voltage axis in hyperpolarizing direction, or by increasing nodal sodium channel density (Halter and Clark, 1991; Nygren and Halter, 1999; Stephanova, 2001; Stephanova and Bostock, 1995; Zhou and Chiu, 2001). Following Halter and Clark (1991), we also used the increased nodal sodium channel density. Was the density of the internodal sodium channels (sufficient for internodal wave propagation) comparable to the literature data in this case? As percent of the nodal density, 0.5% at initial and 0.3% at reduced values for internodal potassium channels, it was certainly low (the literature data by Burke et al. (2001), Chiu and Schwarz (1987), Shrager (1995), Vabnick et al. (1997), Vabnick and Shrager (1998), reports 2–6% of the nodal value). The lower value, expressed in number of channels per mm2, corresponded approximately to 50/mm2, i.e. it was twice higher than the value reported in literature. On the other hand, the increased nodal density of these

4.3. Shape of DAP

4.4. Difficulties in experimental verification of hypotheses Many recordings one would like to make experimentally are difficult or impossible. Measurements of periaxonal potentials (David et al., 1992, 1993; Funch and Faber, 1984) do not show direct signs of internodal activation. The amplitude of the periaxonal AP detected by the authors was approximately twice lower than that of AP detected intra-axonally. If this had been the actual case, the saltatorial AP propagation would have been impossible, as the great portion of the saltatorial AP would have existed over the internodal axolemma. Since there were experimental problems with registering amplitude of a strong saltatorial AP, it is not surprising that a very weak internodal wave was invisible in such detections. This makes also understandable why Barrett and Barrett (1982) suggested and David et al. (1992, 1993, 1995) supported a passive nature of DAP. Up to now, the active processes in axons were associated with saltatorial AP that could be registered in intracellular or

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

peri-internodal (David et al., 1992, 1993) detections. As internodal transmembrane potentials (Fig. 2b) are mainly pronounced in periaxonal space (Fig. 3) and are practically invisible in intracellular recordings (Fig. 5, arrows), it was natural to consider their effects as produced by a passive process. Beside our mechanism for DAP creation, one has to bear in mind that McIntyre et al. (2002) proposed another active process that could contribute to DAP creation. The authors included additional persistent sodium conductance in the nodal membrane. The dynamics of persistent sodium channels is still unknown. Therefore, now it is difficult to predict signs of changes in DAP in accordance with their hypothesis. Looking for onset of the internodal wave’s collision and repolarization of the internodal axonal membrane associated with the collision, could help detecting our active mechanism of DAP creation. 4.5. Advantages of the proposed numerical approach Our results show that ensuring the appropriate physiological parameters is not sufficient to simulate internodal activation. Zhou et al. (1999) and Zhou and Chiu (2001) used even higher internodal sodium channel 1 density (4%, i.e. 25 of the nodal one) but internodes were not activated in their simulations. The problem was that the spatial segmentation of the modeled space has to be sufficiently fine. To assure internodal AP propagation, the spatial segments have to be at most 20 mm long. As the periaxonal space is very narrow, quantitative modeling of internodal activation requires also calculation of diffusive processes at least in the periaxonal space. This means adding diffusive cables to the model, one per each ion at least in this space. The result is a huge total number (N) of segments representing the myelinated fiber. The currently accepted numerical methods (Halter and Clark, 1991; Hines and Carnevale, 1997; Nygren and Halter, 1999) are non-iterative and would require N3 calculations for cases like ours. Reducing the matrix of general type to tridiagonal matrixes, our iterative numerical method thus requires N calculations (Press et al., 1992). With non-iterative methods, previous authors had to use relatively short time steps to achieve the needed accuracy of the solutions, namely: 0.01 ms (Stephanova, 2001), 0.1 ms (Halter and Clark, 1991), 0.2 ms (Nygren and Halter, 1999), 1 ms (McIntyre et al., 2002; Richardson et al., 2000). This additionally increases the computational time and reduces the possibility for a more detailed simulation of ionic processes in myelinated fibers. The time step in our calculations is considerably longer (10 ms). This compensates to a great extent the number of iterations, which vary from 4 to 50 depending on the presence of fast changes in the membrane potential.

457

5. Conclusions The results of the present modeling study show that activation of low-density internodal sodium channels becomes physiologically natural in restricted volumes. Internodal waves could allow activation of a wide variety of internodal voltage-gated channels and pumps. The channels whose existence is neglected may play an active role. Affecting DAP, they could regulate the fiber excitability. Active processes in internodal axolemma make periaxonal space with its large surface area an important site for interaction with Schwann cell and for voltage-driven metabolic processes.

Appendix A. Method development To derive equations that describe excitation of a myelinated axon, we use the continuity equation, which says that change in the quantity of some substance in a given volume for a given time is equal to the total flow of that substance through the surface of that volume. When this substance is electric charge, we obtain equations for electric potentials, and when the substance is some specific ion, we obtain equations for concentrations of that ion. The myelinated axon is represented through 3 cables (Fig. 1) corresponding to intracellular (Cable 1), periaxonal (Cable 2) and extracellular (Cable 3 treated as grounded) spaces. To perform calculations, we apply the continuity equation to a finite-difference approximation of the modeled space. To achieve numerical stability (Press et al., 1992), we have to consider all the currents in the future moment in time (except for derivatives in time), i.e. write implicit (same as backward Euler) system of equations for our Initial Value problem. The procedure for any electrical and diffusive cable is identical; therefore, we show it only for one electrical Cable 2. For any segment
i, the change of electric charge dQ2i =dt can be considered in two ways. On one hand, it is equal to capacitive current through myelin and axonal membrane (Eq. (A.1)). On the other hand, through the continuity equation, it is equal (Eq. (A.2)) to the sum of resistive and diffusive currents from neighboring segments of Cable 2 (I2), as well as currents through axonal membrane (I12) and myelin (I23). dQ2i DðV 2i
V 1i Þ DðV 2i
V 3i Þ þ C23i ,  C12i Dt Dt dt (A.1) dQ2i ¼ I2iþ1
I2i þ I12i þ I23i , dt

(A.2)

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

458

where Dtdt is time step t t
V 1tþ1 DðV 2i
V 1i Þ ¼ ðV 2tþ1 i i Þ
ðV 2i
V 1i Þ, t t
V 3tþ1 DðV 2i
V 3i Þ ¼ ðV 2tþ1 i i Þ
ðV 2i
V 3i Þ.

For definitions of the specific currents, see Appendix B. When we write all the currents from Eq. (A.1) and Eq. (A.2) in a future moment (t+1) in time we will get a system of equations (Eq. (A.3)) for the potentials in Cable 2 in future moment in time. C12i C23i ðV 2tþ1 ðV 2tþ1
V 1tþ1
V 3tþ1 i i Þþ i i Þ Dt Dt tþ1 ðV 2tþ1 ðV 2tþ1
V 2tþ1 iþ1
V 2i Þ i i
1 Þ þ
tþ1 tþ1 RV 2iþ1 RV 2i ðV 2tþ1
V 3tþ1 i i Þ þ
Id2tþ1
Id23tþ1 i i tþ1 RV 23i C12i C23i ðV 2ti
V 1ti Þ þ ðV 2ti
V 3ti Þ þ I12i . ¼ Dt Dt ðA:3Þ The left side of Eq. (A.3) contains terms proportional to potentials in the future moment in time. The right side contains the terms proportional to potentials in the present moment in time, and all nonlinear ionic currents through the axonal membrane. Halter and Clark (1991) obtained similar equation (diffusive terms excluded). To get a complete set of equations, it is necessary to do the same for all cables (for Cable 1 and also for diffusive cables, if present). To solve equations, let us pay attention to the interactions between cables. One could look at the interactions as previous authors (Halter and Clark, 1991; Nygren and Halter, 1999) did. They treated all cables as one system of equations represented by one matrix. This means that if V1, V2, and the diffusive cables (when present), form the corresponding parts of the solution vector, then interactions between segments of the same cable in the left (linear) part of that system, form a tridiagonal subsystem. However, interactions between segments of tþ1 the two cables (i.e. terms proportional to V 1tþ1 i ; V 3i ; tþ1 tþ1 also Id2i and Id23i for Cable 2) form off-tridiagonal elements. As Halter and Clark (1991) used a noniterative approach, they had to use a general matrix inversion method to solve this system. We have noticed that only the tridiagonal part is responsible for the numerical stability. Therefore, we have moved the off-tridiagonal elements to the right (nonlinear) part of the system. This manipulation is correct, i.e. it does not change the solutions of the initial system of equations, because we calculate the currents in the right side of the equations iteratively. Thus, when the solutions converge, all the currents in the right side of the equation—linear or nonlinear, are correctly included into the equations at the correct moment in time and therefore, are correctly reflected in the

solutions. Thus, the matrix of general type (used by all previous authors) disintegrates into tridiagonal matrixes whose number is equal to that of the cables (electrical and diffusive). Therefore, we have the following set of equations for a two-cable model: tþ1 þ V 1tþ1 V 1tþ1 iþ1 PV 1iþ1
FV 1i V 1i i
1 PV 1i ¼ DV 1i ,

(A.4) tþ1 V 2tþ1 þ V 2tþ1 iþ1 PV 2iþ1
FV 2i V 2i i
1 PV 2i ¼ DV 2i ,

where DV1i and DV2i are the nonlinear parts. For the coefficients we have FV 1i ¼ ðk112 C12i =dtÞ þ PV 1iþ1 þ PV 1i tþ1

(A.5) tþ1

DV 1i ¼ ðC12i ðV 10 i ð1
k112 Þ
V 1ti
V 20 i þ V 2ti Þ=dt þ I12i
Id10 i
IstimÞ FV 2i ¼ ððC12i þ C23i Þ=dt þ PV 23i Þ þ PV 2iþ1 þ PV 2i tþ1

DV 2i ¼ ððC12i ð
V 2ti
V 10 i tþ1

ð
V 2ti
V 30 i

þ V 1ti Þ þ C23i

þ V 3ti ÞÞ=dt tþ1


I12i
Id2_230 i
PV 23i V 30 i Þ k112 is an important optimization coefficient. If k112 ¼ 1, then Eq. (A.5) directly follows from Eq. (A.3) but the system converges very slowly, with its speed determined by the relatively large capacity of axonal membrane. We have found that the system converges fast enough, if the capacity associated with the axonal membrane in Cable 1 is decreased in the left (tridiagonal) portion of the system. To preserve the equation, we transfer a large part of the capacity current to be added iteratively instead of directly, to our equation (what k112 actually does). This is because only the direct part of the capacity current (namely, that included in FV1 of Eq. (A.5)) affects convergence. Once again, when the equations converge, the capacity current is represented correctly. It is necessary to perform the above with care, because moving too much capacity current to the nonlinear part of the equations could cause system instability. Eqs. (A.3–A.4) form a fully implicit, first-order in time system. To obtain a numerically stable system that is second-order in time, one has to transform the equations into an explicit form (i.e. to write them for the present moment in time) and then to sum up the explicit equations and the corresponding implicit ones (Press et al., 1992). Both types of systems are used for modeling of a myelinated axon (Halter and Clark, 1991; Nygren and Halter, 1995). We used the second-order like in the model by Halter and Clark (1991).

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

Appendix B. Full system of equations

ðVX ti
VX ti
1 Þ , RLeX ti X ILdX ti , IdX ti ¼

ILeX ti ¼

In our calculations, we used the following second order in time set of equations for potentials (L ¼ V) and ions (L ¼ Na, K, Cl, Ca) included in the model (for specific notations see Table 1): tþ1 tþ1 tþ1 LX tþ1 þ LX tþ1 ¼ DLX i iþ1 PLX iþ1
FLX i LX i i
1 PLX i

L¼Na;K;Cl;Ca

X

IX ti ¼

ILeX ti þ IdX ti ,

L¼Na;K;Cl;Ca

I12ti ¼ ðIna12ti þ Ikf 12ti þ Iks12ti þ Ilk12ti ÞðS12i Þ,

where L ¼ V; Na; K; Cl; Ca; X ¼ 1; 2. For the coefficients FLX, DLX, in the case of L ¼ V we have Eq. (B.1) FV 1i ¼ ð2k112 C12i =dtÞ þ

PV 1tþ1 iþ1

þ

PV 1tþ1 i ,

FV 2i ¼ ð2ðC12i þ C23i Þ=dt þ PV 23tþ1 i Þ tþ1 þ PV 2tþ1 þ PV 2iþ1 i ,

DV 1i ¼

459

where Ina12ti ; Ikf 12ti ; Iks12ti and Ilk12ti are sodium, fast and slow potassium, and leakage current densities. For resistivities expressed through the ions’ concentrations, we have

þ V 3ti ÞÞ=dt
I12tþ1 i tþ1

For the coefficients FLX and DLX in the case of L ¼ Na, K, Cl, Ca, we have Eq. (B.2) tþ1 FL1i ¼ ð2zL F Vol1i =dtÞ þ PL1tþ1 iþ1 þ PL1i ,

LX ti Þ

 dxX i dxX i
1 þ =2, AX i uL zL F LX ti AX i
1 uL zL F LX ti
1



PLe23ti ¼

S23i uL23i zL F 2 L2ti L3ti 1 t t  dx23i ðL2i þ L3i Þ RLe23ti

P2Le23ti ¼

tþ1 þ PL2tþ1 iþ1 þ PL2i ,

DL1i ¼ ð2zL F Vol1i ð
L1ti Þ=dt þ IL12tþ1 þ IL12ti i tþ1
ILd1tiþ1 þ ILd1ti
ILe1iþ1
ILe1tiþ1

þ

L¼Na;K;Cl;Ca uL zL F

S23i uL23i zL F 2 L3ti ðV 2ti
V 3ti Þ dx23i ðL2ti þ L3ti Þ ðV 2ti
V 3ti Þ  PLe23ti , L2ti

FL2i ¼ ð2zL F Vol2i =dt þ PL23tþ1 þ P2Le23tþ1 i i Þ

þ

dxX i

 dxX i dxX i
1 þ =2, AX i uL RT AX i
1 uL RT   dx23i dx23i t þ RLe23i ¼ =2, S23i uL23i zL FL2ti S23i uL23i zL FL3ti RLX ti ¼


PV 23ti

ILe1ti Þ,

P



RLeX ti ¼


Id230 i
Id23ti
PV 23tþ1 V 30 i i t t ðV 3i
V 2i ÞÞ.

ILe1tþ1 i

AX i ð

ðB:1Þ

t tþ1
I12ti
Id2tþ1 þ Id2ti iþ1
Id2iþ1 þ Id2i tþ1

ðL3ti
L2ti Þ , RL23ti L¼Na;K;Cl;Ca

! dxX i
1 P þ =2, AX i
1 ð L¼Na;K;Cl;Ca uL zL F LX ti
1 Þ

þ V 1ti Þ DV 2i ¼ ð2ðC12i ð
V 2ti
V 10 tþ1 i tþ1

X

RV X ti ¼

tþ1 tþ1 ð2C12i ðV 10 i ð1
k112 Þ
V 1ti
V 20 i þ V 2ti Þ=dt þ I12tþ1 þ I12ti
Id1tþ1 i iþ1 t
Id1tiþ1 þ Id1tþ1 þ Id1
IstimÞ, i i

þ C23i ð
V 2ti
V 30 i

Id23ti ¼

X

PV 23ti ¼ ðB:2Þ

DL2i ¼ ð2zL F Vol2i ð
L2ti Þ=dt
IL12tþ1
IL12ti i

L¼Na;K;Cl;Ca

RL23ti ¼



tþ1
ILd2tiþ1 þ ILd2ti
ILe2iþ1
ILe2tiþ1

PLe23ti 

1 , RV 23ti

 dx23i dx23i þ =2 S23i uL23i RT S23i uL23i RT

þ ILe2tþ1 þ ILe2ti
PL23tþ1 L3tþ1 i i i

ðL ¼ Na; K; Cl; CaÞ.


PL23ti ðL3ti
L2ti Þ þ P2Le23ti L2ti Þ.

We calculate cross-sectional and surface areas AXi, SXi, and volumes VolXi as if segments were cylindrical in shape. Following Halter and Clark (1991), we used different radii for their calculation to follow the complex axon geometry:

For the corresponding diffusive (d), electrical (e), total and ionic currents we have ILdX ti ¼

ðLX ti
LX ti
1 Þ , RLX ti

A1i ¼ pr2i ,

ARTICLE IN PRESS 460

A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

Table 1 The structure of notations ðMXb0 ti Þ; possible values and meanings Notation

Possible Values

0

t i X M ¼ NLf

RL

t or t+1 i
1, i, or i+1 1, 2, 3 or 12, 23 N ¼ I, R, P, G, u L ¼ V, L ¼ Na, K, Cl, Ca L ¼ Vol L ¼ na, kf, ks, lk f ¼ d or e L¼V L ¼ Na, K, Cl, Ca

A S dx a r da dr uL zL nl mt

Meaning The parameter is obtained iteratively present or future moment in time previous, current or next fiber segment number of the cable or of the transition between neighboring cables current (I), resistivity (R), conductivity (P), channel dynamics (G), mobility (u) voltage concentration of the corresponding ion volume channel type (sodium, fast potassium, slow potassium, leakage) d—diffusive; e—electrical electrical resistivity specific ion diffusive resistivity segment cross-sectional area membrane or myelin area segment length radius of axon (represents membrane area and periaxonal volume) radius of axon (represents axonal cross-section area) periaxonal space width (represents periaxonal volume) periaxonal space width (represents periaxonal conductivity) mobility of the ion L electric charge of the ion L number of myelin lamellae myelin lamella thickness

A2i ¼ 2pai dri , S12i ¼ 2pai dxi , S23i ¼ 2pðai þ dri þ nl mt=2Þdxi ;

i 2 internode,

Vol1i ¼ pr2i dxi , Vol2i ¼ 2pai dai dxi , where ai and ri represent internal radius of the axon and differ only in FLUT that has non-cylindrical geometry (Fig. 1). dai and dri represent the width of periaxonal space; they differ only in MYSA where the periaxonal volume corresponds to da ¼ 40 nm (Berthold and Rydmark, 1983), while periaxonal conductivity is computed with dr ¼ 1 nm (Halter and Clark, 1991). Values for the radii at the ends of the main sections of axon are represented in Table 2. We assume that they change linearly between those values. C12i ¼ c S12i C23i ¼ ðc=2ÞS23i =nl;

i 2 internode:

Note that single lamella is formed from two membranes. For the nodal gap capacity, we took the capacity of a nodal gap Schwann cell membrane. For calculating the nodal gap conductivity we divide the extracellular nodal gap volume by nodal gap height to obtain nodal gap cross-sectional area.

For ionic mobilities, we took their values for free aqueous solutions, everywhere (including nodal gap) but inside myelin. In myelin, we scaled all mobilities (Dobos, 1975; Koshkin and Shirkevich, 1988) by (2 10
5)/nl and achieved conductivity of about 1.566 mS/cm2 per myelin lamella, which is close to the 1.56 mS/cm2 used by Halter and Clark (1991). For the axonal membrane current densities we have ðzL F Þ2 ðV 1ti
V 2ti Þ RT ðL2ti
L1ti expðzL F ðV 1ti
V 2ti Þ=ðRTÞÞÞ , 1
expðzL F ðV 1ti
V 2ti Þ=ðRTÞÞ

Il12ti ¼ uL12i Gl12ti

Ilk12ti ¼ glk12ti ðV 1ti
V 2ti
Vlk12ti Þ, where l ¼ na, kf, ks; Gl12 represents channel dynamics; L represents corresponding ion. Gna12ti ¼ ðmti Þ3 hti , Gkf 12ti ¼ ðnf ti Þ4 , Gks12ti ¼ ðnsti Þ4 . The general equation for gating variables is dy ¼ ay ð1
yÞ þ by y; dt

where y ¼ m; h; nf ; ns.

ARTICLE IN PRESS A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

461

Table 2 Parameters used in simulations Parameter d

Number of segments Length (mm) r (mm) a (mm) dr (nm) da (nm) Na3 (mM) Na2 (mM) Na1 (mM) K3 (mM) K2 (mM) K1 (mM) Pna12 (cm/s)a Pkf12 (cm/s) Pks12 (cm/s) glk12 (mS/cm2) c (mF/cm2) gm (mS/cm2) nl mt (mm) Vlk12 (mV) uNaX (m2/(s V))b uKX (m2/(s V))b uClX (m2/(s V))b uCaX (m2/(s V))c k112d dt (ms)d T(1C)

Node

MYSA

FLUT

FLUT

STIN

1 1.5 2.5 2.5 30 30 154 154 8.71 5.9 5.9 155 33.9 10
3 0.375 10
3 0.0937 10
3 156 1 — — —
77.9625

4 4 2.2–2.4 2.2–2.4 1 40 154 154 8,71 5.9 5.9 155 1.13 10
3 0.0625 10
3 0.0312 10
3 1 1 1.56 175 0.018
76.48

5 5 2.4–5.1 2.4–9 3 3 154 154 8,71 5.9 5.9 155 1.13 10
3 0.0625 10
3 0.0312 10
3 1 1 1.56 175 0.018
76.48 6.21 10
8 9.23 10
8
9.38 10
8 7.31 10
8 0.02 10 37

70 70 5.1–6.25 9–6.25 3 3 154 154 8.71 5.9 5.9 155 1.13 10
3 0.0625 10
3 0.0312 10
3 1 1 1.56 175 0.018
76.48

180 1800 6.25 6.25 4 4 154 154 8,71 5.9 5.9 155 1.13 10
3 0.0625 10
3 0.0312 10
3 1 1 1.56 175 0.018
76.48

1 of the nodal value, bKoshkin and Shirkevich All values in Table 2 are as in Halter and Clark (1991), except: aInternodal sodium channel density is 30 c d (1988), Dobos (1975), Introduced by us.

Rate constants ðT
T0Þ=10

am ; anf ; ans ¼ Q10

ah ; bm ; bnf ; bns ¼ ðT
T0Þ=10

bh ¼ Q10

Table 3 Parameters used for calculation of rate constants

AðV
BÞ , 1
expððB
V Þ=CÞ

ðT
T0Þ=10 Q10

AðB
V Þ , 1
expððV
BÞ=CÞ

A , 1 þ expððB
V Þ=CÞ

where V ¼ V1
V2
Vrest; Vrest ¼
78 mV (Schwarz and Eikhof, 1987), constants A, B and C (Halter and Clark, 1991; Ilyin et al., 1980; Schwarz and Eikhof, 1987) are represented in Table 3. The general form of the gating variable approximation is ytþ1 ¼

yt þ fay ðV tþ1 Þ þ ½ay ðV t Þ
ðay ðV t Þ þ by ðV t ÞÞyt g dt=2 1 þ ðay ðV tþ1 Þ þ by ðV tþ1 ÞÞ dt=2

,

where y ¼ m; h; nf ; ns: We used a second-order in time method for integrating gating variables as described by Cooley and Dodge (1966),

Parameter

Q10

T0 (1C)

A (m s
1)

B (mV)

C (mV)

am bm ah bh anf bnf ans bns

2.2 2.2 2.9 2.9 3.2 2.8 3.2 2.8

20 20 20 20 20 20 20 20

0.49 1.04 0.09 3.7 0.0078 0.0214 0.00092 0.0042

25.41 21
27.74 56
20.17 17.88
25.75 0.24

6.06 9.41 9.06 12.5 38.85 0.91 15.73 8.54

because the method described by Halter and Clark (1991) is based on linear approximation of nonlinear currents and, therefore, requires too small time steps. Our criterion for sufficient accuracy of potential calculation was as follows. We chose the difference between estimates of potentials and concentrations in consecutive iterations to be less than 10
5 mV and 10
6 mM, respectively. One could refer to Press et al. (1992) or Cooley and Dodge (1966) for the entirely technical issue on solving tridiagonal equations.

ARTICLE IN PRESS 462

A.G. Dimitrov / Journal of Theoretical Biology 235 (2005) 451–462

References Baker, M.D., 2000. Axonal flip-flops and oscillators. TINS 23, 514–519. Barrett, E.F., Barrett, J.N., 1982. Intracellular recording from vertebrate myelinated axons: mechanism of the depolarizing afterpotential. J. Physiol. 323, 117–144. Berthold, C.H., Rydmark, M., 1983. Anatomy of the paranode–node–paranode region in the cat. Experientia 39, 964–976. Berthold, C.H., Rydmark, M., 1995. Morphology of normal peripheral axons. In: Waxman, S.G., Kocsis, J.D., Stys, P.K. (Eds.), Axon. Oxford University Press, New York, pp. 13–48. Blight, A.R., 1985. Computer simulation of action potentials and afterpotentials in mammalian myelinated axons: the case for a lower resistance myelin sheath. Neuroscience 15, 13–31. Blight, A.R., Someya, S., 1985. Depolarizing afterpotentials in myelinated axons of mammalian spinal cord. Neuroscience 15, 1–12. Bostock, H., Bergmans, J., 1994. Post-tetanic excitability changes and ectopic discharges in a human motor axon. Brain 117, 913–928. Bostock, H., Sears, T.A., 1978. The internodal axon membrane: electrical excitability and continuous conduction in segmental demyelination. J. Physiol. 280, 273–301. Burke, D., Kiernan, M.C., Bostock, H., 2001. Excitability of human axons. Clin. Neurophysiol. 112, 1575–1585. Chiu, S.Y., Schwarz, W., 1987. Sodium and potassium currents in acutely demyelinated internodes of rabbit sciatic nerves. J. Physiol. 391, 631–649. Cooley, J.W., Dodge, F.A., 1966. Digital computer solutions for excitation and propagation of the nerve impulse. Biophys. J. 6, 583–599. David, G., Barrett, J.N., Barrett, E.F., 1992. Evidence that action potentials activate an internodal potassium conductance in lizard myelinated axons. J. Physiol. 445, 277–301. David, G., Barrett, J.N., Barrett, E.F., 1993. Activation of internodal potassium conductance in rat myelinated axons. J. Physiol. 472, 177–202. David, G., Modney, B., Scappaticci, K.A., Barrett, J.N., Barrett, E.F., 1995. Electrical and morphological factors influencing the depolarizing after-potential in rat and lizard myelinated axons. J. Physiol. 489 (Pt 1), 141–157. Dobos, D., 1975. Electrochemical Data (a Handbook for Electrochemists, Industry and Universities). Akademiai Kiado, Budapest. Funch, P.G., Faber, D.S., 1984. Measurement of myelin sheath resistances: implications for axonal conduction and pathophysiology. Science 225, 538–540. Halter, J.A., Clark, J.W., 1991. A distributed-parameter model of the myelinated nerve fiber. J. Theor. Biol. 148, 345–385; Errata. J. Theor. Biol. 152, 569–571. Hines, M.L., Carnevale, N.T., 1997. The NEURON simulation environment. Neural Comput. 9, 1179–1209. Hines, M.L., Shrager, P., 1991. A computational test of the requirements for conduction in demyelinated axons. Restor. Neurol. Neurosci. 3, 81–93. Ilyin, V.I., Katina, I.E., Lonskii, A.V., Makovski, V.S., Polishchuk, E.V., 1980. The Cole–Moore effect in nodal membrane of the frog

Rana Ridibunda: evidence for fast and slow potassium channels. J. Membrane Biol. 57, 179–193. Koshkin, N.I., Shirkevich, M.G., 1988. Handbook on Elementary Physics. Nauka, Moscow, main editorial board for physical and mathematical literature (in Russian). Mata, M., Fink, D.J., Ernst, S.A., Siegel, G.J., 1991. Immunocytochemical demonstration of Na+, K+-ATPase in internodal axolemma of myelinated fibres of rat sciatic and optic nerves. J. Neurochem. 57, 184–192. McIntyre, C.C., Richardson, A.G., Grill, W.M., 2002. Modeling the excitability of mammalian nerve fibers: influence of afterpotentials on the recovery cycle. J. Neurophysiol. 87, 995–1006. Nygren, A., Halter, J.A., 1999. A general approach to modeling conduction and concentration dynamics in excitable cells of concentric cylindrical geometry. J. Theor. Biol. 199, 329–358. Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Numerical Recipes in C: The Art of Scientific Computing, second ed. Cambridge University Press, Cambridge. Richardson, A.G., McIntyre, C.C., Grill, W.M., 2000. Modelling the effects of electric fields on nerve fibres: influence of the myelin sheath. Med. Biol.Eng. Comp. 38, 438–446. Roper, J., Schwarz, J.R., 1989. Heterogeneous distribution of fast and slow potassium channels in myelinated rat nerve fibers. J. Physiol . 416, 93–110. Salzer, J.L., 2003. Polarized domains of myelinated axons. Neuron 40, 297–318. Shrager, P., 1995. Action potential conduction recorded optically in normal, demyelinated and remyelinated axons. In: Waxman, S.G., Kocsis, J.D., Stys, P.K. (Eds.), Axon. Oxford University Press, New York, pp. 341–352. Schwarz, J.R., Eikhof, G., 1987. Na currents and action potentials in rat myelinated nerve fibres at 20 and 37 1C. Pflugers Arch. 409, 569–577. Stephanova, D.I., 2001. Myelin as longitudinal conductor: a multilayered model of the myelinated human motor nerve fibre. Biol.Cybern. 84, 301–308. Stephanova, D.I., Bostock, H., 1995. A distributed-parameter model of the myelinated human motor nerve fibre: temporal and spatial distributions of action potentials and ionic currents. Biol. Cybern. 73, 275–280. Vabnick, I., Shrager, P., 1998. Ion channel redistribution and function during development of the myelinated axon. J. Neurobiol. 37, 80–96. Vabnick, I., Messing, A., Chiu, S.Y., Levinson, S.R., Schachner, M., Roder, J., Li, C., Novakovic, S., Shrager, P., 1997. Sodium channel distribution in axons of hypomyelinated and MAG null mutant mice. J. Neurosci. Res. 50, 321–336. Waxman, S.D., Dib-Hajj, S., Cummins, T.R., Black, J.A., 2000. Sodium channels and their gens: dynamic expression in the normal nervous system, dysregulation in disease states. Brain Res. 886, 5–14. Zhou, L., Chiu, S.Y., 2001. Computer model for action potential propagation through branch point in myelinated nerves. J. Neurophysiol. 85, 197–210. Zhou, L., Messing, A., Chiu, S.Y., 1999. Determinants of excitability at transition zones in Kv1.1-deficient myelinated nerves. J. Neurosci. 19, 5768–5781.