Key Bifurcations of Bursting Polyrhythms in 3-Cell Central Pattern Generators

We identify and describe the key qualitative rhythmic states in various 3-cell network motifs of a multifunctional central pattern generator (CPG). Such CPGs are neural microcircuits of cells whose synergetic interactions produce multiple states with distinct phase-locked patterns of bursting activity. To study biologically plausible CPG models, we develop a suite of computational tools that reduce the problem of stability and existence of rhythmic patterns in networks to the bifurcation analysis of fixed points and invariant curves of a Poincaré return maps for phase lags between cells. We explore different functional possibilities for motifs involving symmetry breaking and heterogeneity. This is achieved by varying coupling properties of the synapses between the cells and studying the qualitative changes in the structure of the corresponding return maps. Our findings provide a systematic basis for understanding plausible biophysical mechanisms for the regulation of rhythmic patterns generated by various CPGs in the context of motor control such as gait-switching in locomotion. Our analysis does not require knowledge of the equations modeling the system and provides a powerful qualitative approach to studying detailed models of rhythmic behavior. Thus, our approach is applicable to a wide range of biological phenomena beyond motor control.


Introduction
A central pattern generator (CPG) is a circuit of neuronal cells whose synergetic interactions can autonomously produce rhythmic patterns of activity that determine vital motor behaviors in animals [1][2][3]. CPGs have been found in many animals, where they have been implicated in the control of diverse behaviors such as heartbeat, sleep, respiration, chewing, and locomotion on land and in water [4][5][6][7]. Mathematical modeling studies, at both abstract and realistic levels of description, have provided useful insights into the operational principles of CPGs [8][9][10][11][12][13][14]. Although many dynamic models of specific CPGs have been developed, it remains unclear how CPGs achieve the level of robustness and stability observed in nature [15][16][17][18][19][20][21][22].
A common component of many identified CPGs is a half-center oscillator (HCO), which is composed of two bilaterally symmetric neurons reciprocally inhibiting each other to produce an alternating anti-phase bursting pattern [23]. There has been much work on the mechanisms giving rise to anti-phase bursting in relaxation HCO networks, including synaptic release, escape and post-inhibitory rebound [24,25]. Studies of HCOs composed of Hodgkin-Huxley type model cells have also demonstrated the possibility of bistability and the coexistence of several in-phase and anti-phase bursting patterns based on synaptic time scales or delays [26][27][28].
We are interested in exploring the constituent building blocks -or ''motifs'' -that may make up more complex CPG circuits, and the dynamic principles behind stable patterns of bursting that may co-exist in the circuit's repertoire of available states [13,20,29]. We will refer to such multi-stable rhythmic patterns as ''polyrhythms.'' We consider the range of basic motifs comprising three biophysical neurons and their chemical synapses, and how those relate to, and can be understood from the known principles of two-cell HCOs. We will study the roles of asymmetric and unique connections, and the intrinsic properties of their associated neurons, in generating a set of coexisting synchronous patterns of bursting waveforms. The particular kinds of network structure that we study here reflect the known physiology of various CPG networks in real animals. Many anatomically and physiologically diverse CPG circuits involve a three-cell motif [30,31], including the spiny lobster pyloric network [1,32], the Tritonia swim circuit, and the Lymnaea respiratory CPGs [33][34][35][36].
An important open question in the experimental study of real CPGs is whether they use dedicated circuitry for each output pattern, or whether the same circuitry is multi-functional [37,38], i.e. can govern several behaviors. Switching between multi-stable rhythms can be attributed to input-dependent switching between attractors of the CPG, where each attractor is associated with a specific rhythm. Our goal is to characterize how observed multistable states arise from the coupling, and also to suggest how real circuits may take advantage of the multi-stable states to dynamically switch between rhythmic outputs. For example, we will show how motif rhythms are selected by changing the relative timing of bursts by physiologically plausible perturbations. We will also demonstrate how the set of possible rhythmic outcomes can be controlled by varying the duty cycle of bursts, and by varying the network coupling both symmetrically and asymmetrically [17,20]. We also consider the role of a small number of excitatory or electrical connections in an otherwise inhibitory network. Our greater goal is to gain insight into the rules governing pattern formation in complex networks of neurons, for which we believe one should first investigate the rules underlying the emergence of cooperative rhythms in smaller network motifs.
In this work, we apply a novel computational tool that reduces the problem of stability and existence of bursting rhythms in large networks to the bifurcation analysis of fixed points (abbreviated FPs) and invariant circles of Poincaré return maps. These maps are based on the analysis of phase lags between the burst initiations in the cells. The structure of the phase space of the map reflects the characteristics the state space of the corresponding CPG motif. Equipped with the maps, we are able to predict and identify the set of robust bursting outcomes of the CPG. These states are either phase-locked or periodically varying lags corresponding to FP or invariant circle attractors (respectively) of the map. Comprehensive simulations of the transient phasic relationships in the network are based on the delayed release of cells from a suppressed, hyperpolarized state. This complements the phase resetting technique and allows a thorough exploration of network oscillations with spiking cells [39]. We demonstrate that synapticallycoupled networks possess stable bursting patterns that do not occur in similar motifs with gap junction coupling, which is bidirectionally symmetric [40].

Results
Our results are organized as follows: first, we describe our new computational tools, which are based on 2D return maps for phase lags between oscillators. This is a non-standard method that has general utility outside of our application, and we therefore present it here as a scientific result. We then present maps for symmetric inhibitory motifs and examine how the structure of the maps depends on the duty cycle of bursting, i.e. on how close the individual neurons are to the boundaries between activity types (hyperpolarized quiescence and tonic spiking). Here, we also examine bifurcations that the map undergoes as the rotational symmetry of the reciprocally coupled 3-cell motif is broken. This is followed by a detailed analysis of bifurcations of fixed point (FP) and invariant circle attractors of the maps, which we show for several characteristic configurations of asymmetric motifs, including a CPG based on a model of the pyloric circuit of a crustacean. We conclude the inhibitory cases with the consideration of the fine structure of a map near a synchronous state. We then discuss the maps for 3-cell motifs with only excitatory synapses, which is followed by the examination of mixed inhibitory-excitatory motifs, and finally an inhibitory motif with an additional electrical synapse in the form of a gap junction.

A computational method for phase lag return mappings
We first introduce the types of trajectories we focus on and how we measure them. The reduced leech heart interneuron can demonstrate many regular and irregular activity types, including hyper-and de-polarized quiescence, tonic spiking and bursting oscillations. We focus on periodic bursting, and Figure 1 shows a trajectory (dark gray) in the 3D phase space of the model. The helical coils of the trajectory correspond to the active tonic spiking period of bursting due to the fast sodium current. The flat section corresponds to the hyperpolarized quiescent portion of bursting due to the slow recovery of the potassium current. In Fig. 1, two snapshots (at t~0 and t~10 s) depict the positions of the blue, green and red spheres representing the momentarily states of all three interneurons. The coupling between the cells is chosen weak so that network interactions should only affect the relative phases of the cells on the intact bursting trajectory, i.e. without deforming the trajectory.
V shift K2 is a model parameter that measures the deviation from the half-activation voltage V 1=2~{ 0:018 V of the potassium channel, m ? K2~1 =2. We use V shift K2 as a bifurcation parameter to control the duty cycle (DC) of the interneurons. The duty cycle is the fraction of the burst period in which the cell is spiking, and is a property known to affect the synchronization properties of coupled bursters [16,17]. The individual cell remains bursting within the interval V shift K2 [½{0:024235,{0:01862. At smaller values of V shift K2 , it begins oscillating tonically about the depolarized steady state, and becomes quiescent at greater values of V shift K2 . Therefore, the closer the cell is to either boundary, the DC of bursting becomes longer or shorter respective: the DC is about 80% at V shift K2~{ 0:0225 V and 25% at V shift K2~{ 0:01895 V. For 50% DC we set V shift K2~{ 0:021 V, in the middle of the bursting interval (see Fig. 2).
When an isolated bursting cell is set close to a transition to either tonic spiking or hyperpolarized quiescence, its network dynamics become sensitive to external perturbations from its presynaptic cells. For example, when the post-synaptic cell is close to the tonic-spiking boundary, excitation can cause the post-synaptic cell to burst longer or even move it (temporarily) over the boundary into the tonic spiking (TS) region. In contrast, inhibition shortens the duty cycle of the post-synaptic neuron if it does not completely suppress its activity (Fig. 3).
The return map for phase lags. We reduce the problem of the existence and stability of bursting rhythmic patterns to the bifurcation analysis of fixed points (FPs) and invariant curves of Poincaré return maps for phase lags between the neurons. In this study, we mostly consider relatively weakly coupled motifs, but our approach has no inherent limitation to weak coupling. Here, the weakly coupled case is a pilot study that lets us test our technique and also uncover all rhythms, both stable and unstable, that can possibly occur in the network. Detailed scrutiny of the return maps is computationally expensive: an exploration of one parameter set can take up to three hours on a state-of-the-art desktop workstation depending on the accuracy of the mesh of initial conditions and length of the transients computed.
The phase relationships between the coupled cells are defined through specific events, t (n) 1 ,t (n) 2 ,t (n) 3 n o , when their voltages cross a threshold, H th , from below. Such an event indicates the initiation of the n th burst in the cells, see Fig. 4. We choose H th~{ 0:04 V, above the hyperpolarized voltage and below the spike oscillations within bursts. We define a sequence of phase lags by the delays in burst initiations relative to that of the reference cell 1, normalized over the current network period or the burst recurrent times for the reference cell, as follows: An ordered pair, M n~D w (n) 21 ,Dw (n) 31 , defines a forward iterate, or a phase point, of the Poincaré return map for the phase lags: , yields a forward phase lag trajectory, M n f g N n~0 , of the Poincaré return map on a 2D torus ½0,1)|½0,1) with phases defined on mod 1 (Fig. 5). Typically, such a trajectory is run for N~90 bursting cycles in our simulations. The run can be stopped when the distance between several successive iterates becomes less than some preset value, say DDM n {M nzk DDv10 {3 and k~5. This is taken to mean that the trajectory has converged to a fixed point, M Ã , of the map. This FP corresponds to a phase locked rhythm and its coordinates correspond to specific constant phase lags between the cells. By varying the initial delays between cells 2 and 3 with respect to the reference cell 1, we can detect any and all FPs of the map and identify the corresponding attractor basins and their boundaries.
We say that coupling is weak between two cells of a motif when the convergence rate to any stable FP of the return map is slow. This means that the distance between any two successive iterates of a trajectory of the return map remains smaller than some -like, resp., as well as an electrical connection through the gap junction between interneurons 1 and 2. (B) Bursting trajectory (gray) in the 3D phase space of the model, which is made of the ''active'' spiking (solenoid-like shaped) and the flat hyperpolarized sections. The gap between the 2D slow nullcline, m 0 K2~0 , and the low knee on the slow quiescent manifold, M eq , determines the amount of inhibition needed by the active pre-synaptic cell above the synaptic threshold, H syn , to either slow or hold the postsynaptic cell(s) at a hyperpolarized level around {0:06 V. The positions of the red, green and blue spheres on the bursting trajectory depict the phases of the weakly-connected cells of the CPG at two instances: the active red cell inhibits, in anti-phase, the temporarily inactive green and blue cells at two time instances. doi:10.1371/journal.pone.0092918.g001 is varied. Burst duration increases as V shift bound, e.g. max DDM n {M nz1 DDv0:05. Therefore, we can say that coupling is relatively strong if a remote transient reaches a FP of the map after just a few iterates. We point out that the convergence can be quick even for nominally small g syn provided that an individual cell is sufficiently close to either boundary of bursting activity (tonic spiking or quiescence).
We now make some technical remarks concerning computational derivations of the map, P. A priori, the initial period (recurrence time) of the motif's dynamics is unknown due to the unknown outcome of nonlinear cell interactions; furthermore, it varies over the course of the bursting transient until it converges to a fixed value on the phase locked state. Thus, we control the initial phases between the reference cell and the others releasing the latter from inhibition at various delays. To do this, we first estimate the initial phase lag with a first order approximation,   0ƒDw j1 v1 È É . Thus, we can set the initial phase lags by releasing the reference cell and keeping the others suppressed for durations t 12~D w (0) 21 T synch and t 13~D w (0) 31 T synch from the same initial point on the synchronous bursting trajectory, given by To complete a single phase lag map we choose the initial phase lags to be uniformly distributed on a grid of at least 40|40 points over the ½0,1)|½0,1) torus. The initial guess for the phase lag distribution is based on knowledge of a trajectory for a synchronized motif, whose period is already known. This guess will therefore differ from the self-consistent phase lag distribution once the networked cells begin to interact, especially with coupling strength variations. Similarly, the estimated network period, T synch , will differ from the network's actual self-consistent period. In computations, this may result in fast jumps from the set of guessed initial phases from n~0 to n~1. These jumps are artifacts of our setup and not relevant to our study of the attractors, and so we begin recording the phase lag trajectory settled from the second bursting cycle. Due to weak coupling, transients do not evolve quickly, and we connect phase lag iterates of the map by straight  lines in order to demonstrate and preserve the forward order, making them superficially resemble continuous-time vector fields in a plane. Lastly, we unfold the torus onto a unit square for the sake of visibility.
Note on interpreting phase lag diagrams. We use a consistent labeling convention to make our diagrams of the phase lag maps easy to interpret. In the first presentation of such a map ( Fig. 6 in the next section), we annotate the diagram with arrows to ð Þ(first episode), 3\f1E2g ð Þ(second episode) and 2\f1E3g ð Þ(third episode) in the short duty cycle motif with 25% DC with +5% random perturbations applied to all inhibitory connections with g syn~0 :0005. Switching between rhythms is achieved by the application of appropriately-timed hyperpolarized pulses that release the targeted cells. doi:10.1371/journal.pone.0092918.g007  show the directions of the map on successive forward iterates. We also label the position of FPs with colored dots. Larger dots are used for the stable FPs than the saddle points. The colors of the stable FPs (red, green, and blue) correspond to the colors of the computed phase lag trajectories that approach them (thereby depicting the basins of attraction). In Fig. 6, we indicate the directions of forward iterates of that map to assist the reader in their interpretation. However, all subsequent figures depict dynamics with these same essential interpretable features: all colored trajectories flow towards their color-coordinated FP, and small dots indicate saddle points. We also note that the origin in our maps has a complex fine structure but acts globally as a repeller. As such, we do not depict those FPs explicitly in the fullscale phase lag diagrams. A later section explicitly examines the fine structure near the origin.

Multistability and duty cycle in homogenous inhibitory motifs
We first examine three homogeneous (permutationally symmetric) configurations of the network with nearly identical cells and connections. We demonstrate that these symmetric network motifs are multistable and hence able to produce several coexisting bursting patterns. The homogeneous case allows us to reveal the role of the duty cycle as an order parameter that determines what robust patterns are observable. We suggest a biologically plausible switching mechanism between the possible bursting patterns by application of a small hyperpolarized current that temporarily blocks a targeted cell.
Short duty cycle motif. We begin with a weakly coupled with g syn~5 |10 {4 , homogeneous motif with 25% DC and V shift K2~{ 0:01895 V, which is close to the transition boundary between bursting and hyperpolarized quiescence. Proximity to the   boundary means that even weak inhibition is able to suppress a postsynaptic cell that is near the hyperpolarized quiescent state ( Figs. 1 and 3). Figure 5A shows the transient behaviors of the iterates of the phase lags Dw (n) 21 and Dw (n) 31 (shown in blue and gray colors) arising from initial conditions distributed uniformly over the unit interval. The phase lags exponentially converge to phase-locked states near 0 and 1 2 . Using Eq. (2), we compute the map P that is shown in Fig. 6B. The projection of the map onto the unit square is an efficient way to represent the synchronized evolution of the phase lags and facilitates easy identification of phase-locked states. These states are identified by three coexisting stable FPs or attractors of the system to which all forward iterates converge. Here, the FPs are: The robustness of a rhythm to perturbations is related to the size of its attractor basin. Similarly, FPs that have much larger basins than others can be thought of as ''dominating'' the phase plane in terms of likelihood of becoming the active state for a random initial condition or perturbation. Two triplets of saddles surround two more unstable FPs located at approximately 2 3 , 1 3 À Á and 1 3 , 2 3 À Á . The immediate neighborhood of the origin has a complex structure involving several FPs packed closely together, but globally it acts as a repeller (see later section for a detailed analysis of this).
Let us interpret the role of a stable FP, for example the red one, in terms of phase-locked bursting patterns. Since the phases are defined modulo one, the coordinates (Dw 21 ,Dw 31 )~0, 1 2 À Á , imply that the corresponding rhythm is characterized by two fixed conditions w 1~w2 and w 1 {w 3~1 2 . In other words, the reference cell fires in-phase with cell 2 and in anti-phase with cell 3. Symbolically, we will use the following notation for this rhythm: 3\f1E2g ð Þ , in which in-phase and anti-phase bursting are represented by (Dw 12~0 , or E) and (Dw 13~1 2 or \), respectively.
Following this notation, the stable FP (blue) at 1 2 , 1 2 À Á corresponds to the robust 1\f2E3g ð Þpattern, while the stable (green) FP 1 2 ,0 À Á corresponds to the 2\f1E3g ð Þpattern. These coexisting bursting rhythms are shown in Fig. 7. The motif can be made to switch between the polyrhythms by applying external pulses of appropriate duration to the targeted cells.
Two FPs around Dw 21 ,Dw 31 ð Þ 1 3 , 2 3 À Á and 2 3 , 1 3 À Á correspond to clockwise and counter-clockwise traveling waves (respectively) that Here, the period of either traveling wave is broken into three episodes in which each cell is actively bursting one at a time. For example, in Fig. 4 for the clockwise bursting, 1[2[3 ð Þ , the cell ordering is 1-2-3 before the pattern repeats. The co-existence of these two waves originates from the rotational symmetry of the homogeneous motif. However, both such traveling bursting waves are not robust and therefore cannot be observed in the motif with a short duty cycle because the corresponding FPs are repelling, so that a small perturbation will cause the phase lags of such a traveling rhythm to transition to those corresponding to one of three ''pacemaker'' states, as shown in Fig. 7.
Long duty cycle motif. Next, we consider the bursting motif that has a longer duty cycle of 80%, set by V shift K2~{ 0:0225 V. This brings the cells closer to the boundary separating bursting and tonic spiking activities (Fig. 3). The corresponding return map for the phase lags is shown in Fig. 8A Figure 8B illustrates the waveforms, as well as the bistability of the motif initially producing the counter-clockwise, 1[3[2 ð Þ , traveling wave that reverses into the clockwise one, 1[2[3 ð Þ , after a 10 second inhibitory pulse ended and released the blue reference cell to initiate a burst.
Medium duty cycle motif. To complete the examination of the influence of duty cycle on the repertoire and robustness of bursting outcomes of the homogeneous motif, we now consider the case of the medium-length duty cycle, 50%, set by V shift K2~{ 0:021 V (the middle interval shown in Fig. 3). Similarly to Figure 5A, Figure 9A illustrates the evolution of Dw 21 and Dw 31 (shown in blue and gray colors) from initial conditions uniformly distributed over the unit interval. One can  observe transients ultimately converging to multiple constant phase locked states. The corresponding map P is presented in Fig. 9B. In contrast to the case of short and long duty cycle motifs, the map for the medium duty cycle motif with weak homogeneous connections reveals the coexistence of five stable FPs: the red one at 0, 1 2 À Á , the green one at 1 2 ,0 À Á , the blue one at Þ . By externally applying a current pulse to a targeted cell we can deliberately switch between the co-existing bursting patterns (Fig. 10).

Asymmetric inhibitory motifs
In this section, we elucidate how and what intrinsic properties of the individual bursting cells affect the multistability of the 3-cell inhibitory motif. The answer involves an interplay between the competitive dynamical properties of individual neurons and the cooperative properties of the network. More specifically, it relies on how close an isolated cell is to the boundary between bursting and hyperpolarized quiescence and how sensitive the post-synaptic cell is to the (even weakly) inhibitory current generated by the pre- synaptic cells. We investigate these ideas by introducing asymmetries into the coupling of our homogeneous network motif. We focus on several representative cases of asymmetrically coupled motifs with one or more altered synaptic strengths, and we will elaborate on their bifurcations as we vary the asymmetry.
From multistability to the 1[3 [2 ð Þ pattern. In this subsection, we analyze bifurcations occurring en route from the homogeneous 3-cell motif to a rotationally-symmetric one, during which all clockwise-and counter-clockwise-directed synapses are simultaneously increased and decreased, respectively. In the limiting case of a clockwise, uni-directionally coupled motif there is a single traveling wave. The question is: in what direction will the wave travel?
We use a new bifurcation parameter, , which controls the rotational symmetry as the deviation from the nominal coupling strengths, g syn~5 |10 {4 , such that g syn (1+ ) and 0ƒ ƒ1.

The limit
?1 corresponds to the unidirectional case. Recall that initially, at~0, both the traveling waves Þare unstable in the short duty cycle motif with 20% DC. Then, the network can only generate the 1\f2E3g ð Þ , 2\f1E3g ð Þ , 3\f1E2g ð Þpacemaker rhythms. Figure 11A depicts P at a critical value of~0:41, and reveals that the FP (Dw 21 ,Dw 31 )~2 3 , 1 3 À Á is stable. Thus the counter-clockwise traveling wave, 1[3[2 ð Þ , is now observable in the asymmetric motif. The value~0:41 is a bifurcation value because further increase make the three saddles and the three initially stable FPs (blue, green and red), merge in pairs and annihilate though three simultaneous saddle-node bifurcations. After that, the FP around 2 3 , 1 3 À Á becomes the global attractor of the network (see Fig. 11B Fig. 12, focus on the area near this point, and shed light onto the intermediates in the bifurcation sequence. Figure 12A depicts an enlargement of P at~0:185. It shows a stable invariant curve near a heteroclinic connection involving all three saddles around the FP 2 3 , 1 3 À Á . In this figure, the FP near the center of the plot is still unstable. This indicates that the invariant curve has emerged from the heteroclinic connection at a smaller value of the parameter . The stable invariant curve is associated with the appearance of slow phase ''jitters'' demonstrated by the 1[3 [2 ð Þrhythm in voltage traces. As is increased further, the stable invariant curve shrinks down and collapses into the unstable FP 2 3 , 1 3 À Á making it stable through a secondary supercritical Andronov-Hopf (otherwise known as a torus bifurcation) as shown in Fig. 12B. Figure 16. Transformation stages of the phase lag maps for a motif with uni-directional asymmetry. Two connections g syn 31 and g syn 12 are strengthened from 1:03g syn through 1:5g syn . Due to the uni-directional symmetry breaking, the map first loses the clockwise, 1 3 , 2 3 À Á , FP (light gray) after it merges with a saddle at 1:05g syn , then the blue 1 2 , 1 2 À Á and the red 1 2 ,0 À Á FPs disappear through saddle-node bifurcations at 1:35g syn and 1:45g syn , respectively. As the counter-clockwise connections remain the same, the presence of the remaining FPs at 2  Bifurcations in the motif with one asymmetric connection. The homogeneous 3-cell motif has six independent connections, due to permutation properties we can limit our consideration of asymmetrically coupled motifs only to a few principle cases without loss of generality. First under consideration is the motif with a single synaptic connection, g syn 31 , from cell 3 to cell 1, being made stronger.
We first consider a perturbation to the homogeneous motif comprised of long duty cycle cells where just a single unidirectional connection, for instance from cell 2 to 3, is strengthened. To do this, we increase the coupling stenght g syn 31 from the nominal value, 5|10 {4 , through 1:5g syn , to 2g syn . This is effectively equivalent to increasing the parameter V shift K2 only for cell 3, thus pushing it toward the quiescence boundary and extending its interburst intervals. The corresponding maps are shown in Fig. 13. We observe that the initial increase of g syn 31 breaks the clockwise symmetry of the motif and makes the stable node at 2 3 , 1 3 À Á and a saddle come together. This motion further shrinks the attractor basin of the 1[3 [2 ð Þpattern. When g syn 31 is increased to 2g syn , both FPs have annihilated through a saddlenode bifurcation. In the aftermath, the unperturbed FP at 1 3 , 2 3 À Á remains the unique attractor of such an map. In turn, the asymmetric motif can stably produce the single bursting pattern, which is the 1[2[3 ð Þtraveling wave.
As our case study throughout the rest of the paper, we use the non-homogeneous 3-cell motifs composed of bursting cells with 50% duty cycle at V shift K2~{ 0:021 V. Figure 14 depicts the stages of transformation of the phase lag maps for the motif with the connection g syn 31 increasing from 1:04g syn and 1:4g syn through 1:6g syn . Inset A of Fig. 14  À Á annihilate through a saddle-node bifurcation. This widens the attractor basin (colored red in the figure) of the most robust FP at 0, 1 2 À Á after the clockwise traveling wave has been eliminated at g syn 31~1 :4g syn , as shown in Fig. 14B. At this value of g syn 31 , the 3\f1E2g ð Þ rhythm dominates over the remaining bursting rhythms because the red cell 3 produces more inhibition than the other two. To justify this assertion we point out that another motif, with weakened clockwise connections (g syn 12~g syn 23~0 :9g syn ) generates the identical Poincaré return map to the one shown in Fig. 14A.
In the 3\f1E2g ð Þrhythm, cell 3 bursts in anti-phase with the synchronous cells 1 and 2 that receive evenly balanced influx of inhibition from cell 3. This is no longer the case after the  Key Bifurcations of Bursting Polyrhythms in CPGs connection g syn 31 is made even stronger, so that the active cell 3 cannot hold both quiescent the postsynaptic cells 1 and 2 due to uneven coupling weights g syn 31~1 :6g syn 32 in the motif. One can see from the corresponding map Fig. 14B that the red FP at 0, 1 2 À Á is approached by a saddle point from the left at g syn 31~1 :4g syn . The map in Fig. 14C reveals that increasing g syn 31 through 1:6g syn causes a drastic change in the motif: the dominant red FP has vanished through a subsequent saddle-node bifurcation and so has the 3\f1E2g ð Þrhythm. With a single asymmetric connection, the structure of the phase lag map remains intact. However, the figure shows that the counter-clockwise wave has become the most robust rhythm, as the corresponding FP at 2 3 , 2 3 À Á has the largest attractor basin in the initial phase distribution.
Pyloric circuit motif. As an example, we examine bifurcation scenarios that occur as we transition to a heterogeneous motif that resembles the crustacean pyloric circuit with one inhibitory connection missing [1,14,22,34]. Such a network can be also treated as a sub-motif of a larger crustacean stomatogastric network [1].
The transformation stages are singled out in Fig. 15, which shows the bifurcations of the FPs in the phase lag maps. As in the previous case, decreasing a single either clockwise or counter-clockwise directional connection removes the corresponding FP at À Á , respectively. In this given case, it is the stable clockwise 1 3 , 2 3 À Á FP that vanishes though a saddle-node bifurcation after g syn 23 is decreased below 0:9g syn . Meanwhile, for g syn 23 v0:86g syn , cell 2 cannot maintain the synchrony between cells 1 and 3 in the 2\f1E3g ð Þrhythms, which is explained by a similar argument. This assertion is supported by the phase lag maps in Fig. 15B-C: one of the saddles shifts toward to the green FP at 1 2 ,0 À Á and annihilates it though a subsequent saddle-node bifurcation as g syn 23 is decreased through 0:85g syn . The principal distinction from the prior case is that one connection, g syn 31 , is made twice as strong as the others in the prior case, while here we completely remove a single connection in the limit g syn 23~0 . A consequence is that the basin of the stable FP at 2 3 , 1 3 À Á breaks down after it vanished through the third saddle node bifurcation that occur with the single connection been taken out, even while the three counter-clockwise connections remain intact. Its ''ghost'' remains influential, however, for some initial phase lags the motif can generate a long transient episode resembling the 1[2[3 ð Þ traveling wave. This wave eventually transitions into the dominant anti-phase 1\f2E3g ð Þrhythm that coexists with the less robust 3\f1E2g ð Þrhythm. In the phase plane, the ''ghost'' is located in a narrow region of transition between two saddle thresholds separating the attractor basins, blue and red, of the remaining stable FPs at 1 2 , 1 2 À Á and 0, 1 2 À Á . Finally, removing the g syn 23connection leaves the red attractor at 0, 1 2 À Á and its basin intact in Fig. 15D.

Two
asymmetric connections: uni-directional case. Here, we examine the motif with two uni-directional connection asymmetries, for example where g syn 12 and g syn 31 are strengthened from the nominal value to 1:5g syn . The bifurcation stages of P are depicted in Fig. 16 During the transformations, the map loses three FPs in sequence through similar saddle-node bifurcations. Because increasing g syn 31 and g syn 12 breaks the clockwise symmetry, the corresponding FP at 1 3 , 2 3 À Á for the counterclockwise wave, 1[2[3 ð Þ , is annihilated first at around 1:05g syn after merging with a saddle. Further strengthening both corrections annihilates the blue FP at 1 2 , 1 2 À Á , followed by the red FP at 1 2 ,0 À Á . As such, the pacemaker 1\f2E3g ð Þand 3\f1E2g ð Þ rhythms eventually are no longer available as neither cells 1 nor 3 are able to hold the post-synaptic counterparts in synchrony, and also because the periods of the unevenly driven cells become too different.
The clockwise symmetry breaking does not affect counterclockwise connections. Thus, in the map for 1:5g syn , two rhythmic patterns persist: the 1[3[2 ð Þ traveling wave with a wide attractor basin and the pacemaker 2\f1E3g ð Þ rhythm. Their associated FPs are at 2 3 , 1 3 À Á and 1 2 ,0 À Á , respectively. It is worth noticing that the same sequence of bifurcations will not occur in the map and the motif if only the connection g syn 23 is weakened instead.
Two asymmetric connections: Unilateral dominance case. Next under consideration is a motif in which cell 1 alone produces stronger inhibitory output due to two strengthened connections, g syn 12 and g syn 13 . Figure 17 depicts two snapshots of the phase spaces of the map after g syn 13 and then g syn 12 have been strengthened. One sees that a 10% increase in inhibition in the counter-clockwise direction breaks the rotational symmetry and therefore makes the stable FP at 2 3 , 1 Þrhythm) disappear through a saddle-node bifurcation as it merges with a saddle. As in the previous cases, the attractor basin of the stable blue pacemaker at 1 2 , 1 2 À Á extends to absorb that of the former FP. As expected, since all counter-clockwise connections have remained equal in this case, the stable FP at  Next, in addition to g syn 13 , the second outgoing connection, g syn 12 , from cell 1 is strengthened thus breaking the clockwise symmetry as well. As expected, this eliminates the FP at 1 3 , 2 3 À Á and the corresponding clockwise 1[2[3 ð Þ traveling pattern from the motif. Figure 17B shows the map for the motif with g syn 12~g syn 13~1 :5g syn . While it retains all three ''pacemaker'' FPs, the one at 1 2 , 1 2 À Á corresponding to the strongly inhibiting presynaptic cell 1 possess the largest attractor basin.
We may conclude that strengthening a single directional connection, or alternatively, a simultaneous and proportional weakening coupling strengths of the two synaptic connections of the same orientation in the motif, controls one of the three saddle points between the FP corresponding to traveling waves and the pacemaker FP corresponding to the stronger inhibiting cell. This will eventually causes the disappearance of either point as soon as the rotational symmetry is broken after the coupling strength is increased over some critical value, which varies depending on the nominal value g syn and the duty cycle of the bursting cells.
Motifs with a stronger coupled HCO: loss of phaselocking. A 3-cell motif with the cells coupled reciprocally by inhibitory synapses can be viewed alternatively as a group of three half-center oscillators (HCO). Each HCO represents a pair of cells that typically burst in anti-phase, when isolated from other cells. When a HCO is symmetrically driven, even weakly, by another bursting cell, it can produce in-phase bursting, instead of out of phase bursting [16].
In this section, we consider transformations of rhythmic outcomes in the motif containing a single HCO with stronger reciprocally inhibitory connections, for example, g syn 12~g syn 21~1 :25g syn (see Fig. 18A). It turns out that a 25% increase in coupling is sufficient to break both rotational symmetries because it eliminates the associated FPs around . In other words, the motif can still produce the coexisting 3\f1E2g ð Þrhythm. The following bifurcation sequence involving the dominant FP differs drastically from the saddle-node bifurcations discussed earlier. Observe from the map in Fig. 18A that two saddles separating two attractor basins, have moved close to the blue and green FPs as the coupling between the HCO cells is increased to 1:5g syn . This is a direct indication that a further increase of the coupling strength between the strongly inhibitory cells 1 and 2 will cause two simultaneous saddle-node bifurcations that eliminate both stable FPs.
A feature of these bifurcations of the map at the critical moment is that there are two heteroclinic connections that bridge the saddlenode FPs on the 2D torus. The breakdown of the heteroclinic connections with the disappearance of both FPs results in the emergence of a stable invariant circle that wraps around the torus [41,42]. The attractor basin of the new invariant curve is bounded away from that of the red FP at 0, 1 2 À Á by the stable sets (i.e., incoming separatrices) of the two remaining saddles. This motif is therefore bi-stable as the corresponding map shows two co-existing attractors.
Further increase in the coupling strength between the stronger inhibitory HCO and cell 3 cannot not qualitatively change the structure of the phase lag map, while it can have only a qualitatively effect on the size of the attractor basins of the invariant circle and the remaining FP (red). So, weakening g syn 31~g syn 32~0 :8g syn makes the separating saddles come closer to the red FP and hence shrink its attractor basin, as seen in Fig. 18B. This is not the case when either connection between cell 3 and the HCO is made sufficiently asymmetric. Depending on the connection's direction of asymmetry, such an imbalance causes either of the two remaining saddles to come close and annihilate with the stable red FP at 1 2 ,0 À Á . Figure 19 presents the map for this motif with weakened reciprocal connections between cells 3 and 1: g syn 13~g syn 31~0 :8g syn . This motif, comprised of three HCOs with strong, nominal and weak reciprocal connections, no longer produces any phase-locked bursting rhythm, including 3\f1E2g ð Þ , as the map no longer has any stable FPs. The resulting motif is monostable with a single attractor for the stable invariant curve. This curve can be characterized with a rational or irrational winding number. The number is a rational if the invariant curve is made of a finite number of periodic points across the torus.
The occurrence of the stable invariant curve wrapping around the torus gives rise to a phase slipping phenomenon observed in voltage traces such as those shown in Fig. 19B. We define ''phase slipping'' as a repetitive rhythm with varying phase lags between the bursting cells of the motif. The period of the invariant circle depends on how far the map with the invariant circle is from the bifurcations of ''ghost'' FPs. The ''ghosts'' make the bursting pattern with varying phase lags appear as it is composed of four sequential episodes and transitions between them.
From the top of the Dw 21 ,Dw 31 ð Þ -unit square, the curve begins with the 2\f1E3g ð Þrhythm continuously transitioning into the clockwise 1[2 [3 ð Þtraveling wave, followed by the 1\f2E3g ð Þ rhythm, and being followed by the counter-clockwise 1[3[2 ð Þ traveling wave and finally returning to the initial 2\f1E3g ð Þ rhythm in nine bursting cycles, which is the period of the phase slipping. Each episode of the phase slipping rhythm can be arbitrarily large as it is controlled by the coupling strength of the specific motif connections near the corresponding saddle-node bifurcation(s). Observe that Dw 21^1 2 on the invariant curve, i.e., while cell 1 and 2 are in anti-phase bursting, cell 3 modulates the rhythm by recurrently slowing down and advancing the HCO to generate continuously all four episodes.
One may wonder about what determines the direction of the invariant curve on the torus and hence the order of the episodes of the shown voltage waveform. Observe that the phase slip occurs in the Dw 31 direction and that the invariant curve, unlike a FP, has no fixed period for the whole network. Indeed, the recurrent times of this network change periodically, approximately every eight episodes. The eight episodes constituting the bursting pattern are determined by a rational ratio of the longer HCO period (due to stronger reciprocal inhibition that extends the HCO interburst intervals) to the shorter period of pre-synaptic cell 3 (due to a weaker incoming inhibition) (see Figs. 19B and 20).
Let us discuss another motif configuration, shown in Fig. 18B, that produces maps with stable invariant curves that wrap around the torus. These are qualitatively identical to the maps for the motif containing the strong HCO formed by cells 1 and 2 (Fig. 18A). In this configuration, cell 3 receives weaker inhibition from pre-synaptic cells 1 and 2 according to g syn 13~g syn 23~0 :6g syn . The de-stabilizing factor 0.6 turns out to be small enough to make sure that neither cell 1 nor 2 can be a pacemaker as the corresponding stable FPs have disappeared because the period of cell 3 has become shorter than the periods of cells 1 and 2. As the result, the map demonstrates the same stable invariant curve that ''flows'' downwards with decreasing Dw 31 phase lags.
The direction of the stable invariant circle flowing across the 2D torus can be reversed by making cell 3 receive stronger inhibition instead of weaker inhibition relative to the other cells. An example is depicted in the phase lag map of Fig. 20, where g syn 13~g syn 23~1 :6g syn . Toward control of multistability. We now elucidate the issues involved in designing inhibitory motifs with predetermined bursting outcomes and how to control them. Let us revisit the motif with a single HCO in Fig. 17A. The map is depicted near the bifurcations that eliminate both blue and green FPs simultaneously as the coupling strength between cells 1 and 2 is increased. The corresponding saddle-node bifurcations are each of co-dimension one, i.e. can be unfolded by a single parameter. This means that increasing either coupling parameter, g syn 12 or g syn 21 , makes the respective FP at 1 2 ,0 À Á (green) or 1 2 , 1 2 À Á (blue) disappear or re-emerge. This suggests alternative ways of perturbing the motif to get the desired outcome. For instance, in the motif with the HCO given by g syn 12~g syn 21~1 :5g syn , cell 2 can be made the strongest on the motif by increasing the outgoing inhibitory drive: g syn 23~1 :5g syn . The green FP at 1 2 ,0 À Á in the corresponding map in Fig. 21A has a largest attraction basin that guarantees the dominance of the 2\f1E3g ð Þ -rhythm over the network. The map in Fig. 21B has the basin of the blue FP at 1 2 , 1 2 À Á largest, after strengthening the coupling from cell 1 to 3 in the motif with two robust bursting outcomes: the 1\f2E3g ð Þ -rhythm dominating over the 3\f1E2g ð Þ -rhythm corresponding to the red FP at 0, 1 2 À Á with a smaller basin formed by initial phases. The above configurations of the inhibitory motif are bistable with two coexisting FPs: dominant blue (or green) with a large attractor basin and red with a smaller basin corresponding to the less robust 3\f1E2g ð Þrhythm. To construct the monostable motif with the single rhythm, for example 1\f2E3g ð Þ , cell 1 must be coupled reciprocally stronger with cell 3 than cell 2. Such a motif has two HCOs that both contain cell 1 due to the strengthened pairs of synaptic connections: g syn 12~g syn 21~1 :5g syn and g syn 13~g syn 311 :5g syn . The corresponding map for the phase lags is shown in Fig. 22. The resulting map demonstrates that both the red and green FPs have been annihilated, as well as the corresponding bursting rhythms. Note that the map still has two saddle FPs in addition to the only attractor at 1 2 , 1 2 À Á . It is a feature of a map on a 2D torus that the number of FPs must be even, in general, for them to emerge and vanish through saddle-node bifurcations. Therefore, the map must possess another hyperbolic FP. This point resides near the origin where all three cells burst synchronously, which we consider next.
Fine structure near the origin. A common misconception concerning modeling studies of coupled cells is that fast, nondelayed inhibitory synapses always foster anti-phase dynamics over unstable in-phase bursting. While being true in general for simple relaxation oscillators, interactions of bursting cells can be incomparably more complex even in small networks including HCOs with fast inhibitory coupling [16,27]. It was shown in [28] that overlapped bursters can reciprocally synchronize each other in multiple, less robust, phase-locked states due to spike interactions. Furthermore, the number of such synchronous steady states is correlated with the number of spikes within the overlapped bursts.
To explore the dual role of inhibition, we now explore nearly synchronous bursting in all three cells of the homogeneous, medium DC motif. Because synchronous steady states are due to spike interactions, we restrict the consideration to a relatively small positive vicinity of the synchronous state, Dw 21~D w 31~0 , in P. A magnified portion of the map is shown in Fig. 23, where green, red and black dots denote the locations of the stable, repelling and saddle (threshold) phase locked states (respectively) for the nearly synchronized bursting outcomes. The map reveals that several overlapping burst patterns can occur where either cell spikes slightly in advance or delayed compared to the reference cell. Unstable FPs surround the outer part of this small region of the map make the origin repelling in the map on the global scale.

Excitatory motifs
In this section, we discuss a variation of a homogeneous 3-cell motif with short, 25% DC at V shift K2~{ 0:01895 V, will all three excitatory synaptic connections. The synaptic current is again given through the FTM paradigm: The synapses are made  more excitatory by increasing the synaptic reversal potential, E syn , from {0:0625 V (corresponding to the inhibitory case) to 0:0 V. E syn~0 guarantees that the voltages of all the cells remain below the reversal potential, on average, over the bursting period. In the excitatory motif, whenever the advanced cell initiatives a new bursting cycle, the synaptic current raises the voltages of postsynaptic cells, thus making it follow the pre-synaptic one, at the hyperpolarized knee-point on the quiescent manifold (Fig. 1B). Figure 6 shows the phase lag map for the original inhibitory motif with three stable FPs (shown in blue, red and green) at À Á . The attractor basins of three stable FPs are separated by the separatrices of six saddle FPs (smaller dots). A small area around the origin is globally repelling. This motif can stably produce three coexisting patterns in which either cell bursts in anti-phase with the two remaining in-phase.
It is often presumed in neuroscience that excitation acts symmetrically opposite to inhibition in most cases, i.e. wherever inhibition tends to break synchrony, excitation fosters it. Figure 24 supports this assertion for this particular kind of network and coupling. It depicts the map corresponding to the homogeneous 3cell motif with reciprocally excitatory connections for same short, 25% DC.
Compared to the map for the inhibitory motif, the map for the homogeneously excitatory motif is the inverse: 21 ,Dw (nz1) 31 ? Dw (n) 21 ,Dw (n) 31 ; here the inverse is the forward map in discrete backward time. As such, the FPs at 1 3 , 2 3 À Á and 2 3 , 1 3 À Á , which used to be repelling in the inhibitory case, become attracting but with smaller basins. This means that the motif can generate traveling waves, albeit with low probability. Meanwhile, the FPs colored blue, green and red, are now repellers, and hence none of the pacemaker rhythms can occur. Reversing the stability does not change the topological type of the six saddles, but their stable and unstable separatrices are reversed. The dominant attractor of the map is now located at the origin, to which nearly all transient trajectories converge. This implies that the reciprocally excitatory motif, whether homogeneous or heterogeneous, will exhibit stable synchronous bursting with all three cells oscillating in-phase.

Mixed motifs
Here we discuss two intermediate configurations of mixed motifs having both inhibitory and excitatory connections. First, we consider the motif with a single excitatory connection from cell 3 to 1. Its coupling strength is regulated by the level of the synaptic reversal potential, E 31 syn . Figure 25 depicts three phase lag maps for the motif with E syn being increased from {0:050, {0:030 through 0:0 V.
Initially, an increase in E 31 syn gives rise to two saddle-node bifurcations in the motif (Fig. 25A): the first one breaks the clockwise rotational symmetry and hence annihilates the stable FP at 1 3 , 2 3 À Á . The second bifurcation annihilates the stable red point at 0, 1 2 À Á , because cell 3, inhibiting 2 and exciting 1, cannot hold both of them at the hyperpolarized quiescent state to generate the 3\f1E2g ð Þ -rhythm as it promotes burst initiation in cell 1 following those in cell 3. On the contrary, excitation applied to cell 1 forces it to follow cell 3 after a short delay in the burst initiation. As the result, the disappearance of the 3\f1E2g ð Þrhythm promotes the 2\f1E3g ð Þ -rhythm and an increase of the attractor basin of the green FP.
Initial elevations of the level of E 31 syn keep the other three FPs intact, while widening the basins of the blue and green stable FPs at 1 2 , 1 2 À Á and 1 2 ,0 À Á . Further increasing E 31 syn increases the duty cycle of the blue cell by extending its active bursting phase. Consequently, the counter-clockwise ring no longer contains identical cells that could orchestrate the 1[3 [2 ð Þpattern. This patten is eliminated with the disappearance of the corresponding FP at 2 3 , 1 3 À Á through a merger with a saddle. The map now has two persistent attractors, blue and green, as shown in Fig. 25B. With E 31 syn increased still further, the blue cell 1 receives strongly unbalanced input: larger excitation influx from the postsynaptic cell 3 and an inhibitory drive from cell 2, acting oppositely. This unbalanced input increases the active phase of bursting of cell 1 and hence its duty cycle and period, and hence breaks cell 1's ability to robustly maintain the 1\f2E3g ð Þ -rhythm by evenly inhibiting the pots-synaptic cells 1 and 2 of the same period. In P, this results in the shrinking of the attractor basin of the blue FP, whereas the basin of the dominating green FP widens. By setting E 31 syn~0 :0 V, the resulting strong imbalance between excitation and inhibition onto cell 1 makes the 1\f2E3g ð Þ -rhythm impossible to occur in the network and the corresponding FP at 1 2 , 1 2 À Á disappears in the map. After this last saddle-node bifurcation, the map has a unique attractor in the green FP. Eventually, regardless of initial phases, the synergetic interaction of inhibitory and excitatory cells in this mixed motif will give rise to the 2\f1E3g ð Þ -rhythm led by cell 2 (Fig. 26). Finally, we consider the mixed motif with two excitatory connections originating from cell 3. Figure 27 depicts the transformations of the P as the reversal potentials, E 31 syn and E 32 syn , are increased simultaneously from {0:050, {0:030 through to {0:020 V. The increase makes postsynaptic cells 1 and 2 more excited compared to cell 3, which consequently receives a longer duration of inhibition.
As in the previous case, increasing E 31 syn~E 32 syn~{ 0:050 V breaks both rotational symmetries, which is accompanied by the disappearance of the corresponding, counter-clockwise and clockwise, FPs. The increased excitability of cells 1 and 2 initiates the active bursting states of those cells soon after that of cell 3 in each cycle. Thus, the 3\f1E2g ð Þ -rhythm led by cell 3 is less likely to occur compared to the 1\f2E3g ð Þand 2\f1E3g ð Þ -rhythms. The corresponding (red) FP at 0, 1 2 À Á loses the attractor basin in the map and then disappears, following the FPs for the traveling waves, at {0:030 V, see Fig. 27B. After that, the blue and green stable points have equal attractor basins corresponding to equal chances of the emergence of the phase-locked rhythms 1\f2E3g ð Þ and 2\f1E3g ð Þ . Examination of the map suggests that besides these phase-locked rhythms, the motif can generate long transients with episodes that resemble those of the 3\f1E2g ð Þ -rhythm transitioning back and forth with in-phase bursting. Such transients are due to regions of the map that are constrained by the incoming separatrices of the remaining saddles, which are forced to curve in complex ways to embed onto the torus with two attractors and an unstable origin. Note that the origin may not longer be a repeller as a whole but has a saddle structure, because it is no longer associated with synchronous bursting.
Setting the excitatory reversal potential to zero annihilates the two remaining phase-locked states (blue and green FPs) in two simultaneous bifurcations. As with the case of the inhibitory motif with the HCO (compare with Fig. 18), these global bifurcations underlie the formation of a closed heteroclinic connection between the saddle-node points at the critical moment. These connections transform into an invariant circle that wraps around the torus. Having settled onto the invariant curve that zigzags over the torus, the network will generate long recurrent patterns consisting of three transient episodes: namely, in-phase bursting that transitions to the 1\f2E3g ð Þ -rhythm, which transitions back to in-phase bursting, then transitioning to the 2\f1E3g ð Þ -rhythm, which then returns to the in-phase bursting, and so forth.
At higher values of reversal potentials, excitation overpowers inhibition and this mixed motif fully becomes the excitatory motif with the single, all-synchronous, bursting rhythm forced by the driving cell 3. In the corresponding return map, this rhythm occurs after the invariant curve terminates at a homoclinic connection to a saddle-node near the origin so that the origin becomes a global attractor again.

Gap junction in an inhibitory motif
Finally, let us consider the role of a single electrical synapse though a gap junction between cells 1 and 2 in the inhibitory motif (Fig. 1A). The difference in the membrane potentials gives rise to a directional ohmic current described by I el~gel (V 2 {V 1 ) in the model (3). Figure 28 depicts the stages of transformation of the corresponding maps as g el is increased from 10 {4 through 3|10 {4 .
The electrical coupling breaks down the rotational symmetry that causes the disappearance of the corresponding FPs at 1 3 , 2 3 À Á and 2 3 , 1 3 À Á . The disappearance of both FPs widens the attraction basin of the red FP at 0, 1 2 À Á compared to the individual basins of the blue and green FPs. The bidirectional electrical coupling tends to equilibriate the membrane potentials of the connected cells, so that cell 1 and 2 are brought closer together to burst in synchrony, rather than in alternation. At intermediate values of g el , inhibitory coupling can still withstand the tendency to synchronize cells 1 and 2, while the red basin widens further due to shrinking basins of the blue and green FP. At larger values of g el , synchrony of bursting cells 1 and 2 overpowers their reciprocal inhibition, and the motif generates only the 3\f1E2g ð Þ -rhythm regardless of the choice of initial phases. Thus, we find that motifs with strong electrical synapses describe a dedicated rather than a multifunctional CPG.

Discussion
Our new computational technique reduces the dynamics of a 9dimensional network motif of three cells to the analysis of the 2D maps for the phase lags between the bursting cells. With this technique, we demonstrated that a reciprocally inhibitory network can be multistable, i.e. can generate several distinct polyrhythmic bursting patterns. We studied both homogeneous and nonhomogeneous coupling scenarios as well as mixtures of inhibition, excitation and electrical coupling. We showed that the observable rhythms of the 3-cell motif are determined not only by symmetry considerations but also by the duty cycle, which serves as an order parameter for bursting networks. The knowledge of the existence, stability and possible bifurcations of polyrhythms in this 9D motif composed of the interneuron models is vital for derivations of reduced, phenomenologically accurate phase models for nonhomogeneous biological CPGs with inhibitory synapses.
The idea underlying our computational tool is inspired by common features of electrophysiological experimentation. As such, it requires only the voltage recording from the model cells and does not explicitly rely on the gating variables. We intentionally choose the phases based on voltage as often this the only experimentally measurable variable. Moreover, as with real experiments, we can control the initial phase distribution by releasing the interneurons from inhibition at specific times relative to the reference neuron. Our analysis of the system only utilizes the qualitative, geometric tools of dynamical systems theory. Thus, although in this example we simulated a system of underlying differential equations to generate the maps, in principle such maps could be generated and analyzed directly from experimental data. In this sense, our method does not require knowledge of explicit model equations.

Summary and interpretation of main results
Rhythmic motor behaviors, such as heartbeat, respiration, chewing, and locomotion are often independently produced by small networks of neurons called Central Pattern Generators (CPGs). It is not clear either what mechanisms a single motor system can use to generate multiple rhythms, for instance: whether CPGs use dedicated circuitry for each function or whether the Here, E 31 syn~{ 0:020 V, which corresponds to the phase lag map having a single attractor at the green FP in Fig. 25C; here stronger g syn~0 :005 was used for the sake of illustration of a quicker convergence. doi:10.1371/journal.pone.0092918.g026 same circuitry can govern several behaviors. A systematic way to explore this is to create mathematical models that use biologically plausible components and classify the possible varieties of rhythmic outcomes.
We performed such a study of CPG networks based on three inter-connected neurons. We systematically varied the strength and functional form of coupling between the neurons to discover how these affect the behavioral repertoire of the CPG. To do this, we created a geometric representation of the simulated CPG behavior of each possible configuration of the network, which greatly simplifies the study of the 9-dimensional system of nonlinear differential equations. We discovered several configura-tions that support multiple rhythms and characterized their robustness. By varying physiologically reasonable parameters in the model, we also describe mechanisms by which a biological system could be switched between its multiple stable rhythmic states.
In the Homogeneous Inhibitory Motifs section, we showed that a weakly coupled, homogeneous motif comprised of three bursting interneurons with reciprocally inhibitory fast synapses can produce a variable number of polyrhythms, depending on the duty cycle of the individual components. The phase lag maps are de facto proof of the robust occurrence of the corresponding rhythmic outcomes generated by such a motif. While the occurrence of some rhythms, Þ , in a 3-cell motif can hypothetically be deduced using symmetry arguments; the existence and robustness of the rhythms can be only verified by accurate computations of the corresponding return maps. Moreover, the observability of these rhythms in even the homogeneous motifs, and the stability of the FPs, are both closely linked to the temporal properties of the bursts.
Recall that the inhibitory current shifts the post-synaptic cell closer to the boundary or can even move it over the boundary into silence while the pre-synaptic cell remains actively spiking (Fig. 3). In terms of dynamical systems theory, this means that the perturbed driven system closes the gap between the hyperpolarized fold and the slow nullcline m 0 K2~0 , eventually causing the emergence of a stable equilibrium state on the quiescent branch. As such, the homogeneous network produces three pacemaker rhythms, 3\f1E2g ð Þ , 2\f1E3g ð Þ , and 1\f2E3g ð Þ -the only rhythms available in the short motif. These strongly synchronized activities imply fast convergence to the phase locked states because of the emergent equilibrium state near the closed gap: compare the time spans in Figs. 6A and 9A.
The gap never gets closed by weak inhibition in the long duty cycle motif as the individual cells have initially remained far enough from the boundary separating the bursting and hyperpolarized quiescence. As the result, this motif can only effectively produce two possible bursting outcomes: the clockwise 1[2[3 ð Þ and counter-clockwise 1[3 [2 ð Þtraveling waves. This bistability results from a weaker form of synchronization, which is confirmed by the rate of convergence to the FPs. This will not be the case when the inhibition becomes stronger due to increasing the nominal coupling strength, g syn .
The intermediate case of the medium motif is far from the above extremes and benefits from a natural optimization between the coupling strength, initial phase distributions, and the spatiotemporal characteristics of unperturbed and perturbed bursting. One such characteristic is the slow passage through the ghost of the stable equilibrium state that guarantees the robust synchrony in the short duty cycle motif. As the result, the motif can produce five stable bursting rhythms: the anti-phase 3\f1E2g ð Þ , 2\f1E3g ð Þ , and 1\f2E3g ð Þ ; and the clockwise 1 Þtraveling waves. In the Asymmetric Inhibitory Motifs section, we described a number of generic bifurcations of the original five polyrhythms that can occur in the homogeneous, reciprocally inhibitory motif with the medium duty cycle. We revealed the basic principles of transformations of such a multi-functional network into ones with fewer rhythms or even with a single pattern bursting pattern.
The rhythmic outcomes of the CPG do not always have to involve phase locking, as there can be a stable pattern of phase slipping bursts that have time varying phase lags.
While each rhythm remains robust with respect to variations of the coupling connections, one can still make the network switch between them by applying an external pulse to the targeted cell that, upon release, appropriately changes the relative phases of the cells, as demonstrated in Fig. 10 for the medium duty cycle motif. In terms of the Poincaré return maps for the phase lag, switching between rhythms is interpreted as switching between the corresponding attractor basins of FPs or invariant circles. This causes the state of the network to ''jump'' over the separating threshold defined by a saddle (more precisely, over the incoming separatrices of the saddle). We stress, however, that the choice of timing in the suppression of a targeted cell to effectively switch between these polyrhythms is not intuitive and requires a detailed understanding of the underlying dynamics.
Although there are alternative ways of creating the 3-cell reciprocally inhibitory motif with predetermined outcomes, the fundamental principles are universal: the pre-synaptic cell that produces stronger inhibition gains the larger attractor basin and therefore the corresponding rhythm led by this cell becomes predominant. In particular, a sufficient increase (or decrease) of the coupling strength of a single connection can break the intrinsic rotational symmetry of the motif to remove the possibility of observing traveling wave patterns.
Of special interest are the motifs with asynchronous phaseslipping patterns that have no dominating phase-locked states. Such patterns result from the synergetic interactions between all contributing cells, and is comprised of four transient episodes, but primarily marked by a continuous transition between two primary sub-rhythms: 2\f1E3g ð Þ and 1\f2E3g ð Þ (for example, see Fig. 20). Both competing sub-rhythms are equally possible, and none can prevail over the other without the weaker inhibiting cell 3 whose reciprocal connections change the existing balance by shifting the phase lags during all four episodes.
In all cases we have considered, inhibition was chosen weak enough to guarantee that the post-synaptic cell remains in a bursting state even when its duty cycle decreases to 20% near the boundary between bursting and silence. This ensures that the phases of the cells, as well as the phase lags, are well defined for the return maps. However, this assumption may fail, for example when phase of a post-synaptic cell is not defined because the incoming synaptic inhibition makes it quiescent [17]. This leads to a phenomenon called sudden death of bursting that occurs when a rhythmic leader (red cell 3 in Fig. 29) becomes suppressed by the  other two cells that form an anti-phase HCO, which alternately inhibit the post-synaptic cell. In our example, the outgoing inhibition from cell 3 is several times weaker than the inhibition that it receives from the HCO formed by cells 1 and 2: g syn 31~g syn 32~0 :4g syn and g syn 13~g syn 23~6 g syn (this ratio does not always have to be as large when g syn is increased).
It should be stressed that the sudden death of bursting co-exists with other bursting patterns in the motif at the same coupling parameters. As such, this example bears a close qualitative resemblance to the experimental recordings from identified interneurons comprising the leech heart CPG in which the socalled switch interneuron alternately leads synchronous patterns and then becomes inactive during peristaltic patterns [43][44][45]. Analogous reversions of direction in the blood circulation, peristaltic and synchronous, are observed in the leech heartbeat CPG and its models [21,46]. The contrast between the patterns is that switching appears to be autonomously periodic in the leech, i.e. occurs without external stimulus, in a way similar to the phaseslipping pattern presented in Fig. 19.
Examination of the complex fine structure of the map near the origin reconfirms that fast, non-delayed inhibition can have stabilizing effects leading to the onset of several nearly synchronous bursting patterns with several overlapping spikes [28]. Such bursting patterns are less robust compared to those corresponding to FPs with large attractor basins in the phase lag maps.
We showed in the Excitatory Motifs section that raising the synaptic reversal potential is equivalent to reversing the time in the inhibitory motif model. In the maps, this action transforms attractors into repellers, while saddles remain saddles but have their incoming and outgoing directions reversed. To fully illustrate this phenomenon, we used the symmetric motif with the short duty cycle, in which the clockwise and counter-clockwise traveling waves are unstable. In the symmetric counterpart with excitatory synapses, these FPs become the attractors along with the dominant fixed point at the origin corresponding to synchronous bursting.
In the Mixed Motifs section we analyzed the transformation of the maps corresponding to the motifs with mixed, inhibitory and excitatory synapses. We showed step-by-step how the the multifunctional motif becomes a dedicated motif. The final example was the mixed motif with two excitatory connections. Unlike the former case, such a monostable motif possesses a bursting rhythm with time varying phase lags, which corresponds to a stable invariant curve wrapping around the 2D torus.

Conclusions and future work
We emphasize that a highly detailed examination of the occurrence and robustness of bursting patterns in the 3-cell motifs would be impossible without the reduction of the complex 9D network model with six algebraic equations for chemical synapses to the 2D maps and the numerical bifurcation analysis of FPs, invariant circles, homoclinic structures in it. Recall that the dimension of the map is determined by the number of the nodes in the network, not the number of differential equations per synapse and per neuron, which can be much greater. With the aid of our computational tools we were able to identify even some exotic bifurcations from the dynamical system theory like the Cherry flow (Fig. 27B) on the torus [41] occurring in this neural network. High accuracy of numerical simulations is required in our analysis. This involves at least a 40|40 mesh of initial conditions run for 100 bursting cycles to generate a single map and identify its structure. This computation takes up to three hours on a multi-core CPU workstation, but future work will take advantage of parallel computing architectures.
A stable FP of the map corresponding to a robust bursting pattern with specific phase lags is also structurally stable, i.e. persists under particular variations of coupling parameters. So by varying the given parameter(s) we can evaluate the boundaries and region of its existence and stability, which we will call a ''synchronization zone.'' A boundary of the region corresponds to a bifurcation, which can be either of saddle-node type, in which the FP vanishes, or of Andronov-Hopf type, through which the FP merely changes stability. Figure 30 sketches a possible arrangement of such zones in a parameter space of the network. They are shown to be nested within each other, because changing monotonically a single parameter can cause a cascade of saddle-node bifurcations in the map, such as those shown in Figs. 18 and 19. Given the number of the reciprocal synapses in the motif, the parameter space of the network is at least six dimensional, which presents a challenge for examining all bifurcations in detail and creating a complete catalogue. However, we have demonstrated that several inhibitory configurations of the 3-cell motif generate phase lag maps of qualitatively the same structure, see exemplary Fig. 1. This observation provides underlying foundations for highly effective reduction tools for studies of multi-component neural networks. It implies that variations of different coupling parameters make the network undergo same bifurcations while transversally crossing the corresponding bifurcation boundaries in the parameter space. This approach provides de facto proof that, without the phase lag maps, it would be impossible to claim and understand why two distinct network configurations produce same rhythmic outcomes.
In general, our insights allow us to predict both quantitative and qualitative transformations of the observed patterns as the network configurations are altered or the network states are perturbed dynamically [47]. The nature of these transformations provides a set of novel hypotheses for biophysical mechanisms about the control and modulation of rhythmic activity. A powerful aspect to our geometric form of analysis is that it does not require knowledge of the equations that model the system, for instance if the maps were generated from an unknown model system (or even from experimental data). For the sake of demonstration, we generated our maps from explicit differential equations. Our computational tools help us explain the fundamental dynamical mechanisms underlying the rhythmogenesis in plausible models of CPG networks derived from neurophysiological experiments [48]. Thus, we believe that have developed a universal approach to studying both detailed and phenomenological models that is also applicable to a variety of rhythmic biological phenomena beyond motor control.

Models and Numerical Methods
We study CPG network motifs comprised of three cells coupled reciprocally by non-delayed, fast chemical synapses and with weak coupling strengths. The time evolution of the membrane potential, V , of each neuron is modeled using the framework of the Hodgkin-Huxley formalism, based on a reduction of a leech heart interneuron model, see [49] and the references therein: The dynamics of the above model involve a fast sodium current, I Na with the activation described by the voltage dependent gating variables, m Na and h Na , a slow potassium current I K2 with the inactivation from m K2 , and an ohmic leak current, I leak : I Na~ g g Na m 3 Na h Na (V {E Na ), I K2~ g g K2 m 2 K2 (V {E K ), I L~ g g L (V {E L ): ð4Þ C~0:5 nF is the membrane capacitance and I app~0 :006 nA is an applied current. The values of maximal conductances are g g K2~3 0 nS, g g Na~1 60 nS and g L~8 nS. The reversal potentials are E Na~4 5 mV, E K~{ 70 mV and E L~{ 46 mV.
The time constants of gating variables are t K2~0 :9 s and t Na~0 :0405 s. The Fast, non-delayed synaptic currents in this study are modeled using the fast threshold modulation (FTM) paradigm as follows [50]: here V post and V pre are voltages of the post-and the pre-synaptic cells; the synaptic threshold H syn~{ 0:03 V is chosen so that every spike within a burst in the pre-synaptic cell crosses H syn , see Fig. 1. This implies that the synaptic current, I syn , is initiated as soon as V pre exceeds the synaptic threshold. The type, inhibitory or excitatory, of the FTM synapse is determined by the level of the reversal potential, E syn , in the post-synaptic cell. In the inhibitory case, it is set as E syn~{ 0:0625 V so that V post (t)wE syn . In the excitatory case the level of E syn is raised to zero to guarantee that the average of V post (t) over the burst period remains below the reversal potential. We point out that alternative synapse models, such as the alpha and other detailed dynamical representation, do not essentially change the dynamical interactions between these cells [28]. The numerical simulations and phase analysis were accomplished utilizing the freely available software PyDSTool (version 0.88) [51,52]. Additional files and instructions are available upon request.

Supporting Information
Movie S1 Polyrhythmic dynamics of 3-cell CPGs: from 3\f1E2g ð Þto 2\f1E3g ð Þpattern. Bursting trajectory (gray) in the 3D phase space of the model, which is made of the ''active'' spiking (solenoid-like shaped) and the flat hyperpolarized sections. The gap between the 2D slow nullcline, m 0 K2~0 , and the low knee on the slow quiescent manifold, M eq , determines the amount of inhibition needed by the active pre-synaptic cell above the synaptic threshold, H syn , to either slow or hold the post-synaptic cell(s) at a hyperpolarized level around {0:06 V. The red, green and blue spheres on the bursting trajectory depict the temporal evolution of the the phases of the [not weekly]-coupled cells of the CPG: active cell(s) above H syn inhibits, in anti-phase, the temporarily inactive cells and visa versa. Inhibitory pulse applied to the blue cell changes the relative phases of the bursting cells so that the CPG pacemaker becomes the green cell after the red cell. Below are shown the corresponding voltage waveforms. This motion picture complements Fig. 10.

(MOV)
Movie S2 Polyrhythmic dynamics of 3-cell CPGs: from 3\f1E2g ð Þto 1[2 [3 ð Þpattern. Bursting trajectory (gray) in the 3D phase space of the model, which is made of the ''active'' spiking (solenoid-like shaped) and the flat hyperpolarized sections. The gap between the 2D slow nullcline, m 0 K2~0 , and the low knee on the slow quiescent manifold, M eq , determines the amount of inhibition needed by the active pre-synaptic cell above the synaptic threshold, H syn , to either slow or hold the post-synaptic cell(s) at a hyperpolarized level around {0:06 V. The red, green and blue spheres on the bursting trajectory depict the temporal evolution of the the phases of the [not weekly]-coupled cells of the CPG: active cell(s) above H syn inhibits, in anti-phase, the temporarily inactive cells and visa versa. Inhibitory pulse applied to the blue cell changes the relative phases of the bursting cells so that the the pacemaking rhythm led by the red cell is replaced by the clockwise traveling wave (peristaltic bursting). Shown below are the corresponding voltage waveforms. This motion picture complements Fig. 10.

(MOV)
Movie S3 Multistable dynamics of an asymmetric 3-cell CPG: from 3\f1E2g ð Þpattern to sudden death. Bursting trajectory (gray) in the 3D phase space of the model, which is made of the ''active'' spiking (solenoid-like shaped) and the flat hyperpolarized sections. The gap between the 2D slow nullcline, m 0 K2~0 , and the low knee on the slow quiescent manifold, M eq , determines the amount of inhibition needed by the active presynaptic cell above the synaptic threshold, H syn , to either slow or hold the post-synaptic cell(s) at a hyperpolarized level around {0:06 V. The red, green and blue spheres on the bursting trajectory depict the temporal evolution of the the phases of the [not weekly]-coupled cells of the CPG: active cell(s) above H syn inhibits, in anti-phase, the temporarily inactive cells and visa versa. Inhibitory pulse applied to the blue cell changes the relative phases of the bursting cells so that the pacemaking red cell, leading initially the rhythm, becomes permanently inhibited by the blue and green cells bursting in alternation. Below are shown the corresponding voltage waveforms. This motion picture complements