The threshold firing frequency of a neuron is a characterizing feature of its dynamical behaviour, in turn determining its role in the oscillatory activity of the brain. Two main types of dynamics have been identified in brain neurons. Type 1 dynamics (regular spiking) shows a continuous relationship between frequency and stimulation current (f-Istim) and, thus, an arbitrarily low frequency at threshold current; Type 2 (fast spiking) shows a discontinuous f-Istim relationship and a minimum threshold frequency. In a previous study of a hippocampal neuron model, we demonstrated that its dynamics could be of both Type 1 and Type 2, depending on ion channel density. In the present study we analyse the effect of varying channel density on threshold firing frequency on two well-studied axon membranes, namely the frog myelinated axon and the squid giant axon. Moreover, we analyse the hippocampal neuron model in more detail. The models are all based on voltage-clamp studies, thus comprising experimentally measurable parameters. The choice of analysing effects of channel density modifications is due to their physiological and pharmacological relevance. We show, using bifurcation analysis, that both axon models display exclusively Type 2 dynamics, independently of ion channel density. Nevertheless, both models have a region in the channel-density plane characterized by an N-shaped steady-state current-voltage relationship (a prerequisite for Type 1 dynamics and associated with this type of dynamics in the hippocampal model). In summary, our results suggest that the hippocampal soma and the two axon membranes represent two distinct kinds of membranes; membranes with a channel-density dependent switching between Type 1 and 2 dynamics, and membranes with a channel-density independent dynamics. The difference between the two membrane types suggests functional differences, compatible with a more flexible role of the soma membrane than that of the axon membrane.
It is now more than 60 years since Alan Hodgkin categorized the firing behaviour in his classical study of isolated axons from the crab Carcinus maenas [1]. In many respects his experiments still form the basis for analysis of firing patterns in nervous systems. Using threshold dynamics and maximum frequency as parameters, he identified two major classes of repetitively firing axons (he also defined a class of axons which only fired with difficulty, Class 3): Class 1 axons start firing with very low frequency at threshold stimulation, yielding a continuous f-Istim relationship, whereas Class 2 axons start firing abruptly with a relatively high frequency (typically 75 Hz) at threshold, yielding a discontinuous f-Istim relationship. On the basis of a similar categorization mammalian cortical neurons have also been separated into main classes [2], [3], one exhibiting Class 1 characteristics (regular spiking neurons) and another Class 2 characteristics (fast spiking neurons). The former class consists primarily of pyramidal neurons and the latter primarily of interneurons. This differential classification of excitability has been shown to correlate with a differential bifurcation behaviour of corresponding dynamical models [4]–[6] and successfully been used in analysing the coding properties of neurons [2]–[7]. To avoid confusion, and in accordance with the notation of Tateno and Robinson [7], we in the following use the terms Type 1 and Type 2 dynamics when referring to continuous and discontinuous f-Istim relationships, respectively. This classification takes the threshold dynamics of the regular and fast spiking neurons, and that of the Class 1 and 2 axons, into account, but not all behavioural aspects of these classes [8]. The intricate interactions between the many factors involved in the dynamical regulation of neuronal firing are poorly understood [8]. The dominant idea is that different combinations of ion channel types explain the different patterns [9]. In a previous study we proposed a complementary explanation [10], [11]. We showed that both Type 1 and Type 2 behaviour can be simulated in a dynamical model of a hippocampal neuron [12] by varying the membrane density of voltage-gated Na and K channels (i. e. the number of channels per unit of membrane area, reflected in the Na and K permeabilities when all channels are open; see Figure 1 and Methods). The model used was four-dimensional and based on a detailed experimental voltage-clamp study, thus comprising experimentally estimated parameters. The choice of ion channel densities as bifurcation parameters was due to their physiological and pharmacological relevance. Many drugs act by specifically blocking channels and thereby reducing ion channel density both at a somatic and at an axonal level. Perhaps the most used local anaesthetic drug, lidocaine, acts by blocking sodium channels in axons and sensory nerve endings [13]. An increasing number of studies suggest a role for physiological regulation of channel densities, even at a relatively short time scale [14]–[19]. Each type of dynamics, i. e. , Type 1 and 2, was found to be associated with distinct regions in the channel density plane (−) or with corresponding surface areas of an oscillation volume in the −−Istim space (Figure 2). In regions with high and low values (region C1) the model exhibits Type 1 dynamics, whereas in regions with higher values (regions A2 and B) the model generates Type 2 dynamics. A bifurcation analysis (see Methods) showed that the Type 1 dynamics of the model is due to saddle-node on invariant circle (SNIC) bifurcations [10], [11]. Figure 3A portrays such a bifurcation in a V-Istim plot, calculated for the model using region C1 values. The Type 2 dynamics was found to be due to either local Andronov-Hopf bifurcations and/or global double limit cycle bifurcations [10], [11]. The dynamics of the model associated with region B values is due to double limit cycle and subcritical Andronov-Hopf bifurcations (Figure 3B), while the dynamics associated with region A2 is exclusively due to double limit cycle bifurcations (Figure 3C). The double limit cycle bifurcation implies an unstable limit cycle, which is part of a separating structure (sometimes referred to as a separatrix [9], [20]) which separates trajectories turning to a central stable point and those approaching a stable limit cycle. However, preliminary calculations suggested that the bifurcation structure at the border between regions B (Andronov-Hopf) and C1 (saddle node) is more complex than previously described. When more bifurcation parameters are changed (in our case channel densities and stimulation current) a more intricate loss of stability occurs (e. g. bifurcations with a co-dimension 2) [21]. Thus, to obtain a better understanding of the processes we reanalysed the hippocampal neuron model in more detail. Furthermore, we extended the analysis to two other well-described excitable membranes, i. e. , the myelinated axon of Xenopus laevis [22] and the giant axon of Loligo forbesi [23]. We found that oscillations associated with a subregion of region C1 of the hippocampal model show Type 2 dynamics, and that the oscillations of both axon models exclusively show Type 2 dynamics. We investigated the mathematical background to these findings, using techniques from bifurcation theory. The results suggest that the hippocampal soma and the two studied axon membranes represent two distinct types of membrane with respect to the excitability pattern; membranes with a channel-density dependent switching between Type 1 and 2 dynamics, and membranes with a channel-density independent dynamics. The difference between the two membrane types suggests functional differences, compatible with a more flexible role of the soma membrane than that of the axon membrane.
The three membrane models analysed here are all based on voltage-clamp data, and on the formalism originally developed by Hodgkin and Huxley [23] to describe the dynamics of the squid giant axon. The hippocampal neuron model used is that developed by Johansson and Århem [12] to describe small-sized interneurons in hippocampal slices of the rat Rattus norvegicus. The myelinated axon model used is basically the same as that developed by Frankenhaeuser and Huxley [22] to describe sciatic nerve fibres from the African clawed frog (Xenopus laevis) [24], [25]. The giant axon model used is that of Hodgkin and Huxley [23], describing the dynamics of the giant axon of the squid Loligo forbesi. All the models assume that the membrane current consists of a capacitive current (IC) and a three-component ionic current (Iionic), consisting of a Na current (INa), a delayed rectifier K current (IK), and a leak current (Ileak). It should be noted that in all three models the description of the K currents are based on experimentally measured currents, which cannot be regarded as homogeneous, but are most likely the sum of currents passing through several types of voltage-gated K channels. The descriptions of the currents differ slightly between the three membrane types. For instance, the Na and K currents of the squid axon are described using the conductance concept, while the corresponding currents of the myelinated axon and the hippocampal somatic membrane use the permeability concept, developed by Goldman [26] and Hodgkin and Katz [27] (see Methods). Istim is equal to the sum of the capacitive current (IC), charging the capacitor, and the ionic current Iionic. Thus, (1) where CM is the membrane capacitance. To obtain the time-course of V (t), we solve this differential equation numerically, using the expressions for Iionic presented in the section of Methods. For the analysis of the mathematical nature of the oscillatory activity we determine the stationary potentials (Vs), i. e. the potentials at which the system is in a stationary state and, consequently, the time derivatives of all variables are zero. This was done by solving the following equation for different values of Istim, and (or and values depending on model; see Methods): (2) where Iss is Iionic at steady state. The stability of the system in the neighbourhood of the stationary potentials was examined by a linearization procedure as described in Methods. Graphical solutions to Equation 2 are presented in Figure 4. A more detailed version of Equation 2 is given by Equation 17 in Methods. As shown previously [4], [5], [10], [11], some region in the − (or −) plane must be associated with a non-monotonic Iss-V curve for the model to produce Type 1 dynamics when entering the oscillatory regime. This is due to the nature of a saddle-node bifurcation on an invariant circle (SNIC), requiring an Istim interval at which Equation 2 yields three solutions (i. e. Vs1, Vs2 and Vs3; see Methods). Thus, Type 1 threshold dynamics occurs only when (or) is of a relatively large magnitude, giving the Iss-V curve a non-monotonic, N-like shape (Figure 4). Hence, a switch from Type 1 to Type 2 firing dynamics takes place when (or) is reduced (and or remains intact), corresponding to Na channels being blocked. It should here noted that the existence of three stationary solutions of Equation 2 does not guarantee Type 1 dynamics, as has been pointed out previously [5] and will be seen in the following. In our previous examination of the hippocampal model (see Figure 2), we defined the C area as the region where the model shows a non-monotonic Iss-V relationship (and Equation 2 yields three stationary potentials at some Istim), with area C1 representing the subregion associated with oscillations. As mentioned above, most of this region is associated with Type 1 threshold dynamics. A more detailed analysis reveals, however, that for density values close to the border of the B area the model demonstrates Type 2 dynamics (Figure 5). This can be shown to be due to an Andronov-Hopf bifurcation when the most negative stationary potential (Vs1) becomes unstable as indicated in the bifurcation diagram of Figure 6 (cf. Figure 3A). Figure 7A depicts an oscillation map in the channel density plane, on which minimum frequencies are indicated. As seen on Figure 7B, the line delineating zero frequency deviates from the C1 border at = 14 µm/s and forms a separate, narrow region below this border at higher densities of the Na channel. Below we will denote this subregion C1b and the remaining, larger, region of C1, associated with Type 1 behaviour, C1a. The oscillation map revised accordingly is shown in Figure 7B, where the border between C1a and C1b represents a curve in the ion channel density plane at which Bogdanov-Takens bifurcations occur [28], [29] (see Methods and Table 1). (For the role of Bogdanov-Takens bifurcations in the Hodgkin-Huxley model, integrate-and-fire models and the Morris-Lecar model, see [30]–[33].) It should be noted that in addition to the C1b region, there is another C1 subregion, a narrow strip along the borders to the B, A1 and C2 regions for values below 14 µm/s, associated with a saddle-node bifurcation that causes Type 2 dynamics (non-SNIC), CIc. However, for reasons that will become clear from the analysis of the behaviour of axonal membranes, we will here focus on the dynamics associated with the C1b region. A summary of the regions in the density plane is given in Table 1. The explanation for the deviant (i. e. Type 2) threshold dynamics in the region associated with non-monotonic Iss-V curves (region C1b) becomes evident when the corresponding Istim− diagrams are considered (similar bifurcation diagrams have been successfully used to analyse comparable models [34]). Figure 8 illustrates two such diagrams for = 40 µm/s, one overview and another highlighting the structure at the cusp of the three-root region (which is part of the non-monotonic Iss-V region), with thick continuous lines. The thin continuous line marks points associated with Andronov-Hopf bifurcation dynamics and the hatched line depicts the double-limit cycle bifurcation. The Andronov-Hopf bifurcation line intersects the three-solution region and collides with the saddle-node bifurcation line in a Bodganov-Takens bifurcation point [29]. Hence, at the cusp of the three-root region, the Andronov-Hopf line forms a small subregion, characterised by unstable Vs1 and Vs3. Thus, at these permeabilities and stimulation currents the threshold dynamics is due to subcritical Andronov-Hopf bifurcations and not to saddle-node bifurcations, and the model will show typical Type 2 behaviour with a minimum non-zero threshold frequency. How general is our present description of neuronal models? And how does the density of channels influence the threshold dynamics and firing patterns in other models? To address these issues, two well-described excitable membranes, i. e. , the node in myelinated axons of Xenopus leavis [22] and the giant axon of the squid Loligo forbesi [23] were examined. Both these membranes are similar to the hippocampal membrane with reference to channel composition and kinetics (see Table 2), as can be inferred from the similar mathematical formalism used (see Methods). Nevertheless, the dynamics of both axon membranes were found to show principal differences from that of the hippocampal neuron membrane.
What is the functional reason for the difference in the flexibility of threshold dynamics between the two membrane types defined above (M1/2 and M2, represented by membranes of the hippocampal neuron and the axons)? Type 1 and 2 dynamics per se most likely have different functional roles; Type 1 dynamics being required for low frequency firing and Type 2 being essential for doublet spiking (in hippocampal interneurons, [40]; in dorsal horn neurons, [41]). and for synchronization of firing in coupled neurons (in the synchronization case due to the fact that both a phase advance and a phase delay are possible; the phase response curve being predominantly positive when the oscillations appear via a saddle-node on invariant cycle bifurcation, but both negative and positive in the case of the Andronov-Hopf bifurcation; see [42]). But what about the difference in flexibility between the two membrane types presently discussed? The membranes of the soma and the proximal portion of the axon, which most likely determine the dynamics of the hippocampal neurons analysed here, can be assumed to show a considerable flexibility in their roles as integrative summing points, requiring (transmitter- or trafficking-) induced switchability between Type 1 and Type 2 dynamics (see e. g. [43], [44]). Such flexibility has also been shown in cortical fast spiking (Type 2) interneurons [45]; Type 2 dynamics being changed to Type 1 dynamics when the K channel density (in soma) is reduced in dynamic clamp experiments. Similarly, fast spiking mesenchaplic V neurons have been shown to belong to the M1/2 class [37]. In contrast to the flexible or plastic soma membrane, the axon membranes form passive information transport chains, requiring reliable triggering mechanisms (i. e. high current thresholds leading to rejection of low stimulus noise, and temporal all-or-none responses, meaning that the first spike always occurs early at threshold stimulation) and, therefore, Type 2 dynamics. It should be pointed out, however, that the main value of the axon type dynamics most likely relates to its action when associated with a trigger zone, which likely is the case of the distal process of the dorsal root ganglion. Features associated with Type 2 dynamics, such as subthreshold oscillations and doublet spiking have been postulated to play an important role in pain modulation [46], [47]. Thus, the Type 2 nature of the axon plays a role both in the information propagation and modulation. Clearly, this discussion, based on an analysis of axons from one amphibian (Xenopus laevis) and one cephalopod (Loligo forbesi), cannot be generalized to axons from all animal phyla. As mentioned above, certain axons from the arthropod Carcinus maenas display Type 1 dynamics [1], [23], suggesting that that their cell membranes are of M1/2 type (for a computational analysis of modifying the dynamics of squid axons, see [48]). What about vertebrate axons in general? The phylogenetic bifurcation between the vertebrate and the arthropod lines occurred more than 500 million years ago, allowing a considerable time for specialization of axon membranes. To get information on this issue, we used the data from a voltage-clamp analysis of myelinated rat axons by Brismar [35] to construct and evaluate a dynamical model. The computations suggest that the rat myelinated axon membrane is of M2 type, exclusively displaying Type 2 dynamics (Figure S1 and Text S1). In conclusion, the present analysis shows that axon membranes of two vertebrate and one mollusc species are of M2 type, and axon membranes from one arthropod species are either of M1/2 or of M2 type. More studies are needed to determine whether vertebrate axons mainly are of M2 type or not. It should here be noted that the phylogenetic distance between present day molluscs and arthropods is considerably shorter than that between present day vertebrates and molluscs or arthropods. Mammalian cortical pyramidal cells have been shown to display both Type 1 (regular spiking) and Type 2 dynamics (fast spiking), with Type 1 in majority [7], [44]. Assuming that the trigger zone dynamics is of critical importance for the dynamics of the neuron in toto, the present analysis suggests that the membrane of the trigger zone of the majority of pyramidal cells is of M1/2 type. This also suggests that the assumed trigger zone of pyramidal cells, the initial segment of the axon [49], [50], is not formed by a M2 membrane, contrary to the presently studied axons. A way to experimentally test the hypothesis of a M1/2 membrane as trigger zone in pyramidal cells could be to analyse the results of introducing K channels with the dynamic clamp technique. Such a test is under way. Contrary to the majority of pyramidal cells, mammalian cortical interneurons mainly display Type 2 dynamics (fast spiking). This suggests that the membranes of their trigger zones are either of M1/2 or of M2 type. In the latter case the trigger zone could be assumed to be located in the axon proper; i. e. in the first node of Ranvier or in an initial segment that is more functionally (and structurally?) axon-like than that of the pyramidal cells. A way to experimentally separate between these two hypotheses (whether the trigger zone in interneurons is of M1/2 or M2 type) could be to analyse the dynamics after blocking K channels. Such a test is also under way. As pointed out previously [10], the possibility to modify the threshold dynamics of neurons suggests novel scenarios for the action of channel active drugs such as general anaesthetics; implying mechanisms where selective blocking ion channels in critical neurons induces a switch from one brain state (e. g. associated with consciousness) characterized by certain frequency patterns to another state (e. g. associated with general anaesthesia) characterized by other frequency patterns. Network modelling has shown that such ideas are feasible. Thus, selectively blocking K channels in critical inhibitory neurons (assuming M1/2 membrane trigger zones) in a network of excitatory and inhibitory neurons, distance-dependently connected, can lead to switches from unsynchronised high frequency to synchronised low frequency mean network oscillations [39]. The mechanisms of synchronisation at the network level are still not well understood, but the mechanisms at a cellular level have been extensively studied and a tight connection between the bifurcational structure and the phase-response curve has been established [44], [51]–[53]. Interneurons with Type 2 dynamics have recently been shown to account for the cortical γ-oscillations (20–80 Hz) [54], which are considered to provide a temporal structure for information processing in the brain [55]. Since two of the eigenvalues always are real and negative in the models here discussed, it suggests that the systems essentially are of a two-dimensional character. The decisive variation of the four variables (i. e. V, m, h, n) may then take place on a two-dimensional surface in the four-dimensional variable space. Hence, a model with reduction of variables (as is done in e. g. Fitzhugh-Nagumo and Morris-Lecar models [20], [56]) can give a relatively good description of an excitable membrane. The emergence of a limit cycle following a stability loss can under these circumstances be understood by the Poincaré-Bendixson theorem [57], since the system remains in a finite domain on a curved plane in the phase space. That the system remains in a finite domain is obvious from analysing the variables; the membrane potential V is limited by the reversal potential of Na+ and K+, as well as by the capacitive properties, and the gating parameters are limited by the values 0 and 1. Particularly, the two-dimensional character of the models eliminates more complex types of solutions, such as irregular, “chaotic” solutions or oscillations with two separate frequencies. Nevertheless, local and highly unstable chaos has been reported in a Hodgkin-Huxley system [36], why the models are unlikely to be two-dimensional in the whole parameter space. A more stable chaos seems to require that more voltage-dependent ion channels are added to the model. We thus added two artificial ion channels to the hippocampal model and found chaotic firing (Figure S2). A rather extensive search for chaotic firing in models with just one added ion channel gave no positive results.
The time evolution of the membrane potential (V) was calculated by solving the following equation (derived from Equation 1) numerically: (3) where INa and IK are functions of the activation parameters m and n, and the inactivation parameter h. Ileak is given by (4) INa and IK for the hippocampal and the myelinated axon model are described by the following expressions, based on the permeability concept of Goldman [26] and Hodgkin and Katz [27]: (5) and (6) where and denote the Na and K permeabilities when all Na and K channels are open, and thus represent the Na and K channel densities. R, T and F denote the gas constant, the thermodynamic temperature and the Faraday constant, respectively, and define. [Na]o, [Na]i, [K]o and [K]i are the external and internal concentrations of Na and K ions. The parameter values for the three models are listed in Table 3. For the squid axon model we use the original expressions by Hodgkin and Huxley [23] based on the conductance concept: (7) and (8) where and denote the Na and K conductances when all Na and K channels are open, thus representing Na and K channel densities. The activation and inactivation parameters (m, h and n) are in all three models described by their time derivatives: (9) (10) (11) where αi and βi denote rate functions. For the hippocampal neuron model they are defined as follows: (12) For the myelinated axon model the rate functions are defined as: (13) For the squid axon model the rate functions are defined as: (14) The stability analysis of the differential equations was performed as briefly described by Århem et al. [11]. The stationary potentials can be calculated with the expression for the gating parameters (m, n and h) at steady state together with Equation 3. The time derivates of the gating parameters are zero at stationary potentials, and hence the stationary values of the parameters (denoted m∞ etc) become: (15) (16) (17) Introducing these expressions into Equation 3, we obtain the following equation, the roots of which yield the stationary potentials (Vs): (18) This equation can be solved numerically and always yields at least one Vs. However, if the Na channel density (or) is large enough, the equation can for a defined stimulation interval give three equilibrium points, a requisite for the system to provide a saddle-node bifurcation (see Figure 3). We investigated the character of the equilibrium points, i. e. r* (V, m∞ (V), h∞ (V), n∞ (V) ), when the stimulation current (Istim) and the permeability or conductance parameters (and or and) representing the density of Na and K channels, were varied. This was done by linearizing the differential equations close to r* and by solving the characteristic equation (19) where denotes the identity matrix, and JM the Jacobian matrix (20) where, , and denote the time derivatives of the parameters. The solution to Equation 18 are the four eigenvalues λi (i = 1,2, 3 or 4) yielding an approximate time evolution of the system. Hence any perturbation δr around the equilibrium point r* can be written as (21) where ci (i = 1,2, 3 or 4) depends on initial conditions and ri (i = 1,2, 3 or 4) is the associated eigenvector. Two of the eigenvalues are in the present system (here called λ3 and λ4) always real and negative. Consequently the remaining two eigenvalues determine the character of the Vs (see Table 4); the two negative eigenvalues will cause its associated terms in Equation 20 to decay to zero. Hence, Equation 20 can be approximated as (22) If λ1 and λ2 are a complex conjugated pair (λ1,2 = a±bi), one can rewrite the equation, using Euler' s formula, as (23) why b/2π correlates with the firing frequency. All computations were done in custom software written in Mathematica 6. 0. 2 (Wolfram Research, Inc.) on a 64-bit IBM compatible computer. All values are given in SI-units.
All activity of the brain is manifested in electrical oscillatory patterns, shaped by the firing dynamics of the many neurons forming the brain networks. The underlying mechanisms of the firing pattern in the single neurons are still not fully understood. The distribution and identity of different channel types have been suggested as critical factors. We have suggested that the density of channels in the membrane is a fundamental complementary mechanism. In a hippocampal soma membrane model study we have shown that altering the ion channel densities can cause the membrane to switch between two qualitatively different firing patterns. Here we extend the analysis to two axon membranes. Unexpectedly, both show that channel density alterations do not cause switches between different firing behaviours. We believe that this is an important property of axon membranes, explaining their limited flexibility.
physiology/neuronal signaling mechanisms
biophysics/theory and simulation
physiology
biophysics
neuroscience
neuroscience/theoretical neuroscience
Ion Channel Density Regulates Switches between Regular and Fast Spiking in Soma but Not in Axons