|Home | About | Journals | Submit | Contact Us | Français|
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 gA, the maximal conductance of the A-current, or other model parameters. We demonstrate that, as the value of gA is varied, the trajectory of the follower neuron goes through a set of bifurcations to produce n:m periodic solutions where the follower neuron becomes active m times for each n cycles of the oscillator. Using a generalized Pascal triangle, each n:m trajectory can be constructed as a combination of solutions from a higher level of the triangle.
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 . 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 .
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 gA of the A-current. The bifurcation diagram of the map using gA as a parameter shows the existence of a surprisingly rich set of periodic solutions. These solutions correspond to different types of n:m periodic solutions of the full model in which the follower neuron becomes active m times for every n cycles of the oscillator. Through our analysis, we conjecture that there is a countable number of compact intervals of gA, in each of which, there exists infinitely many periodic solutions with distinct n:m values that are arranged in a generalized Pascal triangle.
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 vO of O as a square wave oscillating periodically between −50 mV and 0 mV with active duration Tact and inactive duration Tin (Fig. 2(a)). For the follower neuron F, two of the model variables (v and w) are from the Morris-Lecar (ML) model representing the membrane potential and the activation fraction of the potassium current, and one variable h describes the A-current inactivation [2, 29]. The spikes are smoothed over since they do not play an important role in determining the effect of the A-current. The model can be described as
represents the ML terms and
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:
Here kx is negative for activation variables and positive for inactivation variables. The time constants τw (v) and τh(v) determine the speed with which the variables w and h change in different voltage regions. In previous studies [2, 29], the parameters kw, kn and kh were fixed, and kn was set to be small in order to guarantee a steep A-current activation curve n∞(v). In the current study we relax these conditions by considering different values of these parameters.
The synaptic current from O to F is given by
where s∞(vO) 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 LB1, MB1 and RB1, 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 .
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 : 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):
The fast equations are obtained by rescaling t = εξ in Equation (1) and then setting ε = 0:
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 LB1 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 FP1 and FP2 (Fig. 2(c)). Similarly, depending on the shape and position of the w-nullcline, there may be a stable fixed point on MB or MB1 (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 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 . 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 FP1 on LB1 during Tact; it stays near FP1 due to its stability. When the inhibition ends, the v-nullcline is raised as the negative term Iinh disappears in Eq. (6); FP1 no longer exists and therefore the trajectory obeys the fast equations (Eq. (7)) and jumps to MB. The trajectory is also able to jump when it reaches the knee (points where F/v = 0) on any branch. For example, when the trajectory encounters a lower knee of MB, it jumps to RB as determined by the saddle-node bifurcation in the fast equations (Eq. (7)). For the middle branch MB, the upper and lower knees are dependent on h, and are referred to as UK and LK.
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) = wFP1 along LB1 and when the trajectory lands on MB, it has a w value equal to wFP1. We also assume that h∞(v) = 1 on LB and LB1 and h∞(v) = 0 on MB, RB, MB1 and RB1. The time constants of w and h on LB, MB and RB (or their inhibited counterparts) are, respectively, referred to as τwl and τhl, τwm and τhm, and τwh and τhh. The reduced equations for the trajectory moving on each of the six branches are as follows:
When a trajectory leaves LB1, the trajectory may land on MB or RB. It will land on MB if gA is sufficiently large and τh1 and τw1 is sufficiently small. The first and second conditions ensure that the minimum w-value along MB is smaller than wFP1, while the third condition guarantees that F is in a neighborhood of wFP1 at the time the trajectory leaves LB1.
All time constants are assumed to be of the same order of magnitude as Tact (O(Tact)) unless otherwise specified in the Results. In cases where two time constants are assumed to be of different orders of magnitude (i.e. τ1 << τ2), the larger one (τ2 in this case) is assumed to be O(Tact).
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 . 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 kn→0. The consequence of this assumption in the v-w phase plane is that the middle branch (MB) of the v-nullcline becomes vertical. Additionally, we had only considered the effects of strong synaptic inhibition . In the current study, we relax these two simplifying assumptions, which allows for a wider range of interactions between the A-current and the other intrinsic and synaptic parameters, resulting in a variety of behaviors in the follower neuron, some of which are quite unintuitive. By following the trajectory of the follower neuron in two distinct phase planes, each a projection of the full three-dimensional phase space, we analyze these potential behaviors and determine which combinations of synaptic and intrinsic parameters result in each behavior.
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 LB1. During the duration of inhibition Tact, the trajectory moves downwards along LB1 towards the stable fixed point FP1 (Figs. 3(a) and 3(b)). Simultaneously, the inactivation variable h of the A-current increases with time constant τhl. The growth of h does not affect the shape of LB, but does result in a larger MB. Thus, the size of MB is related to how long the follower neuron stays in the inactive state (time spent on LB1). We assume that τwl is small enough on LB1 such that, when inhibited, the trajectory reaches the neighborhood of FP1. After the inhibition has finished, the trajectory in the v-w plane is released from a neighborhood of the fixed point FP1 and jumps horizontally to MB. Once the trajectory lands on MB, h begins to decay with time constant τhm which causes MB to shrink from the lower knee (LK; Figs. 3(a) and 3(b)). At the same time, the trajectory moves up, in the direction of increasing w, with time constant τwm.
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 Tin. In Case I, however, the w-nullcline is above the middle branch, which allows the trajectory to jump left if it moves fast enough to reach UK during Tin. Although these two possibilities are not exhaustive, they indicate how the structure of the nullclines can determine the fate of the trajectory on MB.
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 (wFP1, hstart). 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 w-h phase plane is that the fate of the trajectory can be determined by comparing the kinetics of the two slow variables, w and h, as will be described later in the Results.
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 (Tin) and active (Tact) duration of the pacemaker and the geometry of the nullclines. The analysis can largely be broken up into two distinct cases. Case I occurs when there are no fixed points on MB, while Case II occurs when the curve of fixed points FP exists on MB.
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 , 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 (wFP1, hstart) 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 (τ*wm, τ*hm).
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.
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 tm that the follower neuron spends on the middle branch and, in our previous study , 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 n:m periodic solutions.
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 gA or, alternatively, τhm is sufficiently large. To follow the trajectory in this case, consider the gKw-gAh phase plane where the maximal conductances are explicitly incorporated along each axis (Fig. 6(b)). For large values of gA, the starting position at t = Tact as given by will be large. Thus, it will take longer for the trajectory to decay to LK. Moreover, if the trajectory cannot reach LK in one cycle, the inhibition from the pacemaker moves LK down in the gKw - gAh phase plane for the time duration equal to Tact (labeled LK1 in Fig. 6(b)). During this time, the follower neuron may or may not reach LK1 (it doesn’t in the example of Fig. 6) but it does cross over LK. Once the inhibition is removed, the trajectory finds itself below LK and no longer on the middle branch and thus jumps to RB. Thus, for large gA, the pacemaker may go through multiple cycles (two in Fig. 6) while the follower neuron remains on MB. Such trajectories correspond to the n:1 solutions.
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 hn denote the value of h in cycle n at the moment t = (n−1)P + Tact when the inhibition from the pacemaker is removed. Let tmn denote the time spent on MB in the nth cycle, vθ the v value of the middle branch, and wFPthe w value of FP on MB. If tmn < Tin, the value of hn is calculated by integrating the third equations in Eqs. (11), (9) and (10) on the LB1, MB and RB branches. If, on the other hand, tmn > Tin, then the trajectory stays on the middle branch for at least one more cycle. In this case, the value of hn at the next cycle is given by
The values of tmn are calculated by setting the right-hand side of the first equation in Eq. (1) equal to zero:
using the assumption that the middle branch is vertical with v=vθ and the assumption that τwm is small (w reaches a neighborhood of wFP). The more general case is addressed in the Discussion. The map Π describing the value hn is defined as
If tmn < Tin for all n, then Π follows Equation (16)a in each iteration, and it converges to a steady state h* (the analytic form is given in Zhang et al, ). However, when tmn > Tinfor some n, the trajectory is not able to jump to the right branch during Tin, and therefore Π does not converge to h*.
Using XPP we computed a bifurcation diagram for Π as a function of gA (Fig 7). For small values of gA, Π has a single stable fixed point h* (region 1 in Fig. 7(a)). As gA is increased, Π bifurcates in a very complicated manner as highlighted in the boxed regions of Figs. 7(b) and 7(c). For larger values of gA, Π has stable period 2, 3, 4 … orbits as marked in Fig. 7(a). These correspond to cases where the trajectory spends 2, 3, 4 or more cycles on MB before jumping to the active state and are referred to as n:1 periodic orbits for n = 1, 2… Figure 7(b) is a blow up of Fig. 7(a) corresponding to the portion lying between regions 1 and 2. Figure 7(c) shows a blow up of the boxed region in Fig. 7(b) near h = 0.8. As seen in these figures, as gA is decreased, Π bifurcates to higher order stable periodic orbits, apparently through a period-adding process. To understand this and the other solutions found above, we solved the full set of Equations (1) for the values of gA marked as a, b, c and d in Fig. 7(a). These traces along with the cobweb diagram of the associated map Π are shown in Fig. 8. The first thing to note is that for the shown values of gA, the map is discontinuous. This is because of the definition of Π (Equation (16)) that allows Π to take on different values depending on the value of tmn relative to Tin and will be discussed further below. In Fig. 8(a2), one of the branches of the map intersects the diagonal (with slope less than one) and thus there exists a stable period one solution. This solution corresponds to a 1:1 periodic trajectory of O with respect to F (Fig. 8(a1)). In Fig. 8(b2), a period-two orbit is shown. Here the trajectory spends two cycles on MB before going to the active state resulting in a 2:1 trajectory (Fig. 8(b1)). Figs 8(c2) and 8(d2) show two distinct period-three orbits. In Fig. 8(c2), gA=20 nS and the trajectory spends two and a half cycles away from LB (i.e. two cycles on MB and half a cycle on RB) resulting in a 3:1 orbit. In contrast, in Fig. 8(d2), gA=5 nS, and the trajectory alternates between spending one-half cycle away from LB followed by one and a half cycles away from LB resulting in a 3:2 orbit. This is a typical example of the complicated dynamics arising in the regions of the bifurcation diagram of Figs. 7(b) and 7(c). These complex dynamics depend on the fact that the inactivation time constant of the A-current is slow; i.e., τhh is of the same order as τhm and both are sufficiently large. This allows h to be bounded away from zero when F returns to the silent state, thus allowing the trajectory it to carry information about the history of h forward from prior cycles.
To explain how the map Π captures the complicated dynamics of the trajectory, we consider how the discontinuity in Π depends on gA. Let h* denote the point of discontinuity in Π. We show that the value of h* is a decreasing function of gA. For small values of gA, the map has no discontinuity. This is because no trajectory spends a time larger than Tin on MB for any initial condition (even with h=1) and the map Π is always defined by Equation (16)a. The discontinuity in Π occurs for the values of gA > satisfies
which is Equation (16)c with tmn=Tin and hn−1=1. In this case, for hn−1 values close to 1, the time spent on MB could become larger than Tin and thus provide a choice in the value of hn between Equations (16)a and (16)b. Note that hn−1 values close to 1 would result in Equation (16)b for hn. The discontinuity point h* satisfies the equation
for all gA > . Thus, as gA increases, h* must decrease.
When gA > 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 gA (smaller h*), Π has no fixed points (Fig. 8(b2)) but Π2 acquires two stable fixed points corresponding to the stable period-two orbit of Π. If gA is yet larger, a period-three orbit appears (Fig. 8(c2)). Note that the length of the lower branch of Π increases which allows for two iterations of the cobweb trajectory to fall on this branch. As gA is increased further, higher order orbits arise of the n:1 type, where n−1 iterations of the stable periodic orbit will fall on this branch. These iterations correspond to the number of the cycles of inhibition for which the trajectory of F remains on MB.
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 gA. As h* is increased, the upper branch of Π gets longer. This allows the cobweb trajectory to spend two iterates on the upper branch before returning to the lower branch (compare Figs. 8(b2) and 8(d2)) which, in turn, results in a period-three orbit. This period-three orbit corresponds to a 3:2 solution where the oscillator fires three times for every two firings of the follower neuron. As an example of the complicated dynamics of the map Π in the area between regions 1 and 2, we will describe this period-three orbit (corresponding to a 3:2 solution; Fig. 8(d)) in the w-h phase plane. This trajectory can be thought of as a combination of a 1:1 and a 2:1 trajectory. The projection of one full period-three orbit onto the w-h phase plane is shown in the Fig. 9(a) with the corresponding voltage trace in the inset. This trajectory arrives at MB at two different times (t1 and t5; inset of Fig. 9), with two distinct values of h (respectively, hlo and hhi; filled and open circles in Fig. 9(b); expansion of shaded region of Fig. 9(a)) and the same value of w=wFP1. We will track both paths of this trajectory on MB.
First, start with the trajectory on the left branch at t0. Between t0 and t1 = t0+Tact, h increases while w decreases (blue curve). This trajectory jumps to MB at t1 to the point (wFP1, hhi) (open circle; Fig. 9(b)). The trajectory then moves on MB (green curve) and lies above LK at the onset of the next inhibition (t2 = t1+Tin; grey dashed arrow) but below LK at the end of this inhibition at t3 = t0+P (end of the green curve). Therefore, in this case the time spent on MB is tm=P. Once the trajectory jumps to the active state, h continues to decay on RB between t3 and t4 (end of bottom black trajectory in Fig. 9(a)). At t4 = t0+2P, the trajectory lies on LB1 and h increases (yellow curve). At t5 = t4+Tin, the trajectory lands back on MB at (wFP1, hlo) and continues on MB (red curve) until it reaches LK at t6. Note that in this case tm = t6 − t5 < Tin because the trajectory jumps to RB before the onset of the next inhibition. Once the trajectory jumps back to LB, it will repeat the path from the blue portion. In conclusion, each time the trajectory lands on MB, the values of h alternate between hhi and hlo.
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 gA 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 gA values for which there is a 4:3 orbit. Similarly, we find 5:3 orbits between the 3:2 and 2:1 regions. Further zoom-ins between these two regions reveal more periodic solutions. In fact we conjecture that, in the interval of gA values between the 1:1 and 2:1 region, there is a countably infinite set of sub-intervals of gA, such that in each sub-interval a distinct periodic orbit exists. Further, these periodic orbits can be mapped in a one-to-one manner to elements of a generalized Pascal triangle. This generalized Pascal triangle can be defined using the following addition rule: a:b c:d = (a + c):(b + d). The top level of the triangle consists of 1:1 and 2:1 solutions. The second layer has 1:1, 1:1 2:1 = 3:2 and 2:1 solutions where the 3:2 solutions exist in a gA sub-interval between those of the 1:1 and the 2:1 solutions. Any level of the triangle adds in sub-intervals that are found between adjacent subintervals of the previous higher level. This triangle is schematically shown in Fig. 10. Although we do not prove this conjecture, we numerically verified the existence of many such solutions, two of which are shown for iterates of the map Π in Fig. 10.
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 . 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 gA. For small values of gA, the 1:1 periodic solution is stable as indicated by the slope of the map Π where it intersects the diagonal (Fig. 8(a2)). As the discontinuity h* approaches the diagonal by increasing gA, this 1:1 solution retains its stability; yet, interestingly, as h* passes through the diagonal to the left (i.e. by decreasing gA), an explosion of periodic solutions occurs. This transition is a global bifurcation in the dynamics of Π since the 1:1 solution does not lose its stability; instead, it ceases to exist. In the gA sub-intervals between the 1:1 and 2:2 solutions, we conjecture there exists a family of n:m periodic solutions that can be enumerated using a generalized Pascal triangle. It would be of interest to see if these types of solutions have any significance for the biological system. A primary step would be to determine the length of the individual sub-intervals corresponding to different n:m solutions.
The map Π was derived under the assumption that the middle branch of the v-nullcline was vertical, that is v=vθ, and that τwm is sufficiently small. An obvious question is whether the assumption of a vertical middle branch is too restrictive to describe the dynamics of the full system. Figure 8 shows that this is not the case: in the simulations shown in the left panels, the middle branch is steep but not vertical. Nonetheless, the question remains on how we could extend the analysis to cases where the assumption of a vertical middle branch is relaxed. In this case, the curve FP would no longer be vertical in the w-h phase plane but would have negative slope as in Fig. 5. If the assumption that τwm is small is retained, the trajectory will rapidly approach FP and will be funneled toward the intersection of LK and FP which occurs with values v=vθ and w=wFP. Thus Equation (15) still applies and the analysis remains unchanged. Indeed this is the reason why the simulations in Fig. 8 perfectly match the output of the map.
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 tmn, the time spent on the middle branch, allowing us to use equations (16)(a) and (b) to update hn. We can similarly derive a recursive equation for wn and thereby be able to analyze the ensuing two-dimensional map.
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 , 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).
The XPP code:
### define basic parameters ###
### define the currents ###
# external current
# the leak current
# the Ca current
# the K current
# the A current
# gA =4, 8, 20 and 5.
h_inf(x)=1 − heav(x-vm+5)
# the synaptic current
### define the ODEs ###
### set up the initial values ###
### set up the conditions for the ODE solver ###
Yu Zhang, Department of Mathematical Sciences, New Jersey Institute of Technology Newark, NJ 07102 ; Email: email@example.com.
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.