Application of dynamical systems methods to analyze the role of ionic currents on neuronal membrane resonance

Waves and oscillations in the nervous system are commonly produced in central pattern generators and are categorized through their frequency, shape, amplitude, and location. Central pattern generators are a group of neuronal networks that internally produce rhythmic oscillations and this rhythmogenesis does not require the presence of any external rhythmic stimuli. Chewing, breathing, walking, crawling and flying are some of these rhythmic behaviors which are generated in central pattern generators.

How rhythm generation occurs in a central pattern generator?

Over the past couple of decades, neurophysiologists have branched into two schools in answering this question: How rhythm generation occurs in a central pattern generator? Each school to answer this question utilizes the “half-cycle oscillator” model and “pacemaker-based” model. In the second model, these rhythmic oscillations may be generated by an individual pacemaker neuron or by a population of them. The other bursting neurons in central pattern generators commonly follow the rhythmic frequency of these pacemaker(s) neuron(s). For example, pyloric rhythm is generated in stomatogastric ganglion of crustaceans by a single pacemaker neuron (AB neuron). Other neurons in this central pattern generator, which is responsible for producing the pyloric rhythm, simply follow the frequency of AB neuron.

Resonance is the tendency of a neuron to produce amplified responses to oscillating stimuli if the stimuli have a certain range of frequency. This frequency is called preferred resonance frequency. There are two types of resonance: synaptic resonance and subthreshold membrane resonance. For a neuron to have synaptic resonance, the synapse requires to contain both facilitation and depression mechanisms in its response following a train of stimuli. In the case of subthreshold resonance, the stimulus is applied current and the neuronal response is membrane potential. Therefore, resonance is the maximal impedance. Electrical impedance extends the concept of resistance to AC circuits, describing not only the relative amplitudes of the voltage and current but also the relative phases.

(1) Iapp=I0 +ASin(2πft)

(2) V = V* + V1(f) sin (2π f t + θ(f))
(3) Z(f) = V1(f)/A

Both phase difference and the magnitude of the impedance are functions of the driving frequency.

As discussed by Dover et al (2006), synaptic and membrane resonance interact with each other to produce the combined neuronal resonance.

One of the most famous examples of resonance is in neurons of the Stomatogastric nervous system of crustaceans. In fact because of several advantages of this system, such as a minimal number of neurons in a CPG, relatively large size of neurons and easy identification of neurons and their known properties, this system has been one of the best well-studied CPGs. In this network, some neurons including pyloric dilator (PD) neurons exhibit membrane resonance properties.

One of the key features of oscillations in central pattern generators is their identical bursting frequency. For example, hippocampal theta waves have a range of bursting frequency between 4 Hz and 10 Hz. Many scientists now believe that resonance may be the key factor in governing the bursting frequency of neuronal oscillators. For example, Tohidi and Nadim (2009) demonstrated that the frequency of the pyloric rhythm in the stomatogastric nervous system of crustaceans is almost identical. This frequency is determined by the factors that shape the bursting behavior of the pacemaker neurons. Their results showed that the network frequency is correlated with membrane resonance of pacemaker neurons. In addition, Tseng and Nadim (2010) investigated the effects of other electrical properties of pacemaker neurons (i.e. voltage range and waveform), on the preferred resonance frequency and network cycle frequency. Interestingly, the network frequency and membrane resonance of pacemaker neurons are not solely sensitive to electrical characters of pacemaker neurons (such as some ionic currents), but also depend on voltage ranges and waveforms of oscillations.

What ionic currents are responsible for producing membrane resonance?

During the last decades several ionic currents are proposed to be responsible for the generation of resonance: persistent sodium (NaP), Hyperpolarization activated inward current (Ih) and Transient calcium current (ICaT) are among the ionic currents that may play the key role in the generation and/or amplification of membrane resonance.

As mentioned by Hutcheon and Yarom (2000), in order to produce resonance in a neuron, it is necessary to combine two sets of properties: 1) the passive properties of neuron which depends on the membrane conductance and the membrane resistant. 2) Active ionic currents that produce resonance. In fact, passive properties of neuron act as low- pass filter: it filters out inputs with high frequencies. In contrast, active properties of the neuron may act as a high-pass filter and it filters out inputs with low frequencies. The combination of these two properties may enable the neuron to produce resonance.

Two types of ionic currents play role in resonance: resonance currents and amplifying currents. Resonant currents are currents that by their activation and inactivation may produce resonance. Amplifying currents, amplify the resonance effects of resonant currents. Classically, Ih and K current are considered as resonant currents and NMDA current and persistent sodium currents are considered as amplifying currents. Also, Calcium may demonstrate both behaviors.

Our computational models of PD neuron:

Among the pyloric dilator (PD) neurons of stomatogastric nervous system resonance has been well studied. To clarify, we have previously proposed a computational model of the PD neuron, based on Hodgkin-Huxley equations with experimental-driven parameters. This model consists of Calcium current with activation and inactivation dynamics, Hyperpolarizition-activated inward current (Ih) and leak current. The model is defined by following equations:

(1) CdV/dt=–gCa mh(V–ECa)– ghr(V–Eh)– gL(V–EL)+Iapp (2) Iapp=I0 +ASin(2πft)
(3) dx / dt = ( x∞ (V) – x) /Taux
(4) x∞(V)=1/(1+exp((V–Vx)/kx)

As been expected, this model demonstrated the possibility to produce resonance based on three mentioned ionic currents. Plus in the computational model, many features of resonance (i.e. resonance frequency and impedance profile) are compatible with experimental data acquired from in-vitro recordings of PD neurons in Crab Cancer Borealis in Prof. Nadim’s laboratory.

Impedance profile of three- dimensional computational model of PD neuron
Figure 1- Impedance profile of three- dimensional computational model of PD neuron demonstrate resonance around 1 HZ. Frequencies are ranging from 0.4 HZ to 2.6 HZ. Note that the frequency at which the peak of the impedance occurs corresponds to the preferred resonance frequency.

Because our model has calcium and H current and both of these currents are resonant currents, we searched if we can have resonance with either of these currents. Our results led to two computational models. The first model only contains calcium and leak currents. In this model calcium have both activation and inactivation. The second model has H current and leak. Our results demonstrated that both models produce resonance (figure 3). Interestingly, the second model demonstrated that for producing resonance effects, the presence of leak current and h current is enough. Therefore, we developed a computational model of a PD neuron, based on these two currents.

impedance profile of two-dimensional model with leak and H current demonstrates significant resonance

Figure 2- impedance profile of two-dimensional model with leak and H current demonstrates significant resonance.

Interestingly, the second model consists of two sets of variables: V(t) and r(V). Both variables are dependent on each other and more importantly; we may apply dynamical systems tools to analyze the behavior of the systems. Because our reduced model, only has two ionic currents and consequently two variables (Voltage and activation of H current), we can apply dynamical systems tools more efficiently to understand the role that each variable plays in order to generate membrane resonance.

Z Impedence (Mega Ohm)

Reduction of our model to two variables:

As discussed earlier, we may also observe resonance phenomena in a single neuron containing only Ih and leak current (IL). Therefore, we reduced our model to a two variable model, as following:

(1) CdV/dt=–ghr(V–Eh)–gL(V–EL)+Iapp (2) dr / dt = ( r∞ (V) – r) /Taur(V)

(3) r∞ (V) = 1 / (1 + exp((V + 65)/4)
(4) Taur(V) = 500+ 400/ (1 + exp (-(V + 70)/5)

The parameters of the system are: A=0.05 μA; gl=0.005 mS; gh=0.019 mS; Eh=-20 mV; EL=-60 mV; Cm=1 μF;

Iapp=I0 +ASin(2πft)

The effects of conductance of H current on resonance in the two dimensional reduced model

Figure 3 – The effects of conductance of H current on resonance in the two dimensional reduced model. The value of leak conductance is constant. Units of conductances are mS. As gH increases, preferred resonance frequency shifts to higher values.


The above Hodgkin-Huxley-based system of equations is nonlinear. In order to better understand the behavior of the above system, we linearize it and do linear approximation to functions about the fixed point (V ̄, r ̄). The linear version of above system comes as following:

(*) CdV/dt=–Fv(V–V ̄)–Fr(r–r ̄)+Iapp (**) dr / dt = – Gv (V– V ̄) – Gr (r– r ̄)

Here, (V ̄, r ̄) correspond to the fixed point, where dV = dr =0 ; Fv, Fr, Gv, and Gr , correspond to partial derivatives of (*) and (**) with respect to V, r, V, r respectively.

Fv= –gL–ghr∞ (V ̄) Fr= – gh (V ̄– Eh) Gv= r`∞ (V ̄)/ Taur(V ̄) Gr= –1/(V ̄)

comparing resonance between linear and nonlinear models
Figure 4- comparing resonance between linear and nonlinear models. Linear approximation is more beneficial when the amplitude of applied current is small. Therefore in small amplitudes we may use linearized system instead of nonlinear system.



The nullclines of above nonlinear system are:

V-null: dV/dt=0 r = ( – gL (V– EL) + Iapp )/gh (V– Eh) r-null: dr/dt=0 r= r∞ (V)

Because the applied current (Iapp) is sinusoidal and is time dependent, the V-nullcline is time dependent and varies by changes in applied current. Therefore V-nullcline is moving. In contrast, r-nullcline is fixed. Therefore this system contains a moving nullcline and a fixed nullcline.

Phase plane analysis:

We applied two-dimensional phase plane analyses and draw the nullclines in various frequencies. The green nullcline follows r= r∞ (V) and is not moving. In contrast because Iapp is time dependent and is changing, the red nullcline is moving. Therefore we have a fixed nullcline and a moving nullcline and their relative position governs the path that trajectories follow (fig. 5). Because the derivative of V is zero in V-nullcline and similarly derivative of r is zero in r-nullcline, movement of trajectories in nullclines are vertical and horizontal respectively. This shapes the direction and force of moving trajectories in phase-plane of differential equations.

Trajectories follow the moving nullcline with partial freedom and try to reach to the ultimate in V-axis. This ultimate in V-axis determines the voltage range that a trajectory may visit and consequently determines the impedance of this system of differential equations. Increasing the frequency of applied current causes increase in the speed of the moving nullcline. But in higher speeds, the partial freedom of follower trajectories deprives the trajectories from reaching to “ultimate voltages” which they visited in lower speeds (fig. 6). Consequently in higher frequencies impedance decreases and this’d be the demonstration of low-pass filtering in the language of dynamical systems.

The flow and pattern of trajectories in the two dimensional model

Figure 5- The flow and pattern of trajectories in the two dimensional model. Here, V- nullcline (red line) moves in horizontal direction. The r-nullcline (green line) is fixed. Flow of trajectories in green line region is horizontal, but in red line trajectories move vertically. The junction of two nullclines is the fixed point. Trajectories follow the pattern of flow which is defined by interaction of two nullclines. The frequency of applied current in this case is relatively low (0.1 Hz) and the resulting limit cycle has a boomerang shape.


Figure 6- The effects of changing the frequency of applied current on the reduced two- dimensional dynamical system. Applied current follows the equation: Iapp= A Sin (2 π f t); A=0.1 μA; In above graphs frequencies are as following: 0.1 Hz, 0.2 Hz, 0.4 Hz, 0.7 Hz, 1.2 Hz, 1.7 Hz, 2.5 Hz and 4 Hz. Note that as frequency increases, initially the range of voltage increases (equal to increase in impedance). After a certain frequency (preferred resonance frequency) the range of voltage decreases (equal to decrease in impedance).

This partial freedom also plays the key role in high-pass filter, which is essential in generation of resonance. As demonstrated in figure 7, as frequency increases from the very low-frequencies, the shape of the limit cycle changes. In this example, the limit cycle changes from a boomerang shape to a horizontal ellipse and later. The ellipse shrinks. In fact, in the ranges of frequency that transforms a boomerang to horizontal ellipse, changes in the shape of limit cycle causes increases in the voltage range and subsequently, increase in impedance. On the other hand, once ellipse starts to shrink, impedance decreases. Overall, these trends of increase and decrease in voltage range and impedance generates resonance.

Figure 7- Overlap of different limit cycles with different frequencies of applied current. Changes in frequency of applied current alter the speed of moving nullcline in the phase plane. Consequently the pattern of flow would be changed faster. This induces different destinies for the moving trajectories in the phase plane. As is shown in this figure, the voltage ranges in movement of trajectories is different in various frequencies and initially it increases and then decreases. This is exhibition of resonance in phase plane. X-axis is voltage or V and Y-axis is activation of h current or r. The frequency of applied current varies from 0.1 HZ to 4 HZ.


In conclusion, the resonance behavior of the system is based on both high-pass filter and low-pass filter. These filters are generated by the partial freedom in the movement of trajectories and the relative speed of moving trajectories to moving nullcline in the phase- plane. The parameters of the system that affect this relative speed and partial freedom alter resonance of a neuron.

One important question here is the difference between resonant systems and non-resonant systems. Like resonant systems, non-resonant systems may have a similar pair of fixed and moving nullclines. But the combination of the partial freedom of trajectories and the relative speeds of trajectories to moving nullcline do not form any changes in the shape of the limit cycle in low frequencies. Therefore, these systems only act as low pass filter and impedance decreases as the frequency of applied current increases.

Overall, we may benefit from dynamical systems tools to differentiate between the resonant systems and non-resonant systems. Also, these useful tools may predict an estimation of the preferred resonant frequency.


  1. Izhikevich E., Desai N.S., Walcott E., Hoppensteadt F.C. (2003), Bursts as a unit of neuronal information: selective communication via resonance, Trends in neuroscience, 26:161-167
  2. Tohidi V., Nadim F. (2009), Membrane resonance in bursting pacemaker neurons of an oscillatory network is correlated with network frequency, The journal of neuroscience, 29(20):6427-6435
  3. Tseng H., Nadim F. (2010), The membrane potential waveform of bursting pacemaker neurons is a predictor of their preferred frequency and the network cycle frequency , The journal of neuroscience, 30(32):10809-10819
  4. Cell response to periodic inputs: Horacio Rotstein notes
  5. Richardson M.J.E., Brunel N., Hakim V. (2003), From subthreshold to firing-rateresonance, J Neurophysiol 89:2538-2554
  6. Ulrich D. (2002), Dendritic resonance in rat neocortical pyramidal cells, JNeurophysiol, 87:2753-2759
  7. Mascagni M. (1990), The backward Euler method for numerical solution of the Hodgkin-Huxley equations of nerve conduction, Siam J Numer. Anal., 27 (4):941-962
  8. Koch C., (1984), Cable theory in neurons with active, linearized membranes,Biol Cybern. 50,15-33
  9. Schreiber S., Erchova I., Heinemann U., Hertz A.V.M. (2004), Subthresholdresonance explains the frequency-dependent integration of periodic as well asrandom stimuli in the Entorhinal cortex, J Neurophysiol, 92:408-415
  10. Zhao S., Golowasch J., Nadim F., (2010), Pacemaker neuron and network oscillations depend on a neuromodulator-regulated linear current, Frontiers in behavioral neuroscience, 4:1-9

By: Nima Schei
Advisors: Prof. Nadim and Prof. Rotstein January 2011

To download the PDF version of this paper, click here.