|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: HZ CB PÅ. Performed the experiments: HZ CB PÅ. Analyzed the data: HZ CB PÅ. Contributed reagents/materials/analysis tools: HZ CB PÅ. Wrote the paper: HZ PÅ.
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.
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.
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 . 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 , , 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 – and successfully been used in analysing the coding properties of neurons –. To avoid confusion, and in accordance with the notation of Tateno and Robinson , 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 .
The intricate interactions between the many factors involved in the dynamical regulation of neuronal firing are poorly understood . The dominant idea is that different combinations of ion channel types explain the different patterns . In a previous study we proposed a complementary explanation , . We showed that both Type 1 and Type 2 behaviour can be simulated in a dynamical model of a hippocampal neuron  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 . An increasing number of studies suggest a role for physiological regulation of channel densities, even at a relatively short time scale –.
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 , . 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 , . 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 , ) 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) .
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  and the giant axon of Loligo forbesi . 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  to describe the dynamics of the squid giant axon. The hippocampal neuron model used is that developed by Johansson and Århem  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  to describe sciatic nerve fibres from the African clawed frog (Xenopus laevis) , . The giant axon model used is that of Hodgkin and Huxley , 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  and Hodgkin and Katz  (see Methods).
Istim is equal to the sum of the capacitive current (IC), charging the capacitor, and the ionic current Iionic. Thus,
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):
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 , , , , 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  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 ,  (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 –.)
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 ). 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 . 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  and the giant axon of the squid Loligo forbesi  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.
Figure 9 documents the calculated time-course of the changes in voltage calculated to occur in the model of the myelinated axon at different stimulation amplitudes for low and high densities of K channel (=0 and 40 µm/s). In this case, no combination of channel densities yielded Type 1 dynamics. However, the shape of the individual spikes differed at low and high K channel densities, an afterhyperpolarization being prominent only at a high K channel density. As shown in the figure, this model, in contrast to the hippocampal model, displays repetitive firing at =0, which may reflect the fact that myelinated axons at a later evolutionary stage, i.e. mammals, lack K channels and still fire repetitively .
A stability analysis revealed the mechanisms involved. Figure 10 displays Istim− diagrams, showing an Andronov-Hopf bifurcation line within the three-solution space. As summarised below, the resulting channel-density map shows a region with non-monotonic Iss-V relationship. However, unlike the narrow Andronov-Hopf region C1b of the hippocampal neuron model, the corresponding region of the myelinated axon model covers the whole C1 area. Consequently, saddle-node bifurcation dynamics is missing, explaining the absence of Type 1 dynamics in the myelinated axon model. It can also be noted that the double-limit cycle bifurcation occurs very close to the Andronov-Hopf bifurcation, why great care is needed to distinguish them numerically. That this kind of change can occur very close to each other in a parameter space is a known feature in these kinds of models , .
The findings suggest less dynamic flexibility in this axon membrane than in the hippocampal neuron model discussed above. Since the hippocampal neuron model is based on measurements of the soma membrane properties, the comparison between the myelinated axon and the hippocampal model mainly concerns a comparison between axonal and soma membranes. To get further information about the functional relevance of the found differences between axon and soma membranes we next analysed the classical squid giant axon membrane, using the description given by Hodgkin and Huxley .
Calculations showed that the giant squid axon model in similarity with the myelinated axon did not show saddle-node or Type 1 dynamics (Figure 11). In contrast to the myelinated axon model, however, the squid axon model could not fire repetitively at zero K channel density (=0), possibly related to the much higher leak conductance in the myelinated axon (see Table 2). The firing pattern at low values differed from that at higher values; it never showed damped oscillations.
The reason for the absence of saddle-node dynamics is evident in the Istim− diagram of Figure 12. As seen, there exists a three-root region, but in similarity with the myelinated axon model, there is an Andronov-Hopf bifurcation line within this area. Figure 13 shows the oscillation map for the squid axon model (B) in comparison with that of the myelinated axon (A). Figure 13 also shows the onset frequencies of the two axon models, being typically higher than 30 Hz for the squid axon (D) and higher than 100 Hz for the myelinated frog axon (C). In summary, the present analysis suggests a similar response pattern for models of different axon membranes in spite of considerable differences in mathematical structure. Furthermore, it suggests a functional difference between axon and soma membranes, the latter being more flexible.
The manner in which interactions between the ionic currents in a neuron determine the pattern and dynamics of firing is a multifaceted problem, having many ramifications within theoretical systems biology . In a previous study, we demonstrated that altering the ion channel densities in a four-dimensional dynamical model (comprising experimentally measurable parameters) of a hippocampal neuron could cause switches between Type 1 and Type 2 firing behaviour , . We also suggested that this channel density paradigm may explain the different threshold dynamics of regular and fast spiking cortical neurons , , , as well as that of Class 1 and Class 2 axons in Carcinus maenas . The physiological and pharmacological relevance of channel densities as bifurcation variable has recently been experimentally confirmed , .
In the present analysis we show that corresponding channel density alterations in two well-studied axon models (the amphibian myelinated axon and the squid giant axon) cannot change the firing dynamics, exclusively being of Type 2. This suggests that the hippocampal soma and the two axon membranes represent two distinct types of membrane with respect to the excitability pattern; one more flexible that can switch channel-density dependently between Type 1 and Type 2 dynamics (represented by the hippocampal neuron membrane and for simplicity here denoted M1/2) and one, less flexible, that exclusively shows Type 2 dynamics (represented by the membranes of the two axons and here denoted M2).
The mathematical background to the flexibility of the first membrane type is the existence of two types of bifurcations associated with separate regions in the − plane (regions B and C1). Type 1 dynamics (region C1) is associated with saddle-node on invariant cycle bifurcations (SNIC) and Type 2 dynamics (region B and A2) with double-limit cycle bifurcations (in some cases along with a subcritical Andronov-Hopf bifurcation). A requirement for the occurrence of a saddle-node bifurcation is the existence of three stationary voltages at near threshold stimulation and thus, a non-monotonic Iss-V relationship. In the present analysis we show that this requirement is not sufficient. In all three models we find regions in the − plane with a non-monotonic Iss-V relationship, that are associated with Type 2 dynamics (C1b regions). In the hippocampal neuron model the region consists of a narrow band, and in the axon models they cover the whole C1 region, and thus, the axonal models lack Type 1 threshold dynamics (for the analysis of the bifurcation structures in regions with non-monotonic Iss-V relationships of generic two-dimensional models, see –, ).
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, ; in dorsal horn neurons, ). 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 ). 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. , ). Such flexibility has also been shown in cortical fast spiking (Type 2) interneurons ; 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 .
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 , . 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 , , suggesting that that their cell membranes are of M1/2 type (for a computational analysis of modifying the dynamics of squid axons, see ). 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  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 , . 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 , , 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 , 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 . 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 , –. Interneurons with Type 2 dynamics have recently been shown to account for the cortical γ-oscillations (20–80 Hz) , which are considered to provide a temporal structure for information processing in the brain .
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 , ) 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 , 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 , 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:
where INa and IK are functions of the activation parameters m and n, and the inactivation parameter h. Ileak is given by
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  based on the conductance concept:
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:
where αi and βi denote rate functions. For the hippocampal neuron model they are defined as follows:
For the myelinated axon model the rate functions are defined as:
For the squid axon model the rate functions are defined as:
The stability analysis of the differential equations was performed as briefly described by Århem et al. . 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:
Introducing these expressions into Equation 3, we obtain the following equation, the roots of which yield the stationary potentials (Vs):
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
where denotes the identity matrix, and JM the Jacobian matrix
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
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
If λ1 and λ2 are a complex conjugated pair (λ1,2=a±bi), one can rewrite the equation, using Euler's formula, as
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.
Firing frequency as a function of Na channel density in a rat axon model. The firing frequency at threshold is plotted against Na permeability constant. For Na permeabilities less than 59µm/s the axon model does not show any repetitive firing. The curve is calculated from the equations described in Text S1, based on the voltage-clamp analysis by Brismar . Since this axon exclusively comprises Na and leak channels (it lacks K channels), the curve implies that this rat axon, in agreement with the Xenopus and the squid axon discussed in the main text, exclusively shows Type 2 dynamics for all channel densities. The Iss-V relationship is non-monotonic (N-shaped) for values of higher than 166 µm/s.
(0.54 MB EPS)
Chaotic firing in a modified hippocampal neuron model. Time course and phase-plot of a modified hippocampal neuron model, including an extra Na and an extra K channel. The figure depicts data from an extensive search for parameter values that give the model chaos like dynamics (i.e. aperiodic appearance over the time scale simulated). The extra channes were assumed to be described by the same kinetics as the hippocampal channels, i.e. by equations 4, 5 and 6 in the Methods section. The rate functions of the extra channels were calculated from equation 12 with the numerical values altered randomly ±50%. The permeability constants ( and in equations 5 and 6) were altered randomly in the range 0–100 µm/s, and the stimulating current in the range 0–1000 mA/m2. The search was performed using custom written software in Mathematica 6.0.2. Out of about 100.000 trials, eight parameter combinations yielded chaos like firing. When a similar search, equally extensive (100 000 trials), was performed with an hippocampal model comprising only one extra channel (Na or K), no positive results were found. The results seem to suggest that at least four channels (two Na and two K channels), and thus a seven-dimensional system, are required to give neuron models of the type studied here chaotic dynamics.
(1.63 MB EPS)
A model of the myelinated rat axon
(0.07 MB PDF)
We thank Professor John Guckenheimer, Cornell University, and Associate Professor Yuri A. Kuznetsov, Utrecht University, for valuable discussions.
The authors have declared that no competing interests exist.
This work was supported by The Swedish Medical Research Council (grant No 15083). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.