Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2906914

Formats

Article sections

Authors

Related links

SIAM J Appl Dyn Syst. Author manuscript; available in PMC 2010 July 20.

Published in final edited form as:

SIAM J Appl Dyn Syst. 2009 January 1; 8(4): 1564–1590.

doi: 10.1137/090760994PMCID: PMC2906914

NIHMSID: NIHMS208196

Yu Zhang, Department of Mathematical Sciences, New Jersey Institute of Technology Newark, NJ 07102 ; Email: ude.tijn@gnahz.uy;

Current Address: Department of Pharmacology and Systems Therapeutics, Mount Sinai School of Medicine, New York, NY 10029

See other articles in PMC that cite the published article.

The transient potassium A-current is present in almost all neurons and plays an essential role in determining the timing and frequency of action potential generation. We use a three-variable mathematical model to examine the role of the A-current in a rhythmic inhibitory network, as is common in central pattern generation. We focus on a feed-forward architecture consisting of an oscillator neuron inhibiting a follower neuron. We use separation of time scales to demonstrate that the trajectory of the follower neuron within each cycle can be tracked by analyzing the dynamics on a 2-dimensional slow manifold that as determined by the two slow model variables: the recovery variable and the inactivation of the A-current. The steady-state trajectory, however, requires tracking the slow variables across multiple cycles. We show that tracking the slow variables, under simplifying assumptions, leads to a one-dimensional map of the unit interval with at most a single discontinuity depending on *g _{A}*, the maximal conductance of the A-current, or other model parameters. We demonstrate that, as the value of

In many oscillatory neural networks, the activity time of neurons within each cycle plays an important role in determining the network output [23, 24, 17]. It has been shown that the A-current, a transient outward potassium current, is present in many neurons and is essential in determining the frequency-current response of neurons [5, 11, 12, 28] and important in setting the timing of activity [8, 26]. In neurons that receive inhibitory input, the A-current can play an essential role in determining the response. For example, in Central Pattern Generator (CPG) networks where inhibition is predominant, the A-current acts to delay the onset of action potentials, thus contributing to the activity time of the neuron [10]. Other factors such as the synaptic strength and other intrinsic properties also play significant roles in affecting the activity of neurons in inhibitory networks [1, 16].

We investigate how the A-current affects the activity of a neuron that receives rhythmic inhibition. Our model is motivated by the activity of the follower neurons in the pyloric CPG in the stomatogastric nervous system of the crab *Cancer borealis*. The follower PY neurons, for example, receive rhythmic inhibition from the pyloric pacemaker ensemble and produce bursts of action potentials following the inhibition. As with many CPGs, the pyloric oscillations have little variability and the PY neurons burst at a relatively constant phase within each cycle of oscillation (Fig. 1(a)). In a previous study we explored how the interaction between the A-current and other intrinsic properties of the follower neuron determine the post-inhibition activity phase of this neuron [29].

In this study we show that the presence of the A-current in a neuron receiving periodic inhibitory input from an oscillator may lead to complex and non-intuitive dynamics. Our model consists of an oscillator neuron that inhibits a follower neuron. For simplicity, we model the oscillator using a square-wave formalism. The model of the follower neuron is based on the two-dimensional Morris-Lecar model which is characterized by a cubic voltage (*v*) nullcline and a recovery variable (*w*) leading to oscillatory or excitable behavior [21, 25]. The addition of the A-current to the model introduces an additional dynamic variable (*h*) which results in an additional stable “middle” branch to the *v*-nullcline in the *v*-*w* phase plane [2, 29]. In the three-dimensional phase space, this middle branch forms a surface. Using separation of time scales, it can be shown that the flow on this surface is controlled by the slow variables *w* and *h*. We show that the effect of the A-current on the activity of the follower neuron is determined by the geometry and kinetics on this *w*-*h* slow manifold. In each cycle of oscillation, the A-current can delay the transition of the follower neuron to the active state, make it to return to the inactive state, or force it to remain on the middle branch. There exist two cases that correspond to distinct geometric structures on the *w*-*h* slow manifold. These structures, together with the *w* and *h* time constants, determine which of these possibilities occur.

For various sets of model parameters, we use the cycle-by-cycle analysis on the slow manifold to derive a map of the unit interval whose dynamics capture the behavior of the full set of equations. This map has a discontinuity which is dependent on the maximal conductance *g _{A}* of the

The outline of the paper is as follows. In Section 2, we present the model equations for the follower neuron and discuss how to use geometric singular perturbation theory to derive sets of reduced equations. The *w*-*h* slow manifold and equations governing the flow on it are defined. In Section 3, we present our main results. First, we examine the transient response of the follower neuron following the inhibitory input in one cycle of oscillation. The main emphasis here is to show how to use the *w*-*h* slow manifold to determine the behavior of the trajectory and whether it ultimately moves to the active state, returns to the silent state, or remains on the middle branch. We then consider the steady state response of the follower neuron and derive a one-dimensional map that captures the behavior of this neuron. By analyzing the dynamics of this map, we prove the existence of various types of periodic solutions. The paper concludes with a brief Discussion.

Our three-variable model is meant to capture the envelope of the bursting oscillations of the PY neuron when receiving rhythmic inhibition from the pacemaker neurons. In our model, the pacemaker neuron is represented as a single oscillator (O) and the PY neuron as the follower (F; see Fig. 1(b)). Since our interest is in the activity of the follower neuron, we simply define the membrane potential *v _{O}* of O as a square wave oscillating periodically between −50 mV and 0 mV with active duration

$$\begin{array}{c}\epsilon \frac{dv}{dt}=f(v,w)-{I}_{A}-{I}_{\mathit{\text{inh}}}\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{w}(v)}\hfill \\ \frac{dh}{dt}=\frac{{h}_{\infty}(v)-h}{{\tau}_{h}(v)}\hfill \end{array}$$

(1)

Where

$$f(v,w)={I}_{\mathit{\text{ext}}}-{g}_{L}[v-{E}_{L}]-{g}_{Ca}{m}_{\infty}(v)[v-{E}_{Ca}]-{g}_{K}w[v-{E}_{K}]$$

(2)

represents the ML terms and

$${I}_{A}={g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]$$

(3)

denotes the A-current. ε is a non-dimensional parameter and is taken as a small positive number. *w*_{∞}(*v*), *n*_{∞}(*v*) and *h*_{∞}(*v*) are sigmoidal functions representing the steady-state values of the activation variable for the potassium current, the activation and inactivation variables for the A-current respectively. Note that the activation of the A-current is assumed to be instantaneous function of the membrane potential. Each sigmoidal function has the form:

$${x}_{\infty}(v)=\frac{1}{1+\text{exp}(\frac{v-{v}_{x}}{{k}_{x}})}$$

(4)

Here *k _{x}* is negative for activation variables and positive for inactivation variables. The time constants τ

The oscillation cycle of the model follower neuron F shown in the membrane potential trace and in the *v*-*w* phase plane. **(a)** The membrane potentials of F and the oscillator O (modeled as a simple square wave). *T*_{act} and *T*_{in} represent the active and inactive **...**

The synaptic current from O to F is given by

$${I}_{\mathit{\text{inh}}}={g}_{inh}{s}_{\infty}({v}_{O})[v-{E}_{inh}]$$

(5)

where *s*_{∞}(*v _{O}*) is set as a Heaviside function with threshold of −25 mV lying between O’s resting and active values. As a result, both the onset and decay of inhibition occur quickly.

Although we are dealing with a three-variable model, we will keep track of the model behavior by examining the dynamics in the family of *v*-*w* phase planes in which the variable *h* is treated as a parameter. This enables us to compare the dynamics of our model with the behavior of the standard ML model. Later we extend our analysis by keeping track of the trajectories in the *w*-*h* phase plane and the full phase space. The steady state values of the original ML model are described by a cubic *v*-nullcline and a sigmoidal *w*-nullcline in the *v*-*w* phase plane. In this study we assume that, in the absence of inhibition, there is a stable fixed point on the right branch of the *v*-nullcline. This corresponds to the activity of the follower PY neurons (Fig. 1(a)) which, in the absence of inhibition, remain active (not shown). The presence of the A-current term creates a new, stable “middle branch” on the *v*-nullcline and therefore its shape becomes quintic (Fig. 2(b)). For convenience, we call the left, middle and right descending branches LB, MB and RB respectively. The inhibition from the oscillator neuron O causes the *v*-nullcline to move down a certain distance depending on the strength of the inhibitory synapse. We shall refer to the corresponding branches of the fully inhibited *v*-nullcline as LB^{1}, MB^{1} and RB^{1}, respectively (Fig. 2(b)). The sigmoidal *w*-nullcline may intersect the *v*-nullcline and the inhibited *v*-nullcline on any of these six descending branches depending on the shape and positions of these nullclines. It is straightforward to show that every intersection point on any of these branches yields a stable fixed point [25].

The membrane potential of the follower neuron can be mapped to a trajectory curve in the *v*-*w* phase plane, as shown in Fig. 2(c). For ε (in Eq. (1)) small enough, the system is singular perturbed [20]: in some regions of the phase space, *v* changes very quickly while *w* and *h* remain nearly constant. In other regions, the behavior of *v* can be slaved to that of *w* and *h*. Equations to describe either can be obtained by setting ε = 0 in Equation (1) or a time-rescaled version of Equation (1). These two sets of equations are respectively referred to as the slow and the fast equations and can be obtained as follows. The slow equations are found by setting ε = 0 in Equation (1):

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]-{I}_{inh}\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{w}(v)}\hfill \\ \frac{dh}{dt}=\frac{{h}_{\infty}(v)-h}{{\tau}_{h}(v)}\hfill \end{array}$$

(6)

The fast equations are obtained by rescaling *t* = εξ in Equation (1) and then setting ε = 0:

$$\begin{array}{c}\frac{dv}{d\xi}=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]-{I}_{inh}=F(v,w,h)\hfill \\ \frac{dw}{d\xi}=0\hfill \\ \frac{dh}{d\xi}=0\hfill \end{array}$$

(7)

Based on the facts above, the solution of the three-variable model can be projected as a trajectory in the *v*-*w* phase plane (Fig. 2(c)). In the *v*-*w* phase plane, the intersection points of the *w*- and *v*-nullclines on LB^{1} and RB play an important role in influencing the solution trajectory. For fixed values of *h*, these are stable fixed points in the *v*-*w* phase plane and we refer to these points, respectively, as FP_{1} and FP_{2} (Fig. 2(c)). Similarly, depending on the shape and position of the *w*-nullcline, there may be a stable fixed point on MB or MB^{1} (Figs. 2(c) and 3(b)). When this fixed point exists, we refer to it as FP. If the *w*-nullcline is above MB, there is no fixed point on MB. In this case, the flow on MB is still attracted toward the *w*-nullcline as explained in detail in the Results.

The intrinsic parameters can determine the fate of the follower neuron trajectory. Two qualitatively distinct cases are shown in the phase space based on the shape and relative position of the nullclines which are in turn determined by the relative positions **...**

The slow equations force the trajectory to move along the descending branches of the *v*-nullcline toward *w*_{∞}(*v*) (stable fixed points of Eq. (6)) while the fast equations control the trajectory during the jumps between the branches of the *v*-nullcline. A “singular trajectory” can be pieced together from the solutions of the slow equations (Eq. (6)) that determine the movement on the branches and the fast equations (Eq. (7)) that determine the jumps between these branches. It can be shown that for ε small enough, a solution to the full system (Eq. (1)) exists in a O(ε)-neighborhood of the singular trajectory [20]. An example of the pieced-together trajectory can be seen in Fig. 2(c). In this case, the trajectory obeys the slow equations and moves toward FP_{1} on LB^{1} during *T _{act}*; it stays near FP

We make several assumptions to ease the analysis. However, none of the following assumptions change the qualitative nature of the results. First we assume that *w*_{∞}(*v*) = *w _{FP1}* along LB

On LB:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{FP1}-w}{{\tau}_{wl}}\hfill \\ \frac{dh}{dt}=\frac{1-h}{{\tau}_{hl}}\hfill \end{array}$$

(8)

On MB:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{wm}}\hfill \\ \frac{dh}{dt}=\frac{-h}{{\tau}_{hm}}\hfill \end{array}$$

(9)

On RB:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{wh}}\hfill \\ \frac{dh}{dt}=\frac{-h}{{\tau}_{hh}}\hfill \end{array}$$

(10)

On LB^{1}:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]-{g}_{inh}[v-{E}_{inh}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{FP1}-w}{{\tau}_{wl}}\hfill \\ \frac{dh}{dt}=\frac{1-h}{{\tau}_{hl}}\hfill \end{array}$$

(11)

On MB^{1}:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]-{g}_{inh}[v-{E}_{inh}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{wm}}\hfill \\ \frac{dh}{dt}=\frac{-h}{{\tau}_{hm}}\hfill \end{array}$$

(12)

On RB^{1}:

$$\begin{array}{c}0=f(v,w)-{g}_{A}{n}_{\infty}(v)h[v-{E}_{K}]-{g}_{inh}[v-{E}_{inh}]\hfill \\ \frac{dw}{dt}=\frac{{w}_{\infty}(v)-w}{{\tau}_{wh}}\hfill \\ \frac{dh}{dt}=\frac{-h}{{\tau}_{hh}}\hfill \end{array}$$

(13)

When a trajectory leaves LB^{1}, the trajectory may land on MB or RB. It will land on MB if *g _{A}* is sufficiently large and τ

All time constants are assumed to be of the same order of magnitude as *T _{act}* (O(

In a previous study we showed the effect of the potassium A-current in affecting the activity phase of a follower neuron receiving rhythmic inhibition [29]. In that study, we analyzed the role of the A-current with the simplifying assumption that the steady-state activation curve *n*_{∞}(*v*) is a steep sigmoid, i.e. as *k _{n}*→0. The consequence of this assumption in the

We first analyze the fate of the trajectory following a single cycle of inhibition. This fate is determined by the dynamics on the middle branch (MB) of the *v*-nullcline in the *v*-*w* phase plane. Assume that at *t* = 0 the follower neuron is inhibited and its trajectory lands on LB^{1}. During the duration of inhibition *T _{act}*, the trajectory moves downwards along LB

Figure 3 shows two possible situations related to the shape and relative position of the *v*-and *w*- nullclines. In Case I the *w*- and *v*- nullclines do not intersect on MB and *dw/dt* > 0 along this branch (Fig. 3(a)), while in Case II, the *w*-nullcline crosses the *v*-nullcline at FP on MB (Fig. 3(b)). The trajectory cannot reach UK in Case II therefore it can only jump to the right if LK moves up fast enough during *T _{in}*. In Case I, however, the

Understanding the shrinking of MB in the *v*-*w* phase plane is difficult; therefore we consider the dynamics in the *w*-*h* phase plane, where the curves corresponding to LK, UK and FP, as well as the trajectory on MB can be properly visualized. Figures 3(c) and 3(d) show the *w*-*h* phase plane corresponding to the *v*-*w* phase planes of Figs. 3(a) and 3(b), respectively. The point denoted “start” is the location at which the trajectory lands on MB and has coordinates (*w _{FP1}, h_{start}*). The shaded region in each figure is a positively invariant region and shows the points the trajectory can access in this phase plane. The fact that the curve FP (corresponding to the fixed point on MB) is to the left of UK in Fig. 3(d), for example, means that the trajectory can never reach UK. This contrasts with Fig. 3(c) in which the trajectory may potentially reach UK. The main advantage of analyzing the dynamics in the

In general, within one cycle of the oscillation, once the trajectory is on MB, three possible situations can arise. 1. the trajectory reaches LK and jumps to the right branch (RB); 2. the trajectory reaches UK and jumps back to the left branch (LB); or 3. the trajectory remains on MB during the cycle. The fate of the trajectory can be affected by a number of factors, including the time constants of the *w* and *h* variables, the inactive (*T _{in}*) and active (

In Case I, there is no fixed point on MB. Thus it always holds that *dw/dt* > 0 on MB. Case I is summarized in Figs. 4(a) and 4(b). In this case UK lies above LK until *h* is very close to 0, when they meet and disappear. The trajectory can reach UK and then jump to the left if it moves fast enough in the *w* direction (τ_{wm} is small enough), as in Subcases 1 and 3 (Fig. 4(b)). In contrast, if the trajectory moves fast enough in the *h* direction (τ_{hm} is small enough), it can reach LK and jump to the right, as in Subcases 2 and 4. Note in Subcases 1 and 2, the trajectory spends very little time on MB, in contrast to Subcases 3 and 4 where the trajectory is delayed along MB due to the A-current. For Subcases 1–4, using Fenichel theory [7], it is possible to use the singular flow with ε = 0 to approximate the actual flow for ε > 0 in a neighborhood of a normally hyperbolic manifold. In Subcases 1–4, although normal hyperbolicity is lost at LK or UK, there is no ambiguity in the direction which the fast flow should carry the trajectory, once it leaves MB. However, there is a unique set of time constants (τ*_{wm}, τ*_{hm}) where τ*_{wm} ~ τ*_{hm}, for which the singular trajectory starting at (*w _{FP1}, h_{start}*) reaches the intersection of UK and LK (Subcase 5). Here, there is no way to uniquely determine which direction the fast flow should carry the trajectory away from MB. Thus the ε > 0 fate of this trajectory in Subcase 5 can not be discerned from the ε = 0 singular solution. This will also be true when the time constants lie in a small neighborhood of (τ*

There exist three subcases in which FP exists on MB in the *v*-*w* phase plane (Fig. 5(a), (c) and (e)). Case II.a occurs when FP intersects LK (Fig. 5(a) and (b)), Case II.b when FP intersects UK (Fig 5(c) and (d)), and Case II.c when LK, UK and FP all intersect at a common point. Consider first Case II.a. The region between LK and FP is positively invariant and the trajectory initially starts on MB within this region. Thus the trajectory must eventually leave MB through LK and jump to the active state. A primary question is how the time spent on the MB depends on the two relevant time constants τ_{wm} and τ_{hm}. As shown in Fig. 5(b), when τ_{wm} << τ_{hm}, the trajectory moves rapidly in the *w*-*h* phase plane toward FP in the *w* direction and then tracks the FP curve down until it reaches LK (Subcase 1). The time spent on MB mostly depends on τ_{hm}. In contrast, when τ_{wm} >> τ_{hm}, the trajectory moves vertically down in the *h* direction until it reaches LK with little time spent on MB (Subcase 2). When τ_{wm} and τ_{hm} are of the same order, the trajectory may reach a neighborhood of FP first (Subcase 3), or reach LK directly (Subcase 4). In both Subcases, the A-current delays the trajectory from reaching the active state.

The dynamics of Case II in the *v*-*w* and *w*-*h* phase planes. **(a)** In Case II.a, the trace of FP intersects LK as *h* decreases in the *v*-*w* phase plane. **(b)** In the *w*-*h* phase plane, four subcases are possible for Case II.a based on the relative size of the two **...**

In Case II.b, FP intersects UK (Fig. 5 (c) and (d)). Now the positively invariant region contains a portion of UK. In Subcases 1 and 3, the curve FP channels the trajectory toward UK whereas, in Subcase 5, the trajectory directly reaches UK. For these Subcases, the singular trajectory reaches UK and thus returns to the silent state. Subcases 2 and 4 show the trajectories that reach LK and leave MB by jumping to the active state. Subcase 6 shows the case when the trajectory reaches the intersection of LK and UK (with some set of time constants (τ*_{wm}, τ*_{hm})). As before (Case I, Subcase 5), the singular flow cannot be used to determine the ε > 0 fate of this trajectory. Again this will also be true for a small set of time constants in a neighborhood of (τ*_{wm}, τ*_{hm}).

Case II.c lies between Cases II.a and II.b. Now LK, UK and FP meet at a single point. Since the vector field points down along FP, the singular flow cannot access UK and must leave MB through LK. However in Subcases 1, 3 and 5, the trajectory is channeled to the common intersection point and, as before, the ε > 0 fate of this trajectory cannot be determined by the singular flow. However, unlike the rare cases shown in Case I or II.b which occur for a small neighborhood of (τ*_{wm}, τ*_{hm}), the ambiguity in this case exists for a much larger set of time constants. Suppose again that Subcase 5 exists for (τ*_{wm}, τ*_{hm}). Then for any τ_{wm} < τ*_{wm}, the trajectory will be funneled to a neighborhood of FP and the ambiguity will arise. Figure 5(e) shows one realization of this effect where the voltage trace of F is shown in which the solution sometimes jumps to the active state and other times returns to the silent state.

We now turn our attention to showing the existence of periodic solutions by tracking the trajectory of the follower neuron over cycles. To do so, we use a set of recursively defined equation that provided mathematical expression for the time *t _{m}* that the follower neuron spends on the middle branch and, in our previous study [29], were used to give a description of 1-1 periodic activity of the follower neuron with the oscillator. We now use these equations to describe a 1-dimensional map of the unit interval whose dynamics directly determine the activity of F. In particular, we show that this map has a rich bifurcation structure that gives rise to different types of

Once the trajectory is on the middle branch, three distinct possibilities arise: 1. The trajectory can reach LK or UK prior to the onset of the next cycle of inhibition. 2. The trajectory does not reach either curve prior to this time and, if the inhibition is sufficiently weak, it remains on the middle branch for at least one more cycle. 3. The trajectory does not reach LK or UK and, if the inhibition is strong, it returns to the left branch. In this latter case, which is not explored any further, the follower neuron is never active and simply produces subthreshold oscillations out of phase with the oscillator.

Figure 6 shows a trajectory and the corresponding nullclines for a 2:1 periodic case which can be obtained if *g _{A}* or, alternatively, τ

The map that we derive is for the inactivation variable *h*. This map describes the behavior of F in the first two cases described in 3.2.1 and also includes the possibility that the trajectory does not land on the middle branch. Let *h ^{n}* denote the value of

$${h}^{n}={h}^{n-1}\text{exp}(-\frac{P}{{\tau}_{hm}})$$

(14)

The values of *t _{m}^{n}* are calculated by setting the right-hand side of the first equation in Eq. (1) equal to zero:

$$f({v}_{\theta},{w}_{FP})-{g}_{A}{h}^{n-1}\text{exp}(-\frac{{{t}_{m}}^{n}}{{\tau}_{hm}})[{v}_{\theta}-{E}_{K}]=0$$

(15)

using the assumption that the middle branch is vertical with *v=v*_{θ} and the assumption that τ_{wm} is small (*w* reaches a neighborhood of *w _{FP}*). The more general case is addressed in the Discussion. The map Π describing the value

$$\begin{array}{c}{h}^{n}=({h}^{n-1})=\{\begin{array}{ccc}1+[{h}^{n-1}\text{exp}(-\frac{{T}_{\mathit{\text{in}}}}{{\tau}_{hh}}+(\frac{1}{{\tau}_{hh}}-\frac{1}{{\tau}_{hm}}){{t}_{m}}^{n})-1]\text{exp}(-\frac{{T}_{\mathit{\text{act}}}}{{\tau}_{hl}})\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}{{t}_{m}}^{n}{T}_{\mathit{\text{in}}}\hfill & (\mathrm{a})\hfill \\ {h}^{n-1}\text{exp}(-\frac{P}{{\tau}_{hm}})\hfill & \text{if}\phantom{\rule{thinmathspace}{0ex}}{{t}_{m}}^{n}{T}_{\mathit{\text{in}}}\hfill & (\mathrm{b})\hfill \end{array}\hfill & \begin{array}{cc}{{t}_{m}}^{n}={\tau}_{hm}\phantom{\rule{thinmathspace}{0ex}}\text{ln}\frac{{g}_{A}{h}^{n-1}[{v}_{\theta}-{E}_{K}]}{f({v}_{\theta},{w}_{FP})}\hfill & \hfill \text{(c)}\end{array}\hfill \end{array}$$

(16)

If *t _{m}^{n}* <

Using XPP we computed a bifurcation diagram for Π as a function of *g _{A}* (Fig 7). For small values of

Bifurcation diagram of map Π as the parameter *g*_{A} is varied. **(a)** Steady-state values of *h*^{n} indicate a period-adding process as *g*_{A} is increased. Solutions go from a stable fixed point to period 2, 3, 4, … marked as regions 1–4 in **...**

Numerical simulation of Equation (1) (panels a1-d1) and the analytic solutions of the map Π (panels a2-d2 with arbitrary initial condition of *h*^{n}=0.1) for different values of *g*_{A} corresponding to the points marked as a–d in Fig. 7A. **(a)–(c)** **...**

To explain how the map Π captures the complicated dynamics of the trajectory, we consider how the discontinuity in Π depends on *g _{A}*. Let

$${T}_{in}={\tau}_{hm}\phantom{\rule{thinmathspace}{0ex}}\text{ln}\frac{\hat{{g}_{A}}[{v}_{\theta}-{E}_{K}]}{f({v}_{\theta},{w}_{FP})}$$

(17)

which is Equation (16)c with *t _{m}^{n}=T_{in}* and

$${T}_{in}={\tau}_{hm}\phantom{\rule{thinmathspace}{0ex}}\text{ln}\frac{{g}_{A}{h}^{*}[{v}_{\theta}-{E}_{K}]}{f({v}_{\theta},{w}_{FP})}$$

(18)

for all *g _{A}* > $\hat{{g}_{A}}$. Thus, as

When *g _{A}* > $\hat{{g}_{A}}$ but is not too large, the map Π intersects the diagonal along its upper branch at a stable fixed point (Fig. 8(a2)). For larger values of

When the point of discontinuity *h** passes through the diagonal, the transition to the period-two orbit is not immediate. In fact, for a small interval of *h** values near this point, the map Π goes through a complicated bifurcation process (Fig. 7). One way to describe this bifurcation is by starting with the period-two orbit (2:1 orbit) and increasing *h**, alternatively decreasing *g _{A}*. As

The 3:2 trajectory on the *w*-*h* phase plane. **(a)** The 3:2 periodic trajectory of F projected on the *w*-*h* phase plane shows two periodic loops indicating that the trajectory can be thought of as a combination of a 1:1 and a 2:1 solution. Different colors match **...**

First, start with the trajectory on the left branch at *t*_{0}. Between *t*_{0} and *t*_{1} = *t*_{0}+*T _{act}*,

The above analysis only partially describes the behavior of the system in the region between stable 1:1 and 2:1 orbits (Fig. 7(b)). The bifurcation structure between these two regions is in fact much more complicated. As shown, the 3:2 orbit exists for an interval of *g _{A}* values between the 1:1 and 2:1 regions. If we now zoom in to examine the region between the 1:1 and 3:2 orbits, we find an interval of

The generalized Pascal triangle shows the proposed bifurcation structure of the map Π for the parameter *g*_{A}. The top level shows the stable solutions that cover the largest intervals of *g*_{A} values. These are the intervals marked as 1–4 in **...**

Similar complicated dynamics arise in the regions between the 2:1 and 3:1 orbits and, in general, between the countably infinite number of regions containing *n:*1 and *n*+1:1 orbits (white regions of Fig. 7(a)). The dynamics of the map in these regions can be understood using an analysis similar to that described above, thus resulting in a rich bifurcation diagram that can be schematized as a countably infinite sequence of generalized Pascal triangles.

The activity time of a neuron in an oscillatory network can be affected by a combination of various intrinsic and synaptic properties [22, 9, 13, 3]. In a previous study we used a simplified model to predict the steady state value activity of biological PY neurons [29]. In the current study, we expand this model to a very general setting to explore the factors that determine the activity time of a neuron with A-current following periodic inhibitory input. The results indicate that, even with a simple model neuron, the interaction of the A-current parameters and other intrinsic parameters can be quite complex and lead to distinct and sometimes unintuitive model behaviors. Although the roles of additional factors such as synaptic dynamics and the interaction with other network neurons that can affect the role of the A-current in determining the activity remain to be explored, the basic geometric tools have been provided which could be used for such analysis in more complicated settings.

The main advance in this paper is the derivation and analysis on the (*w*-*h*) slow manifold. By making some simplifying assumptions for the dynamics away from this slow manifold, we were able to show how the behavior of F can be completely determined by its behavior on the slow manifold. In particular, we showed how the geometry of one-dimensional curves corresponding to upper and lower knees and fixed points on the manifold organized the behavior of trajectories on it. In turn, this allowed us to understand the circumstances that gave rise to different types of periodic and also potentially chaotic solutions.

The analysis of the fate of the trajectory in the *w*-*h* phase plane indicates that, in some cases, the behavior of the model neuron following inhibition is non-periodic and can be quite unintuitive and even ambiguous. These ambiguous cases are marked by the intersection between the LK and UK curves on the *w*-*h* phase plane, indicating the possibility of a trajectory reaching such an intersection point in which case its fate is unclear, as in Subcase 5 in Case I (Fig. 4(b)) or Case II.c (Fig. 5(f)). The dynamics of the trajectories on MB suggest how such ambiguous and potentially chaotic solutions of the full system may arise. Depending on the time constants τ_{wm} and τ_{hm}, any trajectory leaves MB through UK and return to the silent state or through LK and go to the active state. Denote the point of intersection in the *w*-*h* phase plane by *Q* (Figs. 4(b) and 5(f)). Let Γ denote the trajectory under Equation (9) that contains *Q*. By uniqueness Γ forms a separatrix in the *w*-*h* phase space separating solutions that leave MB through LK from those that leave through UK. Thus there is no ambiguity of the fate of the trajectory from the behavior of the (reduced) flow on MB. Yet the flow on MB is only the slow component of the full 3-dimensional flow as described on the 2-dimensional slow manifold. In cases where the trajectory does not always land on one side of Γ it is not possible to determine the behavior of the full trajectory from that of the singular trajectory on MB. Thus, if the dynamics of the fast system allow the trajectory, in different cycles, to land on different sides of Γ when it reaches MB, the long-term behavior over multiple cycles can be chaotic (Fig. 5(f)). Even if there is no chaotic behavior in such a case, the global dynamics of the actual biological system can be ambiguous from cycle to cycle due to small perturbations from noise.

Several studies of neuronal systems have used reduction to low-dimensional maps to prove the existence and stability of solutions [6, 15, 14, 4, 19, 18], most by tracking state variables in a low-dimensional phase space. The analysis on the slow manifold in our study made it possible to derive a one-dimensional map Π of the unit interval whose dynamics predicted the behavior of the full set of differential equations (1). Π has a discontinuity *h** that depends on *g _{A}*. For small values of

The map Π was derived under the assumption that the middle branch of the *v*-nullcline was vertical, that is *v=v _{θ}*, and that τ

In the case where τ_{wm} is not necessarily small, as in Cases 2 and 4 of Fig. 5(d), we would need to keep track of the values of both *h* and *w* at the end of each cycle of inhibition. The ensuing map would then be two-dimensional and we would have to track whether the trajectory crosses LK in each cycle. A simple way to make this calculation would be to assume a linear relationship between *w* and *h* along LK as has been done in other contexts for slow manifolds [2, 27]. This would allow us to calculate *t _{m}^{n}*, the time spent on the middle branch, allowing us to use equations (16)(a) and (b) to update

In our model we have smoothed over the spikes of the biological follower neuron (see, e.g., the biological PY neuron in Fig. 1) and only modeled the burst envelope. Including spikes in our analysis by choosing a bursting model for the follower cell F would certainly bring the model closer to the biological context. However, because of our assumption of a fixed point in the active state of F, the addition of spikes would not change the time when F returns to the silent state. This time is still controlled by the inhibition from O. As a result, the *h* value giving the level of the *A*-current would also not change if spikes were included. If we remove the fixed point in the active state of F, then there is the possibility that the nature of the local bifurcation and geometry of the nullclines in the bursting model could affect the return time to the silent state. Even then, the *h*-dynamics would still remain independent of the details of the individual spikes and would depend only on the time that F is active.

The analysis presented here lays the groundwork for several interesting extensions. For example, the feed-forward O-F network described can be considered one part of a two-cell feedback loop that is commonly found in CPG networks. An open question worth exploring is how the feedback inhibition from F to O affects the ensuing dynamics. This would depend in part on whether O is a Type I or II oscillator as the A-current would affect the timing of the feedback to O. A second and related question is how the role of the A-current may change if the inhibition between O and F is subject to short-term synaptic plasticity. In an earlier work, we have already addressed the feed-forward situation [2], so it would be interesting to study the feedback situation. A third possible extension would be to consider a pair of reciprocally coupled Morris-Lecar type cells each with an A-current and include synaptic plasticity to understand the dynamics of the feedback network. In all of these cases, reduction techniques similar to the ones described here that allow the dynamics to be projected onto a lower dimensional slow manifold and on which a low-dimensional map can be derived will be critical in any attempt to understand the dynamics of the full system.

This work was supported by NSF DMS0615168 (AB) and NIH MH-60605 (FN).

- O
- The oscillator neuron
- F
- The follower neuron; inhibited by O
- LB
- Left branch of the cubic
*v*-nullcline of F - LB
^{1} - LB while inhibited
- MB
- Middle branch of the cubic
*v*-nullcline of F - MB
^{1} - MB while inhibited
- RB
- Right branch of the cubic
*v*-nullcline of F - RB
^{1} - RB while inhibited
- LK
- Lower knee of MB
- LK
^{1} - Lower knee of MB
^{1} - UK
- Upper knee of MB
- UK
^{1} - Upper knee of MB
^{1} - FP
- Fixed point on MB
- FP
_{1} - Fixed point on LB
^{1} - FP
_{2} - Fixed point on RB
- Q
- Intersection of LK and UK in the
*w*-*h*phase plane - T
_{act} - Duration of the active phase of O during which it inhibits F
- T
_{in} - Duration of the inactive phase of O
- τ
_{wl} - Time constant of
*w*on LB or LB^{1} - τ
_{wm} - Time constant of
*w*on MB or MB^{1} - τ
_{wh} - Time constant of
*w*on RB or RB^{1} - τ
_{hl} - Time constant of
*h*on LB or LB^{1} - τ
_{hm} - Time constant of
*h*on MB or MB^{1} - τ
_{hh} - Time constant of
*h*on RB or RB^{1}

The XPP code:

### define basic parameters ###

C=1

p dur=500

p period=1000

### define the currents ###

# external current

I_ext=75

# the leak current

gl=2

el=−60

I_L(x)=gl*(x-el)

# the Ca current

eca=120

gca=4

kca=18

vca=−1.2

m_cainf(x)=0.5*(1+tanh((x-vca)/kca))

I_Ca(x)=gca*m_cainf(x)*(x-eca)

# the K current

gk=8

ek=−84

vk=15

kk=5

w_inf(x)=1/(1+exp(-(x-vk)/kk))

tk1=10

tk2=300

w_tau(x)=tk1+tk2*heav(x-10)

I_K(x,y)=gk*y*(x-ek)

# the A current

# gA =4, 8, 20 and 5.

gA =4

p vm=−6

p km=0.5

m_inf(x)=1/(1+exp(-(x-vm)/km))

h_inf(x)=1 − heav(x-vm+5)

th1=495

th2=485

th3=800

th4=500

h_tau(x)=th1-th2*heav(x+30)+th3*(heav(x+20)-heav(x-0))+th4*heav(x-10)

I_A(x,y)=gA*m_inf(x)*y*(x-ek)

# the synaptic current

p g_syn=1.2

vO(y)=−50+50*heav(dur-mod(y,period))

I_syn(x,y)=g_syn*(x+80)/(1+exp(-(vO(y)+10)/0.1))

### define the ODEs ###

v'=(I_ext-I_Ca(v)-I_K(v,w)-I_L(v)-I_A(v,h)-I_syn(v,t))/C

w'=(w_inf(v)-w)/w_tau(v)

h'=(h_inf(v)-h)/h_tau(v)

### set up the initial values ###

init v=−41.885

init w=0

init h=0.5

aux hga=h*gA

### set up the conditions for the ODE solver ###

@ total=12000,dt=0.1,xlo=−100,xhi=50,ylo=−0.2,yhi=1.1

@

nmesh=200,xplot=v,yplot=w,maxstor=1000000,

method=stiff

done

Yu Zhang, Department of Mathematical Sciences, New Jersey Institute of Technology Newark, NJ 07102 ; Email: ude.tijn@gnahz.uy.

Amitabha Bose, Department of Mathematical Sciences, New Jersey Institute of Technology Newark, NJ 07102 ; Email: ude.tijn@esob.

Farzan Nadim, Department of Mathematical Sciences, New Jersey Institute of Technology, Department of Biological Sciences, Rutgers University, Newark, NJ 07102 ; Email: ude.tijn@nazraf.

1. Agmon A, Wells JE. The role of the hyperpolarization-activated cationic current I(h) in the timing of interictal bursts in the neonatal hippocampus. J Neurosci. 2003;23:3658–3668. [PubMed]

2. Bose A, Manor Y, Nadim F. The activity phase of postsynaptic neurons in a simplified rhythmic network. J Comput Neurosci. 2004;17:245–261. [PubMed]

3. Brunel N, Wang XJ. What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J Neurophysiol. 2003;90:415–430. [PubMed]

4. Butera RJ. Multirhythmic bursting. Chaos. 1998;8:274–284. [PubMed]

5. Connor JA, Stevens CF. Prediction of repetitive firing behaviour from voltage clamp data on an isolated neurone soma. J Physiol. 1971;213:31–53. [PubMed]

6. Ermentrout GB, Kopell N. Fine structure of neural spiking and synchronization in the presence of conduction delays. Proc Natl Acad Sci U S A. 1998;95:1259–1264. [PubMed]

7. Fenichel N. Geometric singular perturbation theory for ordinary differential equations. J Diff Eqn. 1979;213:53–98.

8. Gerber B, Jakobsson E. Functional significance of the A-current. Biol Cybern. 1993;70:109–114. [PubMed]

9. Grillner S. Biological pattern generation: the cellular and computational logic of networks in motion. Neuron. 2006;52:751–766. [PubMed]

10. Harris-Warrick RM, Coniglio LM, Barazangi N, Guckenheimer J, Gueron S. Dopamine modulation of transient potassium current evokes phase shifts in a central pattern generator network. J Neurosci. 1995;15:342–358. [PubMed]

11. Herrington J, Lingle CJ. Multiple components of voltage-dependent potassium current in normal rat anterior pituitary cells. J Neurophysiol. 1994;72:719–729. [PubMed]

12. Huguenard JR, Coulter DA, Prince DA. A fast transient potassium current in thalamic relay neurons: kinetics of activation and inactivation. J Neurophysiol. 1991;66:1304–1315. [PubMed]

13. Huguenard JR, McCormick DA. Thalamic synchrony and dynamic regulation of global forebrain oscillations. Trends Neurosci. 2007;30:350–356. [PubMed]

14. Lee E, Terman D. Uniqueness and stability of periodic bursting solutions. J Diff Eqn. 1999;158:48–78.

15. LoFaro T, Kopell N. Timing regulation in a network reduced from voltage-gated equations to a one-dimensional map. J Math Biol. 1999;38:479–533. [PubMed]

16. Manor Y, Nadim F. Synaptic depression mediates bistability in neuronal networks with recurrent inhibitory connectivity. J Neurosci. 2001;21:9460–9470. [PubMed]

17. Marder E, Bucher D. Understanding circuit dynamics using the stomatogastric nervous system of lobsters and crabs. Annu Rev Physiol. 2007;69:291–316. [PubMed]

18. Matveev V, Bose A, Nadim F. Capturing the bursting dynamics of a two-cell inhibitory network using a one-dimensional map. J Comput Neurosci. 2007;23:169–187. [PMC free article] [PubMed]

19. Medvedev G. Reduction of a model of an excitable cell to a one-dimensional map. Physica D. 2005;202:37–59.

20. Mishchenko EF, Rosov NK. Differential equations with a small parameter and relaxation oscillations. New York: Plenum Press; 1980.

21. Morris C, Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophys J. 1981;35:193–213. [PubMed]

22. Narayanan R, Johnston D. The h channel mediates location dependence and plasticity of intrinsic phase response in rat hippocampal neurons. J Neurosci. 2008;28:5846–5860. [PMC free article] [PubMed]

23. O'Keefe J, Recce ML. Phase relationship between hippocampal place units and the EEG theta rhythm. Hippocampus. 1993;3:317–330. [PubMed]

24. Poe GR, Nitz DA, McNaughton BL, Barnes CA. Experience-dependent phase-reversal of hippocampal neuron firing during REM sleep. Brain Res. 2000;855:176–180. [PubMed]

25. Rinzel J, Ermentrout GB. Analysis of neural excitability and oscillations. In: Koch C, Segev I, editors. Methods in Neuronal Modeling: From Ions to Networks. Cambridge, MA: MIT Press; 1998. pp. 251–291.

26. Selverston AI. A neural infrastructure for rhythmic motor patterns. Cell Mol Neurobiol. 2005;25:223–244. [PubMed]

27. Terman D, Lee E. Partial Synchronization in a Network of Neural Oscillators. SIAM J Appl Math. 1997;57:252–293.

28. Wustenberg DG, Boytcheva M, Grunewald B, Byrne JH, Menzel R, Baxter DA. Current- and voltage-clamp recordings and computer simulations of Kenyon cells in the honeybee. J Neurophysiol. 2004;92:2589–2603. [PubMed]

29. Zhang Y, Bose A, Nadim F. Predicting the activity phase of a follower neuron with A-current in an inhibitory network. Biol Cybern. 2008;99:171–184. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |