The myelinated nerve fibre is formed by an axon and Schwann cells or oligodendrocytes that sheath the axon by winding around it in tight myelin layers. Repetitive stimulation of a fibre is known to result in accumulation of extracellular potassium ions, especially between the axon and the myelin. Uptake of potassium leads to Schwann cell swelling and myelin restructuring that impacts the electrical properties of the myelin. In order to further understand the dynamic interaction that takes place between the myelin and the axon, we have modelled submyelin potassium accumulation and related changes in myelin resistance during prolonged high-frequency stimulation. We predict that potassium-mediated decrease in myelin resistance leads to a functional excitation block with various patterns of altered spike trains. The patterns are found to depend on stimulation frequency and amplitude and to range from no block (less than 100 Hz) to a complete block (greater than 500 Hz). The transitional patterns include intermittent periodic block with interleaved spiking and non-spiking intervals of different relative duration as well as an unstable regime with chaotic switching between the spiking and non-spiking states. Intermittent conduction blocks are accompanied by oscillations of extracellular potassium. The mechanism of conductance block based on myelin restructuring complements the already known and modelled block via hyperpolarization mediated by the axonal sodium pump and potassium depolarization.
The myelinated nerve fibre has a complex and intricate structure as required for fast and reliable conduction of spike trains both in the peripheral and in the central nervous system (CNS). This structure is a product of two interacting cell types: the axon and the glial cell. In the CNS, the myelinating cell is the oligodendrocyte and in the peripheral nervous system (PNS) it is the Schwann cell. The structure of central and peripheral myelinated fibres is similar, and below we will address mainly the peripheral nerve, keeping in mind that similar mechanisms work for the CNS as well. Failures in dynamic interaction between the sheathing glial cell and the axon lead to conduction block, paraesthesia and other pathologies.
In some conditions, such as nerve ischaemia, demyelination, early remyelination, disruptions of junctions between the glial cell and astrocytes in the CNS, the sheathing glial cell fails to cope with the uptake of potassium ions, fluxing out through the open K+-channels of the depolarized paranodal axolemma. As a result of this failure, K+ accumulates in the confined extracellular space between the axon and the sheath (periaxonal space), and to a lesser extent around the node of Ranvier and outside of the fibre. The effect of the K+-accumulation can be the following: (i) it shifts the reversal potential for K+-channels and thus depolarizes the membrane, (ii) activates the axonal Na+,K+-pump, which creates hyperpolarizing currents, and (iii) pumping K+ into the Schwann cell leads to its swelling and to alteration of the myelin structure, which in its effect is analogous to focal demyelination, i.e. a reduction of myelin resistance. The effect of potassium depolarization has been modelled e.g. by Bostock et al. [1,2] for the case of ischaemic conditions. The effect of the Na+,K+-pump was studied in human nerves e.g. by Kiernan et al.  and Vagg et al. .
The third aspect of potassium accumulation, i.e. the K+-dependent changes in the myelin properties have not been addressed before. We focused on this subject in the current study and were able to observe the development of K+-dependent excitation block under prolonged repetitive stimulation of the nerve. This block is either intermittent with different patterns of behaviour or persistent depending on model parameters such as stimulation amplitude and frequency and sensitivity of the structural reorganizations of myelin to paranodal submyelin [K+]e .
In this paper, we present a model of a nerve fibre that demonstrates K+-accumulation under the myelin during prolonged repetitive stimulation, and excitation block, caused by K+-dependent structural reorganization in the myelin. We use a space-lumped electrical scheme superimposed on several ionic compartments. Changes in myelin structure mediated by high extracellular potassium are modelled as a concentration-dependent decrease in myelin resistance. The model exhibits oscillatory, chaotic or stationary behaviour with respect to extracellular potassium concentration and myelin resistance and the fibre therefore shows periodic, chaotic or complete excitation block at high frequencies of stimulation.
2. Biological background
2.1. Classical approach
Classical understanding of signal conduction by myelinated nerve fibres summarizes to ‘saltatory conduction’ from one node of Ranvier (short spans of relatively open axonal membrane) to the next. The nodes are separated by long internodal regions, sheathed with myelin, which is produced by the glial cell. The peripheral myelin sheath consists of multiple concentric layers of Schwann cell membranes, as it wraps around the axon multiple times and squeezes the cytoplasm from the resulting roll. Because these layers are densely stapled, the myelin is assumed to have very low electric capacity and very high resistance, thus effectively shielding the internodal axolemma.
Submitting to this angle of view, most of the early works focused on the currents and ion channels of the node, treating the myelin sheath as an ideal isolator and the internodal regions just as longitudinal resistances linking one node to the next. Accompanying the experimental studies, the corresponding mathematical models of myelinated nerve fibres were formulated; first by Frankenhaeuser & Huxley  for the amphibian nerve and later by Chiu et al.  and by Schwarz & Eikhof  for the mammalian nerve. In these models only nodal currents were taken into account. The models provided a basic understanding of excitability and the mechanisms of action potential (AP) generation and propagation in the myelinated fibre. They made it possible to dissect the role of different ion channels on the form of the APs, the refractory period and changes in excitability.
It has gradually become evident, however, that electrical leaks through myelin and in the myelin sheath attachment areas near the nodes (paranodal regions) need to be regarded in order to explain some characteristic features of APs, such as depolarizing afterpotentials (DAP). This led to new electrical schemes for the structure of the models [8,9], where internodal axolemma becomes depolarized in the course of APs. If the classical node-only models were ‘cable models’, these studies (which involved both experiments and modelling) introduced the coaxial (or ‘double cable’) equivalent scheme to model nerve fibres. This scheme turned out to be an important framework as it was recently established that in higher vertebrates, the sodium conductance in the axon is mainly restricted to the node and the potassium conductance is primarily confined to paranodal regions [10,11]. Thus, the depolarizing and hyperpolarizing currents in the myelinated nerves are spatially separated. It has to be noted that voltage-gated potassium channels reside under the paranodal myelin layers and face the limited volume between the apposing axonal and glial membranes.
Generation of APs in the node leads to depolarization of the paranodal and internodal axonal membranes owing to the electrical leaks radially through myelin and longitudinally in the paranodes where the myelin sheath is attached to the axon by so-called paranodal loops of the Schwann cell, filled with cytoplasm. This depolarization, though small, leads to activation of paranodal potassium channels and outflux of potassium from the axon into a confined periaxonal volume. Diffusion from this compartment is limited because the axonal and the glial membranes just near the node come as close to each other as 2 nm and are bound to each other with the help of several special protein families .
2.2. Clearance of potassium
As recently reviewed by Rash , repetitive propagation of APs leads to potassium build-up in the confined periaxonal space. This excess potassium is cleared by Na+,K+-pump back into the axon and in the glial cell. Uptaken by inner layers of the glial cell, the potassium has to be passed to the external layers and then be siphoned away by the syncytium of astrocytes in the CNS or by other mechanisms in the PNS. This trans-myelin transport is possible because of the gap junctions linking the myelin layers. Primary pathways linking inner layers of myelin with outer layers are paranodal cytoplasmic loops flanking the node both in the CNS and in the PNS, and Schmidt–Lanterman incisures in the PNS fibres, which are formed by spiral belts of Schwann cell cytoplasm with adjacent myelin layers connected via gap junctions. The mentioned paranodal loops and Schmidt–Lanterman incisures together with the myelin sheath attachment annuli also constitute the electrical leakage of myelin, regarded by Barrett & Barrett .
Failure at any stage of this potassium clearance mechanism results in a pathology. For example, as reviewed by Rash  mutations in genes coding connexins needed to form gap junctions lead to demyelinating diseases, such as Charcot–Marie–Tooth disease and Pelizaeus–Merzbacher-like diseases. In the first case, the junctions between the myelin layers are disrupted and in the second case the link between oligodendrocytes and astrocytes is broken. Sustained repetitive activity leads to excessive swelling of myelin layers in both cases, because K+ is no more effectively evacuated from the glial cell, and is accumulated especially at inner layers, leading to an osmotic influx of water and swelling. In both cases, the swelling eventually results in segmental demyelination with all layers swollen and separated.
Demyelination, i.e. disruption of the myelin sheath up to a complete loss of it, can also result from trauma, poisoning, malignant bacteria or inflammation. Partial or total lack of normal myelin leads to block of conduction, since sodium channels are located almost exclusively in the nodes, and the stripped internodal regions cannot sustain propagation of pulses, the nodes being unable to recharge the large internodal axolemma. Cases of internodal and paranodal demyelination were addressed in model studies e.g. by Shrager  and by Stephanova & Daskalova [13,14] and where they were able to reproduce clinically observed symptoms of demyelinating neuropathies. However, the role of the myelin sheath in these models was purely passive.
2.3. Activity-dependent changes in myelin structure
It is important to realize that myelin is a dynamic structure and always observed either in one of its steady states or during transition from one state to another [15,16]. There is a balance between the processes of aggregation and splitting of myelin layers, especially in the paranode–node–paranode regions, which together with the Schmidt–Lanterman incisures are the most labile and complex regions of the fibre .
Prolonged repetitive excitation leads to changes in the myelin structure even in normal fibres. Moran & Mateu  have shown that after repetitive excitation at 200 Hz rat fibres become sensitive to 4-aminopyridine (4-AP), which is a blocker for fast voltage-sensitive potassium channels. Emerging sensitivity to 4-AP suggests involvement of the normally sub-myelin paranodal channels into AP generation, which can happen owing to a fall in the Barrett and Barrett resistance, which can be thought of a hallmark of paranodal demyelination.
Wurtz & Ellisman  have shown in frog nerves that repetitive stimulation leads to appearance of vacuoles in the paranodal myelin. Thus, the layers of compact myelin become disturbed and separated, creating cytoplasm-filled vacuoles. The appearance of vacuoles naturally leads to a decrease in myelin resistance, which is in accord with activity-related sensitivity to 4-AP.
The described changes are clearly linked to a rise in extracellular K+. For example, activity-dependent swelling of oligodendrocytes in rat optic fibres observed by MacVicar et al.  was mimicked by artificial increase in [K+]e . This is supported by finding of Jo & Boggs  that K+ leads to main myelin protein (MBP)-mediated vesicle aggregation, resulting in formation of myelin vacuoles. As reported by Gould et al.  shiverer mutant mice, lacking MBP, have more than doubled number of SL-incisures, apparently compensating for the inability of MBP-mediated changes in their structure.
Summarizing, existing evidence leads us to the following picture:
— prolonged repetitive activity leads to potassium accumulation under the myelin, especially in the paranodal region, because it is the place where potassium channels are clustered;
— this extracellular K+ is cleared by the Schwann cell (or the oligodendrocyte). Water follows K+ and the glial cell can swell;
— possibly via MBP-related mechanisms, the structure of the myelin changes: SL-incisures widen and myelin becomes vacuolized. This reduces myelin resistance because compact myelin is now pierced by relatively wide cytoplasm-filled channels;
— K+ taken up by the myelinating cells must be passed outside. Failure at this point leads to excessive swelling and myelin disruption; and
— lower myelin resistance results in more pronounced depolarization of the paranodal axonal membrane, which in turn results in an even higher outflux of potassium and raises the excitation threshold.
Indeed the extent of sub-myelin K+ accumulation and myelin restructuring would depend on nerve fibre type and physiological condition. Clearly, an injury, a genetic disorder or other pathological conditions are likely to make the fibre more susceptible to repetitive stimulation and related changes in myelin structure.
The main idea of this study is to combine this picture into a consistent model framework. In the model discussed below, we address possible effects of coupling accumulation of K+ under the myelin to structural changes leading to an increase in myelin leakage, which in turn results in a raised excitation threshold. We show that depending on the amplitude and frequency of stimulation, the model fibre exhibits different patterns of behaviour such as normal excitation with spike generation and periodic, chaotic or total excitation block.
3. Model description
3.1. Electrical circuit
For simplicity, we use an equivalent circuit for the space-lumped fibre analogous to that of previous authors [1,8,23,24]. For details see figure 1. This circuit layout was first suggested for frog and lizard fibres , but has subsequently been found applicable to mammalian fibres as well .
The general form of the equations to define the nodal Vn and internodal Vi potentials for the given circuit is the following: 3.1and 3.2
Here Cn, Ci and Cm are the electrical capacities of the nodal membrane, internodal membrane and myelin sheath, respectively. Vn and Vi are axonal membrane potentials in the node and internode regions. ∑ Iactcompartment is the sum of all ‘active’ currents through the axonal membrane in the corresponding compartment, including currents through ion channels and through the Na+,K+-pump; Iex is the excitation current. Ril denotes the internodal leakage resistance—through myelin and through the paranodal seal (table 1).
3.2. Currents through ion channels in the axon
Although the main goal of this study is to analyse the implications of the accumulation of K+ under the myelin with its possible influence on the myelin resistance, and not to provide an accurate or species-specific description of excitability of human (or rat) nerve fibres, we still aim at a realistic description of the various ion currents in our model.
All the currents are modelled within the Hodgkin–Huxley formalism, i.e. with the assumption of linear instantaneous voltage–current relationship: 3.3where V is either Vn or Vi where applicable while Vxrest can differ between different concentrational compartments: node, paranode and internode because of different concentrations in them. Vxrest = (RT/zF) ln(Cxe/Cxi) is the Nernst potential for the ionx with extracellular concentration Cxe, intracellular concentration Cxi and charge number z. Here R is universal gas constant and F is Faraday's constant.
Following McIntyre et al.  and Hennings et al. , the node in our model includes two types of sodium currents: transient sodium current (Nat): , and persistent (non-inactivating) sodium current . The transient sodium current is primarily responsible for the rising phase of the AP while persistent sodium current is known to be responsible for part of the DAP  and possibly for breakdown of accommodation . Gating parameters for the Nat current were the same as in Hennings et al.  and Schwarz et al. , for Nap—the same as in Hennings et al. . In the paper by Hennings et al. , the maximal conductivity of the transient sodium current is mS cm−2, whereas it is 3000 mS cm−2 in McIntyre et al. . We use an intermediate value mS cm−2 in our model. Maximal conductivity for the persistent current, is 64.6 mS cm−2 in Hennings et al.  and 10 mS cm−2 in McIntrye et al. . The difference is quite large, but the gating kinetics parameters for these two models also differ. In our model we take mS cm−2. Indeed, the model developed by McIntyre et al.  was made for the rat fibre and the model considered by Hennings et al.  was made for the human fibre. This may be the root of the difference in parameter values. We take more or less intermediate values for these parameters since the aim is not to build a human-specific or a rat-specific model. Internodes are free from N+-channels in our model (table 2).
However, there is also evidence that a small number of K+ channels are present in the mammalian node with kinetic properties, different from those of delayed rectifier type (fast) Kv channels . These ‘slow’ K channels are blocked by tetraethylammonium (TEA), but not by 4-aminopyridine. Mammal nerves actually contain more types of potassium channels [30,31], but as the properties of different channel subtypes are not known in detail, it is a common practice to roughly fit all the channels into the two major subtypes (slow and fast) .
Therefore, we only consider the slow potassium current, in the nodal currents (with mS cm−2) while internodal regions feature both the slow potassium current (IKs) and fast potassium current , but no Na+-currents. The kinetics of the fast K+-current is adopted from Hennings et al.  without modification.
The nodal slow potassium current is assumed to reduce DAPs , to take part in the hyperpolarizing afterpotential  and to play an important role in accommodation, i.e. the damped rhythmic response to long depolarizing stimuli.
To briefly recapitulate, slow K+ currents are confined to nodes of Ranvier whereas fast K+-channels (Kv) are clustered in juxtaparanodal regions. Fast Kv channels are present in other internodal regions, but their density is much lower outside the paranode. Since paranodes are short in comparison to the main internode, K+-conductance in the latter is still significant despite the lower channel density. It seems safe to assume that approximately 60 per cent of all internodal K+-conductance is subserved by juxtaparanodal Kv channels and the remaining 40 per cent by Kv channels in the rest of the internode. This is taken into account in our model in the following way: for 60 per cent of the fast K+-current, the Nernst potential for K+ is calculated from the paranodal [K+]e, and for the remaining 40 per cent from the internodal [K+]e (figure 1). The average fast K+ conductance over the entire internode including the paranodes is mS cm−2.
3.2.3. Channel kinetics details
Kinetics of all the gating variables has the classical form: 3.4
Voltage dependence for each rate constant ξ in all αx, βx except for βh is parametrized as follows: 3.5while for βh it is 3.6
A, B and C parameters used in the model are given in the table 3.
Parameters of αs and βs for the slow potassium current activation, variable s differ significantly from one modelling study to another (compare e.g. [23,24,27]), though most authors refer to the work of Schwarz et al.  for the actual experimental data. The form of dependence used in our model approximately follows that of Schwarz et al. , but is slightly sharper to come closer to that of McIntyre et al. .
3.3. Na+,K+-pump in axonal and glial membranes
Following Truskey et al. , the sodium flux through the pump is modelled as two coupled Michaelis-Menten kinetics of order 3 for Na+and 2 for K+ 3.7 3.8 and 3.9Here Kmx is the Michaelis constant for the ion x, and parameters Lx describe the competition with the other ion for the binding site. The potassium flux is taken to be JK =− 2/3JNa.
As discussed by Ransom et al. , axonal and glial pumps differ in their affinities to K+ and Na+. The axonal pump is mainly regulated by the Na+ concentration in the axoplasm, [Na+]i , whereas the glial pump is regulated by the K+ concentration in the intercellular space, [K+]e. This difference is taken into account by using different KmNa and KmK values for glial and axonal pumps (table 4). We assume maximal specific currents to be the same for glial and axonal membranes, (JNamax = 0.7 × 10−7 mmole cm−2s). Parameters LNa = 10−2 and LK = 10−3 are also assumed to be shared by both pumps. The values for JNamax, LNa and LK were chosen by searching for parameter values that give physiologically relevant steady-state values of [Na+]i and [K+]e, current density, produced by the pump and without affecting normal responses to stimulation. JNamax = 0.15 × 10−6 mmole cm−2s for the Na+,K+-pump of kidney epithelial cells , and it seems natural that it should be lower for nerve and glial cells.
3.4. Ion currents in the Schwann cell
Owing to a large number of different types of K+-channels, the Schwann cell membrane is permeable for K+ in a wide range of potentials. This plays a significant role in the exchange of K+ between axon, periaxonal space and myelin . This fact is taken into account in the following way. We add a variable for the internodal glial potential, Vg, with the following dynamics: 3.10where cig = 1 µF cm−2 is the capacitance of the internodal adaxonal (facing the periaxonal space) glial membrane, gKg = 2 mS cm−2 is the overall potassium conductance of the glial membrane, glkg = 0.01 mS cm−2 is the non-specific leakage conductance through glial membrane, Vlkgrest = −85 mV is the leakage reversal potential, and Ipump,g is the current density generated by the glial pump, which was calculated from JNa,g. Values for parameters gKg, glkg and VKrest are based on general considerations and were chosen so as to keep steady [K+]e in the internode at physiologically reasonable values (approx. 1.9 mM). It has to be noted though, that there is no precise experimental measurements for the periaxonal concentration, so this value is also to some extent a guess. In other compartments, the steady potassium were the following: in the node, 3.7 mM and in the paranode, 5.5 mM. Even at resting conditions, [K+]e in the paranode was elevated in comparison to other compartments, which is not surprising given the high density of potassium channels there.
3.5. Passive exchange of ions
Extracellular potassium [K+]e and intracellular sodium [Na+]i are dynamic variables in our model. Respective differential equations for these ions include currents (I = JF) through ion channels and the Na+,K+-pump and passive exchange between the nodal, paranodal and internodal compartments (intra-axonal for Na+ and periaxonal for K+). This diffusive exchange is determined by the geometry of the compartments and the diffusion coefficients for the corresponding ions in aqueous solutions are taken to be DK = 1.96 × 10−5 and DNa = 1.33 × 10−5 cm2 s−1 .
In the node there is also linear exchange [X] = …+ ([X] − [X])/τ of the [K+]e with the external bath medium with [K+]o = 3.2 mM. Diffusion from the perinodal space to the outside medium is limited by the microvilli of the Schwann cells enclosing the nodal gap. This exchange can therefore be rather slow. The model also considers a linear equilibration term for [K+]e in the internodal regions, which sums up exchange through Schmidt–Lanterman incisures, ion fluxes through Na+,K+,2Cl−-co-transport (e.g. ) and others not considered in detail in our model. We choose the characteristic time constants of these exchange processes for both node and internode to be τ = 2 × 103 ms. Similar exchange rates are used for [Na+]i in the internode, but the characteristic time for Na+ is longer, τ = 5 × 103 ms, as there is no diffusion through the incisures, only fluxes through the exchangers.
3.6. [K+]e-dependent changes in internodal leakage resistance Ril
The model described so far is able to account for the accumulation of K+ under the myelin in paranodal and internodal regions as well as in the restricted space around the node. However, conduction of long spike trains can result in structural changes in the nerve, which can be described as functional paranodal demyelination.
In this study, the demyelination is modelled as potassium-dependent decrease in the internodal leakage resistance Ril. This combines in a simple way several different mechanisms, such as separation of myelin lamellae and swelling of the paranodal loops, possible loosening of the paranodal seal and widening of Schmidt–Lanterman incisures, as all of these processes deteriorate myelin resistance.
We introduce a [K+]e-dependent variable w, a ‘widening factor’ with a linear relaxation process (3.12), characterized by the equilibrium value w∞ (3.13). A resting value of internodal leakage resistance is divided by this variable w to get actual Ril. The time constant τRil is 1000 ms as a somewhat arbitrary but presumably reasonable value: 3.11 3.12 and 3.13where CKp is [K+]e in the paranode and is the normalization parameter equal to 1 mM. The curve in equation (3.13) has the basic S-shape. This choice was made based on the following tentative mechanism: main changes in Ril are likely to be linked to the width of paranodal loops and (in the PNS) Schmidt–Lanterman incisures. Normally these structures are in a resting, ‘closed’ state and they can ‘open’ up to a limited extent following the rise in paranodal [K+]e.
The inflection point CK50 determines the myelin sensitivity to [K+]e, several values are tested, CK50 = 6.5, 7.0, 7.5, 8.0 mM. The steepness parameter γ is equal to −25 and describes a rise of w∞ from approximately 1 to 5 over the range of concentrations from 5 to 10 mM. Resting in our model is 50 MΩ, which is the same as in Bostock et al.  (figure 2).
3.7. Simulation details
All stimulation impulses used in the model were 0.1 ms in duration. Threshold current densities were established for each value of CK50. This was 9.0 µA for CK50 = 6.5 mM and 8.2 µA for CK50 ≥ 7.0, which is consistent with values in Bostock et al. .
Differential equations were numerically integrated by explicit Runge–Kutta–Chebyshev method for stiff systems with adaptive step size . Custom software was developed to describe and solve the equations in the model. The software will be available as open source at http://senselab.med.yale.edu/modeldb/ or other public web site.
The common pattern of model behaviour is illustrated in figure 3: prolonged excitation leads to a rise in [K+]e in the periaxonal space, which in turn produces changes in the internodal leakage resistance Ril according to equations (3.11)–(3.13). Decreasing Ril leads to even higher depolarization of the internodal membrane during each spike and henceforth to more K+ accumulation in the periaxonal space, forming a positive feedback loop. After a significant decrease in Ril, the threshold rises above the stimulation current and excitation fails to evoke APs. During the silent period, the membrane potential does not rise to high positive values, the K+-currents remain mostly deactivated and the paranodal [K+]e cleared from the periaxonal space, thus letting the myelin structure relax back to the high resistance stage. This lowers the threshold and allows for the repetition of the cycle (figure 3).
This is seen in the onset of a batch of spikes shown in more detail in figure 3b at an expanded timescale. While paranodal [K+]e remains relatively stable, Ril rises and only when it reaches approximately 36.7 MΩ, the stimulation impulses exceed current threshold value. Resumed spiking leads to a new rise in paranodal [K+]e. With a delay, the Ril begins to fall again. But now because of the high paranodal [K+]e, the fibre is depolarized and excitation breaks at lower Ril values than it has started (approx. 35.7 MΩ), which is illustrated at an expanded timescale in figure 3c.
When Ril is decoupled from the change in [K+]e (setting wmax = 0 in equation (3.13)), the K+ elevation still occurs and reaches a steady level, but no excitation breaks are observed at the interpulse intervals (IPIs) used.
4.1. Changes with stimulation amplitude and frequency
The decrease in Ril resulting from the rise in paranodal [K+]e only leads to a limited increase in the excitation threshold. Therefore, stimuli of different amplitudes are expected to lead to different patterns of activity.
There is also a relation between stimulation frequency and stimulus amplitude. Since most of the outward K+ current occurs during the repolarization phase of APs, it is obvious that the extent to which K+ accumulates depends on the number of APs per unit time, which is defined by the frequency of stimulation. It seems therefore natural to expect the pattern of excitation block to change with IPIs as well as with stimulus amplitude.
This type of relation is illustrated in figure 4. Figure 4a,c demonstrates the case of low-amplitude stimulation (90% of threshold at resting state) while figure 4b,d represents the case of a higher stimulation amplitude (110% of threshold at resting state). Stimulation amplitudes are taken relative to the minimal current amplitude required for a model to generate an AP from the resting (steady) state. With ‘sub-threshold’ 90 per cent stimuli series of APs are generated because of the DAPs that lower the threshold for the second and subsequent stimuli. Therefore, here we speak of ‘sub-threshold’ stimulation in the sense of excitation from the resting state.
In the case of sub-threshold stimulation intermittent excitation blocks occur already at a 15 ms IPI, whereas the delay before block shortens with the decrease in IPIs. Conduction block is permanent at 2 ms IPIs (figure 4a). In the case of stimulation with 110 per cent of the threshold, conduction block develops only at 8 ms IPIs, the block being permanent at 2 ms too.
Figure c–d show 1 s slices from the tracks displayed in figure 4a,b, providing details of the fine structure of the off-duty (silent) regions. The observed patterns of conduction block can be phenomenologically listed as follows: (i) no blocking, (ii) long batches of spikes intervened by relatively short blackouts, (iii) on-duty (spiking) intervals are shortened, (iv) on-duty intervals are split into clusters (e.g. at 8–6 ms IPIs for 90% excitation), (v) only short batches of APs are present, (vi) total block. Note the rather abrupt changes in patterns in a transition from 5 ms to 4 ms IPI for 90 per cent stimulation and from 6 to 5 ms for 110 per cent stimulation. These transitions will be further discussed below.
4.2. Changes with sensitivity to [K+]e
The effect of sensitivity to [K+]e is similar to that of the stimulation amplitude, which is illustrated in figure 5. With the w∞ shifted to the left (CK50 = 6.5 mM) conduction block develops at longer IPIs than in the case of larger CK50 = 0.75 mM. However, clustering within the batches does not occur as it does for sub-threshold stimulation amplitudes (figure 4a,c). Not only does the intermittent block first occur at shorter IPIs for higher CK50 values, but the total conductance block also follows that pattern. As shown in figure 5, total block develops at 3 ms IPI in the case of CK50 = 0.65 mM and does not develop even at 2 ms IPI in the case of CK50 = 0.75 mM.
Relations between the maximal interval at which conduction block occurs, stimulation amplitude (in terms of the resting state threshold) and sensitivity to [K+]e are summarized in figure 6. The maximal interval versus simulation amplitude graphs are shown for different ‘sensitivities’ to [K+]e, i.e. different values of parameter CK50 in equation (3.13). It can be seen that for lower stimulus amplitudes conduction block develops at larger IPIs. This dependence is steeper for lower values of parameter CK50 and gradually flattens for higher CK50 (lower sensitivity to [K+]e). At CK50 > 7.5 mM, the blocking is limited by the sensitivity of Ril to [K+]e because the level of [K+]e could not reach the high critical value for most physiological IPIs. On the contrary, for low CK50 values (high sensitivity), the [K+]e level easily reaches the critical value for a given stimulation amplitude even at moderate IPIs. All the following results are given for CK50 = 7.0 mM as a compromise between the two types of regulation.
4.3. Interspike intervals: how many stimulation pulses are missed?
More information on the patterns of intermittent excitation block can be gained from the distribution of interspike intervals (ISIs). It seemed natural to normalize these ISI to the IPIs. With this normalization, an ISI equal to 1 means that the model fibre responds to each stimulus, an ISI equal to 10 means that one AP is generated for each 10th stimulus. Such distributions of ISIs for IPIs of 2–7 ms and stimulation amplitudes corresponding to 90, 110 and 150 per cent of the threshold are shown in figure 7.
It seems convenient to start from the case of a well over-threshold stimulation (150% threshold, right column in figure 7). At 7 ms IPIs there is no excitation block—all ISIs are equal to 1. At 6 ms ISIs around 20 can be observed, but the column for ISI = 1 is still around 1. This indicates long batches of spikes rarely intervened by relatively long off-duty intervals. Further decrease in IPIs to 5 ms leads to a slight shift of the ISIs to the left, which corresponds to longer off-duty intervals. Also the second group of columns is higher than for 6 ms suggesting shorter on-duty periods. In both cases, the group of columns corresponding to silent periods is relatively narrow. However, at 4 ms the group of columns becomes wider, indicating instability in silent period lengths. At 3 ms, the group of columns apparently do not change, but a new column for ISIs equal to 2 appears, reflecting rhythm transformation within the on-duty periods—the fibre generates a spike only for every other stimulus within a batch of spikes. The presence of the column for ISIs equal to 1 reflects a transient period before the first conduction block or a stable pattern within batches of spikes, we did not consider this period in detail. At the 2 ms IPIs there are columns for ISIs equal to 1, 2 and 3 and the second group of columns shifted to the right. This histogram describes the regime with short batches of spikes (some impulses failing), intervened by long and variable inter-batch intervals.
The overall picture is similar for the dynamics at near-threshold stimulation (110% threshold, central column in figure 7). In this case, widening of the right group of columns in the ISI histograms occurs at 5 ms IPI. At 4 ms there is a switch to shorter IPIs and a column for ISIs equal to 2 appears. This means that at 5 ms the batches of spikes are not decimated, but the silent intervals between them have very variable lengths, whereas at 4 ms the length of inter-batch intervals is less variable and also shorter, the spikes in the on-duty intervals being decimated (i.e. not all impulses lead to spikes). At 2 ms there is again a wide spectrum of different ISIs, but the columns for ISIs equal to 1 and 2 still dominate. However, at this combination of parameters, the conduction block is total and the histogram only reflects a transient period from spiking to a blocked state.
Stimulation with subthreshold pulses (90% threshold, left column in figure 7) apparently just shifts the same picture even more towards larger IPIs. The set of columns in the histograms for this case is rather complex, and it seems fairly difficult to infer more information from them.
Another angle to look at changes in spike patterns is to study projections of phase trajectories. Such projections on the plane CKp, Ril are shown in figure 8. We recall that CKp is the potassium concentration in the extracellular space of the paranodal region and , where w is a ‘myelin structure’ variable that decreases internodal leakage resistance.
Behaviour of the system passes through the following stages: (i) when the IPIs are wide enough for a given stimulation amplitude, a steady value of [K+]e (e.g. at 9 ms interval, 110% threshold stimulation) is reached, (ii) with decrease in IPI a relatively wide limit cycle emerges (e.g. at 7–6 ms and 110% stimulation in figure 8), (iii) with further decrease in IPI, the area taken by the limit cycle is filled with trajectories (5 ms at 110% stimulation) and after that (iv) the area gradually collapses with the further decrease in IPIs. The total block of conduction develops at 2 ms IPI.
Oscillations of [K+]e and w here reflect intermittent periods of excitation block in the node. A large open cycle in these plots indicates long inter-batch intervals at which the [K+]e could decrease significantly, almost down to a resting value. This is illustrated in more detail in figure 9.
Short inter-batch intervals are on the other hand indicated by smaller cycles as the [K+]e does not have enough time to decrease significantly in these periods. Thus, a large filled cycle type trajectory corresponds to a wide distribution of ISIs. Phase trajectories thus provide a more detailed picture of how the changes in behaviour pattern follow the changes in the IPIs: regular oscillations become chaotic and then the enclosing phase space collapses, leaving only small-amplitude oscillations in [K+]e i.e. short ISIs and eventually, a conduction block.
The details of transition from the large cycle (long inter-batch intervals) via the filled cycle (wide group of ISIs) to the small cycle (short inter-batch intervals) are also examined with the wavelet analysis of the system dynamics. Time evolution of [K+]e was extracted from the model output and resampled to a fixed time step with linear interpolation and then analysed with continuous wavelet transform (CWT). Details of the CWT can be found elsewhere . (In essence, CWT is akin to the windowed Fourier transform and provides with time-resolved spectral information on the analysed signal.)
Spectrograms for the [K+]e realizations in the case of 110 per cent of threshold stimulus amplitude and 6, 5 and 4 ms IPIs are shown in figure 10. It can be seen that for 6 ms there is a stable low-frequency rhythm and at 4 ms there are relatively stable high-frequency rhythms with a slowly changing frequency. Interestingly, at 5 ms the dynamics of [K+]e hops irregularly from a higher frequency rhythm to a lower frequency rhythm.
We introduced a model of extracellular potassium accumulation in the periaxonal space of the myelinated nerve fibre during prolonged repetitive stimulation. The accumulation of potassium was linked to the myelin resistance, which was variable in the model. Thus, excitability of the fibre was affected by potassium not only via the depolarization block, but also through increased electrical leakage through the myelin. The model has an electrical scheme for a space-lumped fibre, which is superimposed on several concentrational compartments.
Nerve excitability is a function of myelin resistance. Therefore, making it dependent on sub-myelin potassium concentration has lead to new patterns of nerve behaviour. At high frequencies of stimulation (greater than 150 Hz for near-threshold stimulation and model parameters considered ‘normal’), the model fibre was periodically or completely blocked, depending on the frequency and amplitude of stimulation. This was not observed if the myelin resistance was uncoupled from the potassium concentration and remained constant.
It is known that axons can transform the frequency of stimulation pulses if it is high. Our model suggests the following mechanism for that: a train of APs leads to potassium build-up under the myelin owing to a limited activation of the paranodal potassium channels. This leads to myelin restructuring and a decrease in its resistance. At first, this forms a positive feedback on potassium accumulation, because the axon in the paranodal region becomes more depolarized during the course of APs. However, myelin resistance gradually becomes small enough for the excitation threshold to exceed stimulation amplitude. During this period excitation pulses only cause low-amplitude depolarizations instead of full APs and potassium is cleared by the axon and the Schwann cell via pumps, K+-channels of the Schwann cell and other membrane transporters. With decrease in extracellular potassium, the myelin structure slowly restores its original electrical properties and the cycle reiterates.
In our model, at relatively low-stimulation frequencies, the intervals of spiking responses were long and those of non-spiking responses were short. With increase in frequency, the former became shorter and the latter wider. At certain frequency, depending on stimulation amplitude, the spiking and non-spiking intervals became chaotic in their duration and with further increase in stimulation frequency the fibre switched to another regime, with short spiking and short non-spiking intervals up to 2 : 1 frequency transformation (one spike per two pulses). Even further increase in stimulation frequency blocked excitation completely. Though the fact of excitation block at high-frequency stimulation is well known, the frequency-dependent patterns of such block predicted by our model have not been experimentally reported.
Similar patterns of intermittent failure have been demonstrated by Krnjevic & Miledi [39,40] in the rat. Oscillations of extracellular potassium during intermittent block can be seen in the study of Smith , figure 4. This study was performed on an un-myelinated nerve, but the conduction failure must be related to a similar K+-mediated mechanism, since it takes place where the sheath is thinner and the periaxonal space is reduced .
Intermittent failure patterns similar to the described in our paper are more investigated in C-type unmyelinated fibres, for example, see [43,44]. It must be noted that in the case of C-type fibres, the failure appears at rather low frequencies (4 Hz), whereas in our model of the myelinated fibre the failure appears at higher frequencies (greater than 100 Hz). The reason why such intermittent patterns have not been studied for myelinated fibres is not clear, but it is likely to be observable in them as well.
Stimulation-dependent conduction block is important in several aspects. One of them is high-frequency stimulation treatment for Parkinson's disease and dystonia in deep brain stimulation [45,46]. Block of axonal conduction during high-frequency stimulation was shown in vivo experimentally by Jensen & Durand , though the mechanisms of such blocks remain unknown. This problem has been addressed in a modelling study of Bellinger et al. . In that model, the block was achieved through potassium-dependent depolarization block. The authors modelled potassium accumulation in the sub-myelin space. However, real potassium uptake from the sub-myelin space can be more efficient than that of the model, e.g. because the pumps of the Schwann cell membrane facing the periaxonal space were disregarded in the model. High potassium concentration (more than 10 mM) was crucial for the block to occur. Our model complements that study with an additional proposed mechanism for the functional block which develops at lower potassium concentrations (6–7 mM).
Conduction block may not take place at physiologically normal conditions without artificial stimulation. However, it may occur in the case of abnormal neuronal activity associated with high-frequency bursts of APs synchronized in a nerve trunk or ischaemic conditions. Such cases facilitate extracellular potassium accumulation and this can lead to excitation block at lower frequencies of spike trains and thus in some cases have a protective function. This kind of block may also be important in the pathological cases of chronic loosening of myelin layers, which would make the myelin sheath more structurally responsive and unstable.
The model discussed in this study has a number of limitations. First, myelin capacitance was a parameter in our model and did not change with changes in myelin structure, though this can be the case in the real fibre. However, from the equation (3.2) this is expected to have an impact very similar to that of just decreased myelin resistance. Second, the model is electrically space-lumped; a spatially extended model with the same channel kinetics and geometry parameters has also been tried with the similar behaviour patterns, but the space-lumped model was chosen for the sake of computational efficiency. Third, the real fibre includes more types of potassium channels (e.g. Ca2+-activated and inward rectifier channels) and transporters (e.g. Na+,K+,2C−-co-transport) that were not included in our model. Though these additional currents are indeed important for fibre function, it seemed safe to suppose that an explicit inclusion of them in the model would not have affected the observed results dramatically and the main mechanism of the excitation block would have remained the same.
We acknowledge support from the European Union through the Network of Excellence BioSim, contract no. LSH3-CT-2004-005137
One contribution to a Theme Issue ‘Advancing systems medicine and therapecutics through biosimulation’.
- Received August 27, 2010.
- Accepted November 5, 2010.
- This Journal is © 2010 The Royal Society