Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Math Biosci. Author manuscript; available in PMC 2017 August 1.
Published in final edited form as:
PMCID: PMC4940261

The role of electrical coupling in generating and modulating oscillations in a neuronal network


A simplified model of the crustacean gastric mill network is considered. Rhythmic activity in this network has largely been attributed to half center oscillations driven by mutual inhibition. We use mathematical modeling and dynamical systems theory to show that rhythmic oscillations in this network may also depend on, or even arise from, a voltage-dependent electrical coupling between one of the cells in the half-center network and a projection neuron that lies outside of the network. This finding uncovers a potentially new mechanism for the generation of oscillations in neuronal networks.

Keywords: electrical coupling, half-center oscillator, phase-plane analysis

1 Introduction

Networks of neurons display a variety of oscillatory behaviors. For example, oscillations in the levels of calcium concentrations, gene expressions and in the membrane voltage across cell membranes are all commonly found in neuronal systems. Often these oscillations are rhythmic in that they display a consistent pattern at a prescribed frequency [1]. Central pattern generating (CPG) networks provide several examples that exhibit rhythmic activity. CPGs refer to networks of neurons in the central nervous system that produce patterned (usually oscillatory) activity in the absence of patterned sensory input. These networks play a critical role in generating a diverse array of motor functions such as digestion, locomotion, respiration and regulation of heartbeat in invertebrates [2]. A central question in the study of neural oscillations is what are the mechanisms that underlie the generation of rhythmic activity and how that activity is regulated. This study will focus on this general question in the context of the gastric mill rhythm (GMR; frequency 0.1 Hz) that arises in the stomatogastric ganglion (STG) in the crustacean central nervous system. In particular, we will show the existence of a new mechanism based on voltage-dependent electrical coupling for generation of oscillations within a neuronal network.

The gastric mill network consists of a small number of neurons in the STG that control muscles that move teeth to provide grinding of food (chewing) within the gastric mill stomach of crustaceans [3]. In the Jonah crab, a pair of neurons, the lateral gastric (LG) and Interneuron 1 (INT1) form a half-center oscillator (HCO) and are primary contributors to the GMR. These neurons are connected by reciprocally inhibitory synapses and, during gastric mill activity, display anti-phase bursting oscillations. They also receive input from various parts of the stomatogastric nervous system (STNS). In particular, INT1 receives rhythmic inhibition from the pacemaker anterior burster neuron (AB) of the pyloric CPG. Because the pyloric rhythm (frequency 1 Hz) is much faster than the gastric mill, the AB to INT1 input produces pyloric timed patterns in the INT1 bursting activity. Both LG and INT1 receive excitatory input from the modulatory commissural neuron 1 (MCN1) with INT1 receiving fast excitation and LG receiving slow modulatory excitation. Additionally, the MCN1 axon terminals are electrically coupled to LG in a manner that is dependent on the voltage of LG [5]. It is the role of this electrical coupling that is of particular interest to us in this paper.

Neurons that lie within a HCO typically utilize reciprocal inhibition to generate oscillations [6]. In particular, in a two cell HCO, when one of the cells is active, its inhibitory synapse suppresses the other. At some later time, the silent cell escapes or is released from inhibition and the roles of the two cells switch [7]. In the gastric mill network, LG and INT1 can oscillate in this manner with the ability to escape inhibition and generate oscillations, but only in the presence of the excitatory input provided by MCN1 [5, 8].

Although a number of modeling studies have explored the generation of oscillations in the gastric mill network [8, 9, 10, 11, 13], the role of the strong electrical coupling between the MCN1 axon terminals and the LG neuron has not been previously explored. In this study, we will show that voltage-dependent electrical coupling can provide an alternative mechanism for the generation of oscillations when the inhibition based HCO mechanism is incapable of doing so. In particular the LG – INT1 HCO can be rendered ineffective if 1) the inhibitory synapse form INT1 to LG is inactivated, or 2) if the excitability property of LG is reduced. In order to fully understand how electrical coupling affects this network, we will first consider a simple model to see how how electrical coupling between LG and MCN1 axon terminals affects the ability of oscillations to be created through the standard HCO inhibition based mechanism. We will discuss how the electrical coupling modulates the rhythmic properties of this oscillation. We will then remove the INT1 to LG synapse and show that rhythmic oscillations can still arise through the electrical coupling between LG and MCN1 axon terminals, but only if this coupling is voltage dependent, as has been reported experimentally [5]. We will then demonstrate the same in a biophysical model based on the Morris-Lecar equations [15]. For both models, we derive conditions on parameters showing why the electrical coupling must be voltage dependent to produce oscillations.

The modeling and analysis in this paper is based on the use of geometric singular perturbation theory. Exploiting inherent differences in timescales, we will derive sets of fast and slow equations that can be studied in the relevant phase space. For the simple model, this can be done on a two-dimensional phase plane and is the focus of Sections 3.1–3.4. The analysis in those sections follows the tradition of using relaxation oscillators with the individual neurons modeled as passive elements. The relaxation oscillations in this case arise due to the method of model reduction that incorporates a slow synaptic variable. In Section 3.5, the fast-slow analysis allows us to project the relevant dynamics onto two different phase planes to facilitate understanding of the model.

2 Model

2.1 Simple passive cell network model

We describe the simple network that we shall initially consider. A key assumption for this model is that INT1 and LG are modeled as passive cells with no active currents or excitable properties. Thus if oscillations are to be generated, they must arise as a direct result of network interactions. By identifying variables that evolve on different time scales and by making a few other assumptions, we can use geometric singular perturbation theory to focus on the analysis of a reduced two-dimensional system of equations. These variables correspond to the voltage of LG and to the synaptic input that LG receives from MCN1 and are shown in solid in Fig. 1. The electrical coupling is also shown in solid in Fig. 1 as it can be defined in terms of the reduced quantities including the voltage of LG. Shown with dotted lines/circles are the other variables that we will incorporate into the solid variables and thus will not need to explicitly track.

Figure 1
Schematic diagram of the modeled network

Let VL and VI denote the voltages of LG and INT1 respectively. We will not model individual spikes but instead keep track of when a cell is above (active) or below (silent) threshold. These voltages will evolve on a fast time scale. Notice that AB and MCN1 do not receive synaptic input from any other cells in the circuit. Thus we do not explicitly model either but instead need only keep track of their synaptic and electrical output. The equations that describe the relevant voltages are:



The intrinsic current Irest,x(Vx) = grest,x[VErest,x] where grest,x and Erest,x are the passive rest conductance and reversal potentials. Notice that in the absence of any other currents, the value V = Erest,x is a stable rest point. For LG, Erest,L < VT while for INT1, Erest,I > VT for a fixed threshold VT. MCN1 is assumed to be tonically active which we model by setting its voltage to a value VM > VT. The synaptic currents obey an equation of the form Isyn,xy = gxysxy[VyEinh] where x and y are the pre-and post-synaptic cells. The variables sABI, sLI and sIL are straight forward to understand and are instantaneous. The synaptic variable sABI provides the input due to AB activity and is modeled using a periodic, half-sine function with an amplitude of 1 and period of 1 second. This synapse takes on the value one when the sine function is greater than a threshold, set here to 0.5, and is zero otherwise. The synapses between LG and INT1 are also instantaneous and we utilize the fact that these cells are always out-of-phase with one another.


The remaining synaptic variable s requires some explanation. In the biological system, MCN1 exerts a slow excitatory effect on LG that is modulated by pre-synaptic inhibition from LG onto the MCN1 to LG synapse. Thus when LG is active, this excitation is slowly removed; when LG is silent, the excitation slowly builds. This is modeled by the variable s that evolves on a slow time scale and is the only slow variable in our model. Equations governing this variable are:


In equation (1), the synaptic current is then given by


Figure 1 shows an electrical coupling between LG and the MCN1 axon terminals. The electrical current is given by


This coupling is dependent on the voltage of LG and MCN1 in two different ways. First, the strength is an increasing function of VL. The dependency of the conductance gelec on VL is incorporated using an increasing sigmoidal function n(VL). Second this strength is dependent on the driving force which is the difference between the LG and MCN1 voltages. In the biological system, the electrical coupling has a minimal effect on the MCN1 voltage[5], almost as if the electrical coupling were rectifying. We model this by simply keeping the MCN1 voltage fixed at VM independent of the value of VL. We define




where vel is the half activation value at which n(VL) = (1 − gmin)/2 and kel is the reciprocal of the slope at that point. The asymptotic value of n(VL) as VL → −∞ is denoted by gmin [set membership] (0, 1) and is the smallest positive value of the electrical conductance.

While, equations (1)(8) govern the flow of the gastric mill circuit, the dynamics can be simplified by exploiting the small parameter ε that demarcates the fast and slow time scales, as was first done by Kintos et. al. [11]. Set ε = 0 in (1)–(2). The latter of these equations can be rewritten in terms of VL and of the independently controlled quantity sABI. Namely, from (2), note that we can solve for VI = h1(VL, sABI); see the Appendix. Thus the set of equations governing the slow flow can be reduced to


Denote the right-hand side of (11) by F (VL, s). The first equation constrains the flow to lie on F (VL, s) = 0, and slaves the evolution of VL to s which is governed by the second equation (12). Rescale t = ετ, then set ε = 0 to obtain the fast equations


Equations (1314) govern the fast jumps that a trajectory in the phase plane makes between different possible (stable) branches of the VL-nullcline. For ε small enough, an actual solution to (1)–(8) lies O(ε) close to a singular periodic orbit which is pieced together from solutions of (11)–(14).

The VL nullcline is the set of points {(VL, s) : F (VL, s) = 0} and can be graphed by explicitly solving for s to obtain


The s-nullcline is simply the Heaviside function given by s = 1 when VL < VT and s = 0 when VL > VT. We could smooth this nullcline out to a sigmoid with no qualitative change in results.

The shape of the VL-nullcline is dependent on our choice of parameters. It is known from prior modeling work of this system [11, 12], and of many others in different contexts, that when one of the nullclines is cubic shaped and the other is linear or sigmoidal that oscillations may occur if the nullclines intersect on the middle branch of the cubic. In the results section below we will show how various parameters related to both the synaptic and electrical coupling affect the shape of the VL nullcline and allow it to be a cubic.

2.2 Biophysical model

In Section 3.5, we will use the Morris-Lecar equations to model both LG and INT1. As a result of the added dimensionality of the model, we will not be able to reduce the analysis to a two-dimensional phase plane. However, similar to our analysis with the simple model, we will be able to show that the projection of the LG trajectory onto two distinct two-dimensional phase planes will be crucial to understanding the role of voltage-dependent electrical coupling. When parameters are chosen in the Morris-Lecar equations to reduce the excitability of LG, the inhibition based HCO becomes ineffective. In that case, as in the case of the simple model, electrical coupling will be able to produce oscillations but only when it is voltage-dependent. Details of the model will be provided in Section 3.5 and the Appendix.

3 Results

3.1 Oscillations that arise through the INT1-LG reciprocal inhibition

For completeness and for ease in explaining the role of the voltage dependent electrical coupling, we begin by reviewing the case when gelec = 0 as described in [10]. Oscillations in this case arise as a direct consequence of the mutually inhibitory pair INT1 and LG. Because of different synaptic strengths between the two and different time constants in the active and silent states of LG, the cells form an asymmetric half-center oscillator (HCO) in that the duty cycle of each cell is not equal to 1/2. They do, however, oscillatate in anti-phase where only one of the cells is active at any moment in time.

First set gABI = 0 meaning that AB inhibition to INT1 is absent. We choose similar parameter values to [10] such that the VL-nullcline is then a cubic shaped curve where the left and right branches are positively sloped; see the left panel of Fig. 2A. Except for the local extrema, points that lie on the left and right branches are stable fixed points of the fast equations (13). The threshold VT is chosen to intersect the middle branch of the cubic nullcline. The solution trajectory for this case is easy to understand. Starting at the local maximum of the left branch, equation (13) is used to make a fast jump to the right branch. Note that this jump is horizontal since ds/dτ = 0 according to (14). Then (11)–(12) are used to evolve the slow flow down the right branch until the trajectory reaches the local minimum. A fast jump back to the left branch under (13)–(14) then ensues, followed by slow evolution under (11)–(12) along the left branch back to the local maximum.

Figure 2
Synaptic and electrical connectivity

When the AB to INT1 inhibtion is present (gABI > 0), then a portion of the VL nullcline moves in phase space. In particular, when the AB to INT1 synapse is active, then VI decreases. In turn, through equation (5), sIL decreases causing the VL nullcline to move down in the phase space. However, since the AB to INT1 synapse is irrelevant when LG is active, only the left branch of the nullcline is affected. The left panel of Fig. 2B shows the LG trajectory when the AB to INT1 inhibition is present. The small depolarizations while the trajectory is on the left branch correspond to periodic disinhibition from the INT1 inhibition to LG that is itself created by the periodic inhibition of INT1 by AB. When the trajetory has evolved sufficiently far up the left branch to above the local maximum of the lower nullcline, the disinhibition allows LG to escape from the INT1 inhibition and become active. In this case, the period of the orbit is reduced since the time spent on both the left and right branches is reduced.

3.2 The effect of non-voltage dependent electrical coupling on the INT-LG generated rhythm

We next investigate the effect of adding electrical coupling to the network. First we consider the case when the electrical coupling is not voltage dependent. To do so, set vel = −100. Since VL > vel in this case, this causes n(VL) = 1 in equations (9)(10). The effect of gelec > 0 is to lower the VL nullcline in the phase space; see Fig. 3A. Note that because VM does not change and, for this case gelec(VL) is constant, the effect of the electrical current on the VL nullcline is largely due to the difference VLVM. This difference is the driving force of the electrical current. Since VM is constant, it acts like the driving force of a synaptic current that drives the voltage towards a constant reversal potential. When gelec > 0, the left branch of the VL-nullcline moves down more than the right branch since the driving force is larger there. That being said, the effect on the left branch is not too much larger than on the right branch. The result of the electrical coupling is simply to increase the burst duration of LG and shorten its interburst duration. The reason for this is readily explained through the phase plane of LG. The slow flow is directly related to the distance of the trajectory from the s-nullcline. When gelec > 0, the right branch of the nullcline moves down toward s = 0 thereby slowing the trajectory down when LG is active. The opposite happens to the left branch; the distance from the s-nullcline increases, thus speeding up the trajectory in the silent state. The period of LG is an increasing function of gelec. In fact, the period tends to infinity when gelec becomes sufficiently large as a saddle-node bifurcation at s = 0 is created.

Figure 3
The effect of non-voltage dependent electrical coupling

Next, observe that electrical coupling and the MCN1 synapse have similar effects on the VL-nullcline. Namely, increases in either gML or gelec lower the VL nullcline. This implies that some amount of the chemical synaptic excitation can be replaced by the metabolically less costly electrical coupling. For instance, begin with gelec = 0 and gML chosen such that the left branch of the VL-nullcline intersects s = 1 creating a stable fixed point; Fig 3B. If gelec is now chosen sufficiently large then the VL-nullcline is lowered enough so that the fixed point on the left branch moves to the middle branch and is unstable. However, if gelec is too large, then the right branch of the VL-nullcline intersects the s-nullcline at s = 0 creating an asymptotically stable fixed point there. Thus, there can exist a range (g*(gML), g*(gML)) of gelec values for which the fixed point lies along the middle branch and oscillations can occur. Note however if gML is too small, then the value gelec needed to move the local maximum below s = 1 would be so large that it would also lower the local minimum to below s = 0, creating a stable fixed point there. In these cases there is no range of gelec values that produce oscillations.

We can get a better understanding of the range of conductance values for which oscillations exist. Figure 3C shows a bifurcation diagram in gML-gelec space for the non-voltage dependent case. The shaded region R1 depicts the range of parameter values for which oscillations exist. Note that this region is bounded on three sides by lines. The lower boundary along gelec = 0 corresponds to the range of oscillations that exist when there is no electrical coupling. For this set of parameters, the boundary begins at roughly (8.91, 0). If gML < 8.91 and gelec = 0, then there are no oscillations as the VL-nullcline has a fixed point on its left branch at s = 1.

The left boundary corresponds to the set of saddle-node values along the local maximum of the VL-nullcline at s = 1. This curve is a line and has negative slope. To see why, consider the equation F (VL, s) = 0 and equation (15) for the VL nullcline in the voltage independent case where gelec(VL) = gelec. We rewrite (15) as follows


where f(VL) refers to the first two terms in the numerator on the left hand side of (15). A saddle-node point occurs when F (vL, 1) = 0 and ds/dVL = 0. The equation F (VL, 1) = 0 implies


Next observe that


The condition ds/dVL = 0 implies that the numerator of the above fraction equals zero which reduces to the relationship,


Let VL(g¯elec) denote the solution of (19) and note it that does not depend on gML. Further, it only weakly depends on gelec in the sense that this term is scaled by the difference VMEexc. Therefore the curve that defines the saddle-node points given in (17) is basically a line with the slope given by the ratio of the driving forces (VL-Eexc)/(VL-VM). Note that if Eexc = VM, then the slope of the saddle-node curve is negative one and the VL value of the local maximum is independent of both gML and gelec.

The top boundary of the oscillation region corresponds to the set of saddle-node points when the minimum of the cubic nullcline is tangent to s = 0. This curve is given by F (VL, 0) = 0 and ds/dVL = 0. From (19), we already know that the solution to the latter are independent of gML. Now from (16), the intersection of the VL nullcline with s = 0 is also independent of gML. Thus the top boundary is simply a horizontal line in the gML-gelec plane.

The region R1 is unbounded on the right. This is precisely because the local minimum of s at s = 0 is independent of gML. As gML → ∞, the oscillations are no longer burst-like. Instead the trajectory spends almost all of its time on the right branch in a neighborhood of the local minimum.

3.3 The effect of voltage dependent electrical coupling on the INT1 – LG generated rhythm

To explore the role of voltage dependence on the electrical coupling in the INT1-LG generated rhythm, we let vel = VT which is a value that lies along the middle branch of the VL-nullcline. The voltage dependence now allows the conductance of the electrical coupling to vary as a function of VL between gmin along the left branch of the VL-nullcline and gelec along the right branch. Thus the voltage-dependent electrical coupling affects the right branch of the VL nullcline much more than the left branch. This is in contrast to the non-voltage dependent case; compare Fig 4A and B.

Figure 4
The effect of voltage dependent electrical coupling between LG and the MCN1 output to LG

Figure 4C shows the regions of oscillations for these cases. For this set of parameters, there are two primary differences between the voltage-dependent (R2) and independent (R1) cases. First, the left boundary is more steeply sloped and the top boundary sits at a higher gelec value compared to the voltage-independent case. Both are easily explained. In the voltage-dependent case, equation (17) becomes


The condition ds/dVL = 0 yields a solutions VL(g¯elec) which is again independent of gML. By definition n(VL) < 1. Thus the prefactor multiplying gML is in fact a slope and is larger in magnitude than in the voltage-independent case. Thus the left boundary is steeper (−5.4 compared to −0.8 for the default parameters).

The intersection of the VL nullcline with s = 0 satisfies


The value VL increases with voltage dependence (specifically with vel from (10)). As a result, the right-hand side of (21) increases since the numerator increases while the denominator decreases. In the voltage independent case, n(VL) [equivalent] 1, whereas in the voltage dependent case n(VL)<1. To compensate, the maximal conductance of the electrical coupling gelec must increase. This allows the top boundary of the region R2 to sit at higher values of gelec (≈ 1.57 compared to 1.2 in the voltage independent case).

The effect of voltage dependence can be amplified by making the n(VL) curve less steeply sloped. For instance, if kel is increased from 5 to 20, then the slope of the left boundary of the oscillation region decreases in magnitude to around 2, which is much closer to the voltage-independent case; see R3 in Fig. 4D. Further, because the change in n is more gradual, larger values of gelec are needed to satisfy (21), so that the top boundary of R3 now sits around 2.02 compared to 1.57 for R2. Other changes of parameters can similarly be explored.

3.4 Oscillations arising through the voltage-dependent MCN1–LG coupling in the absence of the INT1-LG HCO

To this point, we have simply shown how electrical coupling affects the existing oscillations that arise through the INT1-LG HCO. A more important observation that we now make is that oscillations can arise in the absence of this HCO provided that the electrical coupling is voltage-dependent.

Consider equations (11)(14) with gIL = 0. This removes the INT1 to LG inhibition and destroys the HCO mechanism for oscillations. The VL-nullcline now is defined by


In this case, to see why voltage dependence is necessary for oscillations, first take the case where the electrical coupling is non-voltage dependent. Then ds/dVL = [grest,L[EexcErest,L] + gelec[EexcVM ]]/gML[VLEexc]2 > 0 if VM is not too large. In this case, the VL nullcline is a monotone increasing function that asympotes to −[grest,L + gelec]/gML as VL → −∞ and Eexc as VL → ∞; see Fig 5A. In this case, oscillations are not possible as any ensuing fixed point is asymptotically stable.

Figure 5
Oscillations arising through the voltage-dependent coupling between MCN1 and LG when the INT1–LG HCO is ineffective

Now take the case when the electrical coupling is voltage dependent. Then after some algebraic manipulation, the condition ds/dVL = 0 yields


For simplicity, take Eexc = VM in which case the condition reduces to


The left hand side is independent of gelec, while the right hand side increases with it. Further the right hand side has a zero at VL = VM and also tends to 0 as VL±∞. Thus for gelec sufficiently large, there are two solutions of (23), meaning that the graph of (22) has a local maximum and minimum. In this case, the VL-nullcline is again cubic shaped and oscillations are possible; see Fig 5B black trajectory and voltage trace. Therefore, voltage-dependent electrical coupling together with the slow excitation from MCN1, and its subsequent removal, via pre-synaptic inhibition from LG provides an alternate mechanism for the generation of oscillations. Note that the voltage range and period of the oscillation are within the range of the oscillation generated by the INT1 – LG HCO.

Using equation (10), we can derive an estimate on how large gelec needs to be to obtain oscillations. The right hand side of (23) has a local maximum at VL = Vel. Substituting and finding the smallest value of gelec that allows the right hand side to equal the left yields


This condition is fairly straightforward to interpret. Namely, the stronger the passive properties of LG, either through larger leak conductance grest,L or smaller leak reversal Erest,L, or the more gradual the voltage dependence, larger kel or vel, the larger the electrical conductance gelec needs to be.

We next explore the role of INT1 on the MCN1 – LG generated oscillation. We emphasize that, although the inhibition from INT1 to LG is restored, the parameters remain in range where the inhibition based HCO-based mechanism is not capable of producing oscillations. INT1 inhibition to LG raises the LB of VL-nullcline as shown in Fig 5B. Now the trajectory (black) must increase to higher values of s in the phase plane to escape inhibition, thereby increasing the interburst duration. In turn, when LG is active, the trajectory must also traverse through a larger range of s values to reach the local minimum of the cubic, thereby increasing LG’s burst duration. Thus the effect of this inhibition is to increase the oscillation period (and range of voltage values) by increasing both the interburst and burst duration (black voltage traces).

When AB to INT1 inhibition is included, the trajectory is allowed to leave the left branch prematurely at one of the moments in time when INT1 is inhibited by AB. This results in a shorter interburst and burst duration very similar to what was described in Section 3.1. Note that the period is very similar to that obtained when INT1 to LG inhibition is completely absent (gIL = 0); see Fig. 5C. This makes sense as the AB inhibition to INT1 has the practical effect of making gIL = 0 periodically when LG is in its interburst. Thus it is at one of those moments in time when LG is able to escape from inhibition.

3.5 Voltage-dependent oscillations in the Morris-Lecar equations

We now demonstrate that our main findings regarding the role of voltage dependent electrical coupling hold in a model in which LG and INT1 are modeled using biophysical equations. We model each of these cells using the two-dimensional Morris-Lecar equations, which are a commonly used set of equations that are derived in the Hodgkin-Huxley formalism. The voltage equation includes ionic currents for calcium, potassium and a leak current. There is a recovery variable associated with the activation of the potassium current. The equations for each cell are


On the right-hand side of equations (25) and (27), the first three terms are specific to the Morris-Lecar equations, while the remaining terms have the same form as defined in Section 2.1. The specific details of the model and parameter values are provided in the Appendix. Of interest to us here is the shape of the nullclines of the two cells. For INT1, parameters are chosen such that in the absence of input (gLI = 0, gABI = 0), the VI nullcline is cubic shaped and intersects the sigmoidal WI nullcline on its right branch. This high voltage fixed point indicates that INT1 is tonically active in the absence of input.

For LG, we consider two different parameter choices. In one case, in the absence of input, we choose gCa,L = 4.0 which is large enough so that the VL nullcline is cubic shaped. In that case, the VL and WL nullclines intersect along the left branch of VL which models LG being at rest. In this case, LG is excitable in the classical sense that, if it receives the appropriate synaptic input from MCN1, it will fire. With these parameters, and with the LG-INT1 HCO intact, the presynaptic LG to MCN1 inhibition is sufficient to produce oscillations, as was shown in the previous sections. The addition of electrical coupling, either voltage-dependent or not, simply modulates the oscillations in a manner analogous to that found in Sections 3.2 and 3.3. In other words, electrical coupling is not necessary to produce oscillations. Numerical simulations (not shown) in this case yield results that are qualitatively similar to those found in Figs. 3 and and44.

The more interesting situation arises in the second case when we choose gCa,L = 0.5 so that the VL nullcline is monotone decreasing. Now, LG is no longer excitable. As a result, the INT1-LG HCO is not able to produce oscillations, independent of whether gIL is zero or not. Just as in Section 3.4 with the simple model, voltage dependent electrical coupling can provide an alternative mechanism for oscillations. Figure 6 shows results from the biophysical model. First consider the case when gIL = 0. The phase plane in Fig. 6A shows the projection of VL nullcline for four different cases onto the VLWL phase plane when the electrical coupling is independent of voltage. The solid curves are for gelec = 2.2 where brown corresponds to s = 0 and red is s = 1. The dashed curves are their counterparts for gelec = 22. Because the MCN 1 to LG excitation which is governed by s can change slowly, the four nullclines that are shown are only representative snapshots of the VL nullcline. However for all values of s, the VL nullcline is monotone decreasing, precluding the possibility of oscillations.

Figure 6
Oscillations arising through the voltage-dependent coupling between MCN1 and LG in the biophysical model

In contrast, consider Fig. 6B1. Shown is the VL nullcline when the electrical coupling is voltage-dependent for two different values of s (smaller s in brown, larger in red). As can be seen, the voltage dependence creates a cubic shaped nullcline by preferentially affecting the nullcline at higher voltages. As a result oscillations are possible. The VL trajectory is superimposed on the figure. The red nullcline associated with the larger value of s corresponds to those at which the trajectory jumps from the left branch to the right branch signaling LG’s transition to the active state. The brown nullcline is associated with a smaller value of s when the LG trajectory jumps from the right branch to the left branch signaling LG’s transition to the silent state. The dependence on s is seen in panel B2 which shows the projection of nullclines and the trajectory onto the VL versus s phase space; note the parallel to Fig. 5B. Recall that s increases when LG is in the silent phase. This means that in the VLs phase plane, the trajectory moves up along the left branches. However, the left branches themselves are moving down because as s increases, the added excitation from MCN1 produces a greater chance to become active. The jump to the active state occurs from a local maximum of the red nullcline. On the right branch, the trajectory moves down, but the nullcline moves up. The jump to the silent state occurs from the minimum of the brown nullcline. The corresponding voltage traces for both LG and INT1 are shown in panel B3.

In Fig. 6C, we restore the INT1 to LG synapse gIL = 10.The LG interburst length increases, as was also seen in the simple model Fig. 5B. As before, this is because the inhibition from INT1 to LG means that s has to increase to larger values for LG to jump to the active state. This implies a longer interburst duration. Finally, in Fig. 6C, we restore the AB input to INT1 which shortens the LG burst and interburst in a similar manner to Fig. 5C because the periodic inhibition of INT1 by AB provides periodic disinhibition of LG. This provides LG an opportunity to escape the silent state earlier just as with the simple model, thereby shortening LG’s interburst and speeding up the rhythm.

Just as in Section 3.4, we can determine conditions under which voltage-dependence allows the electrical coupling to produce oscillations. Consider the case gIL = 0. For compactness of notation, define f(VL) = −grest,L[VLErest,L] − gCa,Lm(VL)[VLECa], h(VL) = −gelecn(VL)[VLEM] − gss[VLEexc] + Iapp,L. Note that s depends on VL. Let prime denote the derivative with respect to VL. We can solve for the VL nullcline by setting the right-hand side of equation (25) to zero and solving for WL.


The slope of this nullcline is given by


To show that the VL nullcline can be cubic shaped, we need to find conditions under which the derivative (29) changes sign. Observe that f(VL) + h(VL) > 0 thus the second term in the numerator of (30) is negative. Next observe that


if gCa,L is sufficiently small. The derivative


The first term in (32) is non-negative, while the remaining three are all negative (note that s′(VL) < 0). Thus the sign of h′(VL) will be negative unless the first term is sufficiently large. When the electrical coupling is not dependent on voltage, n(VL) [equivalent] 1 and therefore n(VL)=0. Thus the first term will actually be zero in this case. In turn this implies that h′(VL) < 0. Together, these results imply that, in the voltage-independent case, the VL nullcline remains monotone decreasing and that no oscillations are possible. Alternatively, when the electrical coupling is voltage dependent, then n(VL)>0 and is also relatively large in an intermediate range of VL values (roughly −35 to −15 mv). Thus for gelec sufficiently large, the first term in (32) dominates the others and h′(VL) > 0. As a consequence, if gelec is large enough, then voltage dependent electrical coupling allows dWL/dVL > 0 over a range of intermediate VL values. In conjunction with the synaptic input from MCN1, this provides the opportunity for oscillations to exist.

4 Discussion

Neuronal circuits involved in the generation of rhythmic behavior often involve half center oscillators that are composed of sets of reciprocally inhibitory neurons. There is an extensive and ongoing effort to understand the dynamics of half center oscillators in the context of central pattern generation [10, 11, 14, 16, 17]. In many cases, it has been noted that a careful coordination between network elements is necessary to generate and set the frequency of the network [18, 19, 20]. The role of electrical coupling in rhythmic networks has also been studied [21, 22] where the neurons were modeled as intrinsic oscillators. Electrical coupling was not needed to generate oscillations, but rather used to modulate the characteristics of the oscillation.

As part of a larger work on the role of feedback to projection neurons, Kintos and colleagues [10, 12] had shown how to employ phase plane analysis to understand the effect of MCN1 synaptic input on the GMR. In particular, they showed how to analyze MCN1 synaptic input and AB inhibition of INT1 to determine the frequency of the GMR. In this paper, we have extended this analysis to show how to incorporate the effect of MCN1 – LG voltage-dependent electrical coupling to determine the conditions under which electrical coupling in the absence of the LG – INT1 HCO can generate oscillations.

In the presence of an intact LG – INT1 HCO, we first considered the effect of non-voltage dependent electrical coupling. We showed that the non-voltage dependent electrical coupling acts to increase the LG burst duration while shortening its interburst duration. This occurs because the voltage of LG is driven towards the fixed, large voltage of MCN1. If the strength of the electrical coupling is too large, however, LG gets stuck in its burst phase. One advantage of the non-voltage dependent electrical coupling is that it can be used in conjunction with the MCN1 chemical synapse allowing for the generation of the GMR for a smaller amount of the chemical excitation. This is a ”cheaper” way to generate oscillations as it requires less synaptic resources. The bifurcation diagram in Fig. 3C shows the precise relationship between electrical and synaptic coupling needed to create oscillations. We showed that boundaries of this diagram are all roughly linear. In the case of voltage dependent electrical coupling, the right branch of the LG nullcline is affected much more significantly than the left branch. This allows for an increase in the LG burst duration and a larger range of values of gelec for the generation of network oscillations.

A significant finding of our study is that network oscillations can also be generated in the absence of coupling between LG and INT1 simply through the voltage dependent electrical coupling between MCN1 and LG and the slow excitation from MCN1, together with its removal due to the pre-synaptic inhibition of this excitation. We derived a condition on the minimum value of gelec in order for the GMR oscillations to exist in the absence of the HCO. We showed that non-voltage dependent electrical coupling alone is not sufficient for generation of the GMR. When the reciprocal inhibition between LG and INT1 is restored, the period of the oscillations increases due to increases in both the interburst and burst durations of the oscillations. If, in addition, AB periodically inhibits INT1, the interburst duration of LG is shortened. This is a direct result of the disinhibitory effect of LG from INT1 each time AB fires.

Our findings are not limited to the simple model in which LG and INT1 are modeled as passive cells that we first considered. We showed that voltage-dependent electrical coupling played the same role in a model in which these cells were described using the biophysically based Morris-Lecar equations. In order for voltage-dependent electrical coupling to create the mechanism for oscillations, we showed that LG must not be modeled as being excitable. This fact is consistent with the underlying biological properties of the LG neuron, which, in the absence of MCN1 or other modulatory input, shows no active properties (e.g. post-inhibitory rebound, voltage sags or plateaus) that are associated with slow bursting oscillations [4].

There are several natural extensions of this work. In previous work [23], based on experiments of Wood et. al. [24], we showed that AB inhibition toMCN1 provides an alternate mechanism to regulate the gastric mill frequency. In the current work, we did not include the inhibition from AB to MCN1. If the AB inhibition to MCN1 were included, the LG burst would end when AB inhibits MCN1. It would be necessary for MCN1 to be gated when LG is in its active state in order to maintain robust oscillations. Indeed, in the VCN-activated version of the gastric mill rhythm, the AB to MCN1 synapse is gated out during LG active phase [25]. It would be of interest to extend our current model to test whether this gating is truly necessary to maintain oscillations.

Another area that remains to be explored is the role of electrical coupling in the MCN1/CP N2 generated gastric mill rhythm. Kintos and Nadim [10] showed that the LG – INT1 HCO could be replaced by a tri-synaptic pathway that included the projection neuron CP N2. Of interest would be to see whether voltage dependence can replace one or more of those synaptic pathways.

Although the networks under consideration in this, and related papers, are relatively simple and only involve a small number of neurons, it is evident that the dynamics exhibited by them can be quite complicated. Moreover, the neural mechanisms that underlie the existence of oscillations are often hard to separate from those that simply modulate the rhythmic properties of these networks. Minimal modeling and mathematical analysis of small networks plays a critical role in allowing us to discern which inputs generate oscillations versus those that modulate oscillations by providing valuable insights into how these important central pattern generating networks operate.


  • A simple model is used to understand how electrical coupling affects the ability of oscillations to be created through a standard half-center oscillator mechanism is incapable of doing so.
  • Voltage-dependent electrical coupling is shown to provide an alternate mechanism for the generation of oscillations when an inhibition based half-center oscillator mechanism is incapable of doing so.
  • This same result is then demonstrated in a biophysical model based on the Morris-Lecar equations.
  • Conditions on the parameters in both models are derived that show why electrical coupling must be voltage dependent to produce oscillations.


This work was supported in part by the National Science Foundation, DMS 1122291 (AB, FN) and the National Institutes of Health, NIH R01-MH060605 (FN).

List of Abbreviations

central pattern generating
gastric mill rhythm
stomatogastric ganglion
half-center oscillator
lateral gastric
interneuron 1
stomatogastric nervous system
anterior burster
modulatory commissural neuron 1

5 Appendix

Numerical simulations were performed using XPPAUT [26]. For the simple model of passive cells used to produce Figs. 25 the following set of equations was used.


The term VI = h1(VL, sABI) that appears in equation (33) is goverened by


Note here the presence of the function P(VL)=(1+exp(VL-v3k3))-1. This term is used to gate out the effect of AB input to INT1, and its subsequent effect on the VL nullcline when INT1 is in its silent state. The other remaining equations are simply (3–5) sABI(t)=Heav(sin(2πt)1000-0.5),sLI(VL)=(1+exp(v1-VLk1))-1,sIL(VI)=(1+exp(v2-VIk2))-1 and (9) and (10) written as one gelec(VL)=g¯elec[(1-gmin)(1+exp(vel-VLkel))-1+gmin].

Table 1 shows parameter values that were common to all simulations of the simple model. Below that we show specific values used for gelec and gML for each of the relevant figures.

Table 1
Parameter values common to all simulations of the simple model

For Figures 2 to to4,4, we chose gIL = 12. For Figure 2: gML = 10, gelec = 0, gABI = 0 (2A) and gABI = 0.2 (2B). For Figure 3A: gML = 10, gABI = 0 and gelec = 0, 0.5, 1.0, 1.5. For Figure 3B: gML = 8.8 and gelec = 0.0 (upper cubic) and 0.8 (lower cubic). For Figure 4A, the electrical coupling is non-voltage dependent: gML = 8.8 and gelec = 0.088, gelec = 0.155 at the left saddle node point, and gelec = 1.2 at the right saddle node point. For Figure 4B, the electrical coupling is voltage dependent and we chose, gML = 8.8 and gelec = 0, gelec = 0.594 at the left saddle node point, and gelec = 1.57 at the right saddle node point.

For Figure 5A: gML = 0.35, gIL = 0, the three monotone nullclines are when the electrical coupling is non-voltage dependent with gelec = 0, gelec = 0.6 and gelec = 1.3. The cubic is for the voltage-dependent case with gelec = 1.3. For Figure 5B: voltage dependent electrical coupling, gML = 0.35, gelec = 1.24, gIL = 0 (lower cubic) or gIL = 0.2 (upper cubic). For Figure 5C: same as Fig. 5B except gABI = 0.2.

For the simulations shown in Fig. 6, the following set of equations was used:


The synaptic variables are governed by


The remaining terms are given by m(Vx)=(1+tanh(Vx-cv1cv2))/2,τ(Vx)=cosh(Vx-cv32cvv),w,x(Vx)=(1+tanh(Vx-wfxcv4))/2, where the subscript x refers to either L or I. In addition, we used equations (3), (9) and (10). For Fig. 6A–B,gIL = 0, for Fig. 6C, gIL = 10 and for Fig. 6D, gABI = 1. Table 2 shows parameter values for the simulations of the Morris-Lecar model.

Table 2
Parameter values common to all simulations of the Morris-Lecar model


Conflict of Interest

The authors declare that they have no conflict of interest.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Sneyd J, Keizer J, Sanderson MJ. Mechanisms calcium oscillations and waves: a quantitative analysis. The FASEB Journal. 1995;9(14):1463–1472. [PubMed]
2. Marder E, Bucher D. Central pattern generators and the control of rhythmic movements. Current Biology. 11(23):R986–R996. [PubMed]
3. Nusbaum MP, Beenhakker MP. A small-systems approach to motor pattern generation. Nature. 2002;417:343–350. [PubMed]
4. Coleman M, Nusbaum MP. Functional consequences of compartmentalization of synaptic input. Journal of Neuroscience. 1994;14:6544–6552. [PubMed]
5. Coleman M, Meyrand P, Nusbaum MP. A switch between two modes of synaptic transmission mediated by presynaptic inhibiton. Nature. 1995;378:502–505. [PubMed]
6. Calabrese RL. Oscillations in motor pattern-generating networks. Current Opinion Neurobiology. 1995;5:816–823. [PubMed]
7. Skinner F, Kopell N, Marder E. Mechanisms for oscillation and frequency control in reciprocally inhibitory model neural networks. Journal of Computational Neuroscience. 1994;1(1–2):69–87. [PubMed]
8. Nadim F, Manor Y, Nusbaum MP, Marder E. Frequency regulation of a slow rhythm by a fast periodic input. Journal of Neuroscience. 1998;18:5053–5067. [PubMed]
9. Ambrosio C, Bose A, Nadim F. The effect of modulatory neuronal input on gastric mill frequency. Neurocomputing. 2005;65–66:626–631.
10. Kintos N, Nadim F. A Modeling exploration of how synaptic feedback to descending projection neurons shapes the activity of an oscillatory network. SIAM J Applied Dynamical Systems. 2014;13(3):1239–1269. [PMC free article] [PubMed]
11. Kintos N, Nusbaum MP, Nadim F. A Modeling comparision of projection neuron-and neuromodulator-elicited oscillations in a central pattern generator network. J Comput Neurosci. 2008;24:374–97. [PMC free article] [PubMed]
12. Kintos N, Nusbaum MP, Nadim F. A Modeling comparision of projection neuron-and neuromodulator-elicited oscillations in a central pattern generator network. J Comput Neurosci. 2016;40:113–135. [PMC free article] [PubMed]
13. Manor Y, Nadim F, Epstein S, Ritt J, Marder E, Kopel N. Network oscillations generated by balancing graded asymmetric reciprocal inhibition in passive neurons. Journal of Neuroscience. 1999;19(7):2765–2779. [PubMed]
14. Reyes M, Carelli P, Sartorelli J, Pinto R. A modeling approach on why simple central pattern generators are built of irregular neurons. PLoS ONE. 2015;10(3):e0120314. [PMC free article] [PubMed]
15. Morris C, Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophys J. 1981;35:193–213. [PubMed]
16. Weeks J. Neuronal basis of leech swimming. Journal of Neuroscience. 1981;45(4):698–723. [PubMed]
17. Zhang C, Lewis T. Phase response properties of half-center oscillators. Journal of Computational Neuroscience. 2013;35(1):55–74. [PubMed]
18. Cymbalyuk G, Gaudry Q, Masion M, Calabrese R. Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J Neurosci. 2002;22(24):10580–10592. [PubMed]
19. 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(2):169–187. [PMC free article] [PubMed]
20. Szucs A, Selverston A. Consistent dynamics suggests tight regulation of biophysical parameters in a small network of bursting neurons. J Neurobiol. 2006;66(14):1584–1601. [PubMed]
21. Chow C, Kopell N. Dynamics of spiking neurons with electrical coupling. Neur Comp. 2000;12:1643–1678. [PubMed]
22. Lee E, Terman D. Stable antiphase oscillations in a network of electrically coupled model neurons. SIAM J Appl Dyn Syst. 2013;12(1):1–27.
23. Ambrosio-Mouser C, Farzan N, Bose A. The effects of varying the timing of inputs on a neural oscillator. SIAM Journal on Applied Dynamical Systems. 2006;5(1):108–139. [PMC free article] [PubMed]
24. Wood D, Stein W, Nusbaum M. Projection neurons with shared cotransmitters eleicit different motor patterns from the same neural circuit. J Neuroscience. 2008;20(23):8943–8953. [PubMed]
25. Blitz DM, Nusbaum MP. State-dependent presnyaptic inhibition regulates central pattern generator feedback to descending inputs. J Neurosci. 2008;28:9564–9674. [PMC free article] [PubMed]
26. Ermentrout B. Simulating, analyzing, and animating dynamical systems : a guide to XPPAUT for researchers and students. Society for Industrial and Applied Mathematics; Philadelphia: 2002.