The defect effects on the signal transport of an excitable soft cable

How a local perturbation affects a propagating wave traveling in a homogeneous medium is a general physics question widely investigated in condensed materials. Intuitively, one might expect that a perturbation would suppress the transport ability of the medium if it is quasi one dimensional. This is generically true as defects and impurities influence numerous non-excitable systems such as carbon nanotubes, nanowires and DNA double helixes. However, if the system is excitable, such as a neuron, a defect may generate a highly non-trivial dynamical behavior. In this paper, using the Hodgkin–Huxley model, we explored this diversity generated by locally non-uniform ion channel densities caused by toxins, diseases, environmental disorders or artificial manipulations. These channel density defects could induce several exotic behaviors, in contrast with the normal destructive role of defects in solid-state physics. They may behave as an electric signal generator exhibiting spontaneous or stimulated emissions, as well as trap, reflect, rectify, delay or extinguish propagating signals or be switched to different functions by a signal. Nonlinear analysis and phase diagrams were used to quantify this dynamical complexity. The results may contribute to research on signal manipulation in biotechnology, neuronal diseases and damages, channel distribution-related cell functions and defect dynamics in general excitable mathematical models.


Introduction
Ion channels are crucial for the propagation of electric signals on excitable cells [1]. In the landmark work of Hodgkin and Huxley [2], these channels were postulated as 'molecules with a charge or dipole moment' to account for the voltage-dependent ion permeability of an excitable membrane. Although the types of ion channels and their mechanisms were unknown at the time of that work, these authors had speculated that the density of ion channels is relatively low. After the existence of the ion channel was directly confirmed in the 1970s [3], an increasing number of natural non-uniform channel distributions were uncovered, most of which are believed to play a role in specific neural functions. For instance, aggregations of Na and K channels at the axon initial segment generate action potentials (APs) [4,5] integrate slow subthreshold signals and control neuronal activity [6,7]. The gradient distributions of the same channels in dendrites control back propagation, regulate dendritic excitability and modify synaptic plasticity [8][9][10][11][12]. Even small clusters containing several tens of Ca and K channels are thought to optimize the response to small stimuli during cell signaling [13][14][15][16][17]. Undoubtedly, intrinsic channel density heterogeneity is common in neurons and, in addition to channel type, essential for controlling neural behavior [18].
These heterogeneous channel distributions are created through a delicate channel targeting machinery that includes channel trafficking, retention and endocytosis pathways. Although the number of molecules involved in this machinery is still unclear, several channels are known to anchor to the cytoskeleton through protein scaffolds, specific channel motifs, auxiliary subunits, as well as ankyrin G, βIV-spectrin and filamin A [18][19][20][21][22]. Thus, toxins, mutations, cytokines, abnormal phosphorylation, actin filament depolymerization, damage or diseases that disrupt this targeting machinery or immobilization proteins can readily perturb ion channel densities [10,[23][24][25][26][27][28]. Some of these perturbing factors can be implemented artificially to obtain the desired channel distributions. For instance, latrunculin B and cytochalasin B have been used to depolymerize local actin filaments in order to remove different types of channels simultaneously [28]. Toxins such as tetrodotoxin and tetraethylammonium can be applied through micropipettes to cripple specific types of channels locally on a neuron [10]. Under these perturbed conditions, the distribution of the remaining functioning channels is referred to as an extrinsic channel distribution. If a homogeneous channel density is perturbed locally, its extrinsic distribution will contain a local channel density defect.
To gain insight into the influence of this defect on the electric properties of neurons, one may consider applying the techniques of phase space analysis in nonlinear science. This approach has been instrumental in extracting the key behavior of several neural models from the vector field structure of their reduced dynamics. A typical example is the FitzHugh-Nagumo reduction that combines the fast-and slow-varying variables of the Hodgkin-Huxley (HH) model to reduce the number of variables from 4 to 2 [29,30]. However, the simplicity of this analysis largely rests on the requisition that the reduced system is low dimensional. The inhomogeneous channel distributions discussed below belong to an infinite-dimensional problem, which does not conform to that requisition, and a dimension reduction there usually leads to a completely irrelevant dynamics. On the other hand, one might seek clues from a variety of reported studies that aimed at finding appropriate conductances for the neuronal systems of one's interest. Unfortunately, these studies mainly pertain to generic neurons of moderate channel density variations and are incapable of predicting the influence of the following abrupt local channel density perturbations. Thus, although the complexity of defect dynamics is qualitatively to be expected, it lacks a quantitative analysis and remains elusive. To our knowledge, even basic questions such as those on the critical defect size that can lead to qualitatively different neural dynamics and the number of different dynamics that can emerge under a perturbation have not been answered. These questions and their related neural properties are the central issues this work is intended to analyze and classify.

The shapes of channel density defects
The defect effect caused by the aforementioned extrinsic channel non-uniformity was studied using the HH model. The original constant Na and K conductances therein were multiplied by a weight w(x) to characterize channel non-uniformity along the neuronal axis x (see the appendix). If w(x) was uniform except around x 0 , the channel density would have a defect at x 0 . For instance, the Gaussian weight w(x) = 1 − γ exp[−(2x − 2x 0 ) 2 /σ 2 ] and the square weight w(x) = 1 − γ [1 + (2x − 2x 0 ) n /σ n ] −1 described a density defect at x 0 with the defect width σ and the defect ratio (or relative depth) γ , where 0 γ 1 (figure 1). Both w(x) were analytical functions lying between 0 and 1. While the former had smooth boundaries, the latter could have arbitrarily sharp boundaries for large n. For γ → 0, both weights approached w(x) ≡ 1 and the system returned to the original defect-free HH equations. The neural dynamics caused by these defects was calculated by NEURON, 4 with the parameter values of squid in [2], omitting the leakage current owing to its negligible effect.

Individual channel density defects
Na and K channels underlie the depolarization and repolarization of the neuronal membrane, respectively. For uniform Na and K channel distributions, the equilibrium membrane voltage V m remains at the position-independent resting potential. If Na channels are slightly blocked around x 0 , the resting V m will be hyperpolarized and form a dip at that site. If K channels are inactivated locally instead, the resting V m will be depolarized and form a hump there. With the increased width or depth of the defect, the dip or hump of V m will grow. This deformed V m will remain stable during the Na channel density defect growth, but may abruptly become unstable and begin to oscillate when the K channel density defect reaches a critical size. Beyond that size, the K channel density defect will behave like a pulse generator, eliciting APs regularly (figures 1 and 2(a)). Below but close to that size, some deformed resting potentials can be stimulated by an oncoming AP to fire regularly ( figure 2(b)).
These properties indicate three kinds of phase space structures for the K defect-containing HH equations: (I) a stable fixed point (resting potential); (II) a stable limit cycle (spontaneous emission), i.e. V m (x, t) = V m (x, t + T ) for all x and t with a constant T ; and (III) a stable fixed point plus a stable limit cycle (stimulated emission). For a Gaussian K channel density defect characterized by (γ , σ ), these three phases are divided by the two red solid curves in figure 3(a). These curves reveal the critical defect sizes beyond which one can see qualitatively different neural dynamics. According to the numerical results, the transitions from (I) to (II), (III) to (II) and (I) to (III) look like the supercritical, subcritical Hopf bifurcations and the saddle node bifurcation of limit cycle, respectively. Along the σ -axis at γ = 0.97, the bifurcation through these three phases seems to undergo a Mexican hat scenario illustrated in figure 3(b), whose inset shows the attracting (red) and repelling (green) solutions from the side view. For an Na channel density defect, the long-term behavior is much simpler and only exhibits phase (I). However, such defects may extinguish an oncoming AP, which never occurs on a K channel density defect. If the Na channel density defect has a Gaussian shape, the (γ , σ ) values to the right of the blue solid curve in figure 3(a) indicate the extinction phase of this defect, and to the left of that curve, the transit phase.
The minimum γ to find phases (III) and (II) of a Gaussian K channel density defect is approximately 0.94 and 0.72 and that to find the extinction phase of a Gaussian Na channel  density defect is approximately 0.89. That is, at least 70% of ion channels must be damaged at the defects to generate qualitatively different neural behaviors. These high ratios reflect the fact that ion channels are more abundantly expressed than required for normal neural functions [31]. The pattern of the phase boundaries for Gaussian defects (red and blue solid (a) Phase diagram for individual Na and K channel density defects (blue and red, respectively) of Gaussian (solid) and square (dash) weights. The red curves divide three kinds of long-term neural behaviors on a K channel density defect: (I) a stable fixed point (resting potential); (II) a stable limit cycle (spontaneous emission); and (III) a stable fixed point plus a stable limit cycle (stimulated emission). The blue curves divide two kinds of short-term neural behaviors on an Na channel density defect: the oncoming AP can pass through or be extinguished by the defect. (b) For a Gaussian K channel density defect, the system bifurcates through three phases along the σ -axis at γ = 0.97. The unstable and stable solutions are plotted in red and green, respectively, both on the Mexican hat in the main plot and in its side view in the inset. curves in figure 3(a)) is rather typical of other defects. Even the square defect with extremely sharp boundaries (n = 50) has a similar structure, as shown by the dashed curves in figure 3(a).
Since the square defect has larger areas in phases (II) and (III), it is more excitable than the Gaussian defect. For general defect shapes, there might exist several equally appropriate definitions of defect width. Although the phase diagram will vary with the width definition, its main structure is essentially the same as that in figure 3(a). This insensitivity to defect shape highlights the influence of the diffusion term in the HH equation (see (A.1) in the appendix). This crucial term blurs the defect identity and leads to the 'universal' feature in the phase diagram.

Blended channel density defects
If both Na and K channels are deficient around the same site, perhaps owing to local cytoskeletal disruption, a blended defect is formed. This defect may become a stimulated pulse generator, which will generically elicit APs on the same side of the stimulation (figure 2(c)). However, some defects can elicit APs on the opposite side (figure 2(d)) or with different frequencies on different sides (figure 2(e)). Under incommensurate frequencies, the latter dynamics will become a quasi-periodic motion distinct from the three phases in figure 3(a). Some exceptional defects cease firing after a particular number of pulses and behave like a bursting generator (figure 2(f)). Moreover, some defects will trap an approaching AP and freeze the depolarization at the defect site ( figure 2(g)). This dynamical system has at least two fixed-point attractors, corresponding to two different resting potentials. Some defects resemble an elastic scatter (right in figure 2(h)), where the time lag during the collision is reminiscent of a soft wall collision. Some defects resemble a signal rectifier, allowing a leftward-moving AP to pass through, but blocking a rightward-moving AP (left in figure 2(h)). Some defects resemble a switch, with a low V m hump state (ON) that can be turned to the high V m hump state (OFF) by a signal. The ON state allows an AP to pass, whereas the OFF state blocks all APs (left in figure 2(i)). That is, such defects allow at most only one signal to pass through.
In these patterns, one has to distinguish (A) the pulses induced by the defects from (B) the pulses generated externally to stimulate and probe the defects. In figures 1 and 2(a), the pulses belong to (A) and any initial neural state will evolve to such a kind of cyclic motion, irrespective of whether an external pulse is launched to stimulate the defect or not. In figures 2(b)-(f), the pulses before and after around 10 ms belong to (B) and (A), respectively. If the defects in figures 2(b)-(f) are absent, an external pulse is not able to trigger any ceaseless firing and any initial neural state will converge to the resting state. If they are present, an initial state will converge either to the resting state or cyclic pulse emission, depending on the basins where the initial state is located. The initial states with incoming pulses (those before about 10 ms) in figures 2(b)-(f) are examples that evolve to cyclic pulse emissions. In figures 2(g)-(i), the defects are not pulse generators. After the external pulses of type (B) pass through those defects, the systems converge to different resting states. Note that an initial state here can have two different interpretations. Firstly, it can be a state around the defect area and the externally generated pulse is not a part of the state. Secondly, it can be a state of the whole neuron and different external pulses belong to different inertial states. The first interpretation has a vivid physical picture, while the second one is mathematically more rigorous.
The contour patterns demonstrated in figure 2 reveal several basic types of defect dynamics generated by combined Na and K channel density heterogeneities. Some of these patterns are easy to understand. For example, our intuition might tell us that the one-sided pacemaker in figure 2(c) is created by a two-sided pulse generator with a pulse suppressor on its right-hand side to extinguish rightward pulses. Indeed, an Na defect is closely adjacent to a K defect, as shown at the bottom of figure 2(c). However, a slight modification on each of these defects could largely change the pattern. It is generally hard to predict beforehand which pattern will come out of a given defect combination. Some of the basic defect dynamics in figure 2 are reminiscent of the fundamental neurocomputational properties generated, for instance, by Izhikevic's simple model [32]. Examples include the typical bursting and pacemaker-like dynamics, which can be found in both the simple model and our systems with channel density defects. Nevertheless, while the former is controlled by the input current, the latter is tuned by the defect shape. Note that the neuron lengths considered in figure 2 are sufficiently long, so that most patterns in the panels are robust against small variations of the parameters γ Na , σ Na , γ K , σ K and d. Exceptions are the defect in (f) and the reflectors in (h) and (i), whose parameters still slightly depend on the selected neuron lengths.
The blended defects in figures 2(c)-(i) can be characterized by the widths σ Na and σ K , the ratios γ Na and γ K of the Na-and K-deficient regions and the position d of the K depletion center relative to the Na depletion center. For γ Na = γ K and σ Na = σ K , this five-parameter space is projected onto a three-dimensional diagram in figure 4. The red curves on the (γ , σ ) plane correspond to the Hopf bifurcations at different d's. The (γ , σ ) values to the left of each red curve have a resting potential solution, whereas those to the right indicate the spontaneous emission of APs. The blue curves denote the transit-extinction phase boundary in analogy to that in figure 3(a). The (γ , σ ) values to the left of each blue (green) curve allow a leftward-(rightward-)moving signal to pass, whereas those to the right block transmission. The existence of the non-zero narrow area enclosed by the blue and green curves at the same d, say d = 5 × 10 4 µm, allows us to find the rare rectifier defect (left in figure 2(h)).
Apparently, the spontaneous pulse generator regime denoted by phase (II) in figure 3(a) shrinks to a smaller area in figure 4 and the stimulated pulse generator regime represented by phase (III) in figure 3(a) even disappears in figure 4. This excitability suppression is due to the involvement of the Na channel deficiency in the blended defect. The suppression extent is significant for Gaussian weights with long tails (red solid curves in figure 4), but weak for square weights with sharp boundaries (black dashed curves in figure 4). Lifting the blue and red curves to the corresponding d's results in the formation of an outer cyan extinction funnel and an inner orange Hopf bifurcation funnel plotted in the three-dimensional space (figure 4).
While the former converges to the rightmost blue curve on the (γ , σ ) plane at d → 0, the latter disappears in the same limit. The disappearance at d ≈ 0 is because the Na-and K-deficient regions almost overlap each other and the depolarization effect of the K deficiency is nearly completely compensated by the hyperpolarization effect of the Na deficiency. In this limit, the blended defect cannot become a pulse generator and the system looks like it is defect-free. But it differs from the defect-free case in that this defect can extinguish an approaching signal. The size of the extinction funnel indicates the suppression ability of the blended defect or, more explicitly, how likely a defect chosen from figure 4 can extinguish an oncoming signal. Moreover, the size of the Hopf funnel represents the defect excitability or how likely it is that a chosen blended defect will be a pulse generator.
If the restriction σ Na = σ K is lifted, the shapes of the two funnels in figure 4 will vary with σ Na and σ K independently. For σ K < σ Na , the blended defect becomes less excitable and the whole Hopf funnel under a given σ Na will shrink to zero at σ K → 0. This limiting case corresponds to a pure Na channel density defect, which will never generate any pulses. In contrast, for σ Na < σ K , the defect becomes more excitable and the Hopf funnel under a given σ K will deform to a d-independent perpendicular funnel at σ Na → 0. In that case, only a pure K channel density defect is present. The behavior caused by this defect is obviously independent of its distance d to any virtual Na density defect of zero width. The projected curve of that d-independent funnel on the (γ K , σ K ) space is exactly the red solid boundary of phase (II) in figure 3(a). Interestingly, the stimulated emission denoted by phase (III) in figure 3(a) does not exist in the phase diagram of blended defect in figure 4, but appears and grows up gradually in the limit σ Na → 0 when the suppression effect of the Na density defect diminishes. Finally, the funnel deformation trends under different defect relative depths, γ Na = γ K , are similar to those under σ Na = σ K .

Asymmetric channel density defects
The asymmetry of a defect can affect other neural transport properties such as the signal propagation speed. This effect can be readily detected by measuring the speed of an AP running on a repeated asymmetric defect like a ratchet, as the triangular example of period L and ratio γ in the inset of figure 5. The main plot of figure 5 shows that this speed will rise slightly above the standard AP speed on a defect-free neuron (green line) for the K ratchet and decline considerably below that for the Na ratchet. For L larger than the typical AP width, around 2 × 10 4 µm, of the squid giant axon, the AP stays mostly on a rising or falling ramp during its running time. When L exceeds the values marked by the crosses, the AP will be extinguished because the channel densities there are too dilute to support signal propagation. The crosses appear only for those channel density defects with a large γ Na close to 1. It reflects again the fact of abundantly expressed ion channels, as the large γ for observing phases (II) and (III) in figure 3(a) indicates. A clear trend illustrated in figure 5 is that a leftward-moving signal (blue) always propagates faster than its rightward-moving (red) counter signal on an Na ratchet. The speed difference can reach 1 m s −1 for γ Na = 1. It is comparable with the conduction velocity difference 12-18 m s −1 under a temperature variation of 6.3-18.5 • C [33]. For infinitely dense ratchet zigzags at L → 0, the ramp period shrinks to zero. The traveling AP will behave like a shape-invariant soliton running on a defect-free neuron. Nevertheless, its speed will differ from the defect-free speed (green line), since the total channel amount at L → 0 is less than that of the defect-free case if γ > 0. How the speed approaches its limiting value at L = 0 cannot In the inset, a rightward-moving signal (red) and a leftward-moving signal (blue) travel on a ratchet defect of ratio γ and period L. These traveling signals are externally generated and not induced by the ratchet defects. The main plot depicts how the rightward-moving (red) and leftward-moving (blue) speeds vary with γ and L. Therein, γ Na and γ K marked on the curves denote the defect ratios of the corresponding Na and K ratchets, respectively, and the green straight line represents the defect-free speed.
be solved numerically, owing to the increasing non-smoothness of the differential equation at L → 0. Consequently, the curve information is missing at small L in figure 5.

Discussion
A neuron is an intriguing biological structure, which is frequently integrated into and/or manipulated using manmade solid devices for applications in bionanotechnology and neural engineering [34]. For instance, nanowire arrays have been put on neurons for non-invasive and highly sensitive detection, stimulation and inhibition of neuronal signal propagation [35]. Carbon nanotubes have been integrated with neurons to promote neuronal electrical activity [36]. Furthermore, reconstitution methods have been developed for inserting voltagegated ion channels into cell-sized giant unilamellar vesicles, allowing for the study of the electric dynamics of different ion channel distributions [37]. Clearly, the effects of inhomogeneous channel density demonstrated above may further enrich the manipulation versatility in these newly emerging fields. Experimentally, it should not be difficult to observe the neural dynamics predicted in this study, since the defects considered here are about a centimeter wide and located on a squid giant neuron, which is typically several centimeters long. A defect of this scale can be precisely positioned and controlled, for instance, by micropipettes filled with solutions of channel blockers or filament depolymerization agents [10]. The injected solutions will diffuse and most likely generate a density defect close to a Gaussian form as depicted in figures 1 and 2. Such an experiment should also be feasible for other thinner neurons, as long as the channel density perturbation can be controlled locally.
In addition to applications in engineering, the present results might provide information on some neuronal diseases. It is well known that improper channel localization can cause communication problems in neuronal networks [18], and channel dysfunction can lead to diseases in many tissues [38]. However, most reports on such diseases are limited to macroscopic descriptions and seldom discuss their origin at the subcellular level. In fact, the channels in some diseases could be non-uniformly crippled, blocked or transported during the channel trafficking process. Thus, some abnormal behaviors might be a combined effect of the several basic dynamical patterns shown in figure 2. For instance, recent studies have demonstrated that a reduced expression of A-type K channels in primary sensory neurons led to mechanical hypersensitivity [39] and that genetic elimination of the K channel Kv4.2 in the dorsal horn neurons enhanced sensitivity to tactile and thermal stimuli [40]. Unless the channels on neurons are uniformly eliminated, this reduction might cause hypersensitivity phenomena similar to the irregular firing occurring after the K channel density defect demonstrated in figure 1. To compare other patterns in figure 2 with possible channel density-related neural diseases, subcellular measurements of the distributions of channels, or their anchoring proteins [41], are crucial. Upon comparison, we still need to note that the patterns in figure 2 are for unmyelinated squid giant axons. If a channel density disease occurs on a myelinated axon, its dynamics will be more complex, since the system contains more degrees of freedom. Nevertheless, some complex dynamics might share common properties with and could be more easily understood through the simpler patterns in figure 2.
Another intriguing comparison is between the current findings and those in other excitable systems. It is widely known that heterogeneity is a big issue of concern in excitable systems such as cardiac tissues. Therein, heterogeneous cellular types in a multicellular system can cause active activities. For instance, modeling studies have shown that passive defects, such as non-excitable fibroblasts, coupled to excitable ventricular myocytes can cause spontaneous oscillations [42][43][44], which has been demonstrated in experiments as ectopic excitations [46,52]. Normal ventricular myocytes coupled to ischemic myocytes cause spontaneous oscillations [47,48] due to elevated resting potential of the ischemic cells. In the present study, reducing the local density of K channel in the HH model also elevates the local resting potential. Thus, the pacemaker patterns of neurons in figures 1 and 2(a) can be regarded as a subcellular version of the oscillatory behaviors of heterogeneous tissues. However, to our knowledge, other patterns in figure 2 have not been reported in tissues. It might be ascribed to two differences between the current neuron system and usual tissue models. Firstly, whereas the voltage dynamics in figure 2 evolves in a continuous space along the axon, that in usual tissue models lives in discrete networks. Secondly, usual tissue models do not contain cell types playing the role of an excitability suppressor such as the Na channel defect in neurons. From this aspect, the dynamics of excitable neurons seem to be more versatile than those of excitable tissues.
an inhomogeneous position deviating from a uniform phase. This slight deviation usually has a dramatic impact on the transport properties and wave propagation of a system, especially a quasi-one-dimensional system. Intuitively, a defect would be expected to suppress the mass and energy flow and to monotonically reduce the conductance or information transfer through a system; this is indeed true for normal solid materials, which are usually non-excitable [49][50][51][52][53]. However, how a defect will change the behavior of an excitable one-dimensional cable, such as a neuron, is highly non-trivial and hard to predict. To investigate this, we focused on a family of basic local ion channel density perturbations in this study. A systematic survey of the defect parameter spaces revealed several fundamental short-and long-term defect dynamics in the HH model and yielded the critical defect size for qualitatively different neural dynamics. Our results may (i) serve as a simple example illustrating the typical defect dynamics in general excitable cables, (ii) give hints on how to utilize defects for manipulating neuronal signals for biotechnological applications, (iii) quantitatively elucidate disease-related irregular firings at the subcellular level and (iv) enrich our understanding of the challenging issue of channel density heterogeneity-induced neural behaviors [18]. We expect that the present study will stimulate further theoretical and experimental works in these areas.
Na ion channel, as well as the gate of the K ion channel, change with V m , respectively. They follow the first-order kinetic dX dt = α X (1 − X ) − β X X (X = m, n, h) (A. 5) with the empirical V m -dependent rate constants