2θ-Burster for Rhythm-Generating Circuits

We propose a minimalistic model called the 2 θ -burster due to two slow phase characteristics of endogenous bursters, which when coupled in 3-cell neural circuits generate a multiplicity of stable rhythmic outcomes. This model offers the benefits of simplicity for designing larger neural networks along with an acute reduction in the computation cost. We developed a dynamical system framework for explaining the existence and robustness of phase-locked states in activity patterns produced by small rhythmic neural circuits. Several 3-cell configurations, from multifunctional to monostable, are considered to demonstrate the versatility of the proposed approach, allowing the network dynamics to be reduced to the examination of 2D Poincaré return maps for the phase lags between three constituent 2 θ -bursters.


INTRODUCTION
Neural networks called central pattern generators (CPGs) [1][2][3][4][5][6][7][8] produce and control a great variety of rhythmic motor behaviors, including heartbeat, respiration, chewing, and locomotion. Many physiologically diverse CPGs involve 3-cell motifs, such as the spiny lobster pyloric network [6], the Tritonia swim circuit [4], and the Lymnaea respiratory CPGs [3]. Pairing experimental studies and modeling studies has been proven to be the key for disclosing basic operational and dynamical principles of CPGs [9][10][11][12][13][14]. Although various circuits and models of specific CPGs have been developed, the mystery of how CPGs gain the level of robustness and adaptation observed in nature remains unsolved. It is also not evident what mechanisms a single motor system can use to generate multiple rhythms, that is, whether CPGs need a specific circuitry for every function or whether it can be multifunctional to determine several behaviors [15][16][17]. Switching between multistable rhythms can be attributed to input-dependent switching between attractors of the CPG, where each attractor is associated with a specific rhythm. The goal of this article is to characterize how observed multistable states can emerge from the coupling using simple neural models on small networks. This article, based on our original work, reemphasizes some basic principles well established in the characterization of 3-cell networks made of detailed Hodgkin-Huxley-type models of endogenous bursters [18][19][20] and the Fitzhugh-Nagumo-like neurons [21]. We use a bottom-up approach to show the universality of rhythm-generation principles in such 3-cell circuits regardless of the cell model selected, which can be the HH-type model of the leech heart interneuron [22,23], the generalized Fitzhugh-Nagumo (gFN) model of neurons [24], or the minimalistic 2θ-burster suggested in this article, provided that all three meet some simple and generic criteria. We are convinced that one should first investigate the rules and mechanisms underlying the emergence of cooperative rhythms in basic neural motifs, as well as the role of coupling in generating a multiplicity of coexisting rhythmic outcomes in larger networks s [25].
The predecessor of a 2θ-burster proposed and examined below is the so-called spiking θ-neuron [26]. It is described by a phase differential equation with a specific term cosθ. The θ-neuron is meant to demonstrate a slow "quiescent" phase followed by a fast "spiking" transition. Mathematically, its equation is normal for a saddle-node bifurcation on a unit circle through which two equilibrium states, stable and repelling, merge and disappear. After the equilibriums are gone, the phase point keeps revolving on a unit circle (see Section 3 below). That is why this bifurcation is referred to as a homoclinic saddle-node bifurcation on an invariant circle or SNIC in short. The notion of the θ-neuron capitalizes on the feature of the saddle-node bifurcation, causing the well-known bottleneck effect, which results in slow quiescent and fast spiking time-scale dynamics in such systems.
The concept of the new model, called the 2θ-burster due to the driving term cos2θ in its ODE description, is inspired by the dynamics of endogenous bursters (like ones shown in Figure 1) with two characteristic slow phases: depolarized tonic spiking and hyperpolarized quiescent. These phases are also referred to as "on" or active and "off" or inactive depending on whether the membrane voltage is above or below some synaptic threshold. During the active phase, the presynaptic cell releases neurotransmitters to inhibit or excite other cells on the network, whereas during the inactive phase, the cell takes a pause from "communicating." This is a feature of chemical synapses that contrasts electric one or gap junctions, allowing cells interact all the time regardless of the voltage values. In contrast to interact the θ-neuron, there are two slow transient states, active and inactive, in the 2θ-burster due to two saddlenode bifurcations that alternate with fast progressions in between. We recall that a similar saddle-node bifurcation controlling the duration of the tonic-spiking phase, and hence the number of spikes is associated with the codimension-one bifurcation known as the blue-sky catastrophe [23,[29][30][31][32]].

RETURN MAPS FOR PHASE LAGS
We developed a computational toolkit for oscillatory networks that reduces the problem of the occurrence of bursting and spiking rhythms generated by a CPG network to the bifurcation analysis of attractors in the corresponding Poincaré return maps for the phase lags between oscillatory neurons. The structure of the phase space of the map is a unique signature of the CPG as it discloses all characteristics of the functional range of the network. The recurrence of rhythms generated by the CPG (represented by a system of coupled Hodgkin-Huxley-type neurons [23]) allows us employ Poincaré return maps defined for phaselags between spike and burst initiations in the constituent neurons [27,28], as illustrated in Figures 1, 2, and 6. With such return maps, we can predict and identify the set of robust outcomes in a CPG with mixed, inhibitory, excitatory, and electrical synapses, which are differentiated by phase-locked or periodically varying lags corresponding, respectively, to stable fixed points and invariant circles on the return map.
Let us introduce a 3-cell network ( Figure 1A) made of weakly coupled HH-like bursters; see the equations in the Appendix. Here, "weakly" indicates that coupling cannot disturb the shape of the stable bursting orbit in the 3D phase space of the individual HH model ( Figure 1A). Weak interactions, inhibitory (mainly repulsing) and excitatory/gap junction (mainly attracting), can only affect the phases of the periodically varying states of the neurons, represented by the color-coded spheres (blue/green/red, respectively, for cells 1/2/3), on the bursting orbit in the 3D phase space of the given interneuron model. As such, weak coupling can only gently alter the phase differences or phase lags between the networked neurons ( Figure 2A). We also note that coupling remains weak as long as individual cell models stay away from bifurcations, such as a saddle-node bifurcation typical for 1D 2θ-bursters. Making coupling stronger will make the convergence to the phase-locked states faster. However, in this study, we would like to demonstrate the reduction approach through which the 2D maps appear because they were defined analytically and not computationally. Otherwise, the trajectories will lose smoothness and look jagged and tangled.
Being inspired by neurophysiological recordings performed on various rhythmic CPGs, we employ only voltage traces generated by oscillatory networks to examine the time delays, τ 21 and τ 31 , between the burst upstrokes on each cycle in the reference/blue cell 1 and in cells 2 (green) and 3 (red). Next, we will show that like the biologically plausible HH-type networks, 3cell circuits of coupled 2θ-bursters can stably produce similar phase-locked rhythms. They include, for example, peristaltic patterns or traveling waves, in which the cells burst sequentially one after the other (see Figures 1 and 3C,E), and the so-called pacemaker rhythms, in which one cell effectively inhibits and bursts in antiphase with the other two bursting synchronously. The symmetric connectivity implies that such 3cell networks can produce multiple rhythms due to cyclic permutations of the constituent cells (see Figure 3 below).
To analyze the existence and the stability of various recurrent rhythms produced by such networks, we employ our previously developed approach using Poincaré return maps for phase lags between constituent neurons. We introduce phase lags at specific events in time when the voltage in cells reaches some threshold value, signaling the burst initiation (see Figure 1B). The phase lag Δϕ (n) 1j is then defined by a delay between nth burst initiations in the given cell and the reference cell 1, normalized over the bursting period: Sequences of phase lags {Δϕ (n) 12 , Δϕ (n) 13 } defined on module one represent forward trajectories M n on a 2D phase torus ( Figure 2B). The specific phase-lag values such as 0 (or 1) and 0.5 represent, respectively, in-phase and antiphase relationships of cells 2 and 3 with the reference cell 1. We examine the (Δϕ 12 , Δϕ 13 ) phase-lag structure of the 2D Poincaré return maps (such as one shown in Figure 3A) of the 3-cell networks by initiating multiple trajectories with a dense distribution of initial phase lags (50 × 50 grid) and by following their progressions over large numbers of cycles. In the long run, these trajectories can eventually converge to some attractors, one or several. Such an attractor can be a fixed point (FP) with constant values Δϕ * 12 and Δϕ * 13 in Eq. 1, which correspond to a stable rhythmic pattern with phase lags locked ( Figure 2A). All phase trajectories converging to the same fixed point are marked by the same color to reveal the attraction basins of the corresponding rhythms ( Figures 2B and 3A). This reduces the analysis of rhythmic activity generated by a 3-cell network on the examination of the corresponding 2D Poincaré map for the phase lags. For example, the map shown in Figure 3A reveals the existence of pentastability with the symmetric circuit generating three pacemaker (PM) rhythms and two, clockwise and FIGURE 2 | (A) The demonstration of the slow exponential convergence of initial states of Δϕ 21 (yellow curves) and Δϕ 31 (purple curves) to four phase-locked states: 0 ≡ 1, 1 3 , 1 2 , 2 3 in the inhibitory 3-cell motif (4) with weak coupling β i j 0.003. (B) Poincaré return map defined on a unit 2D torus, S 2 S 1 ⊗S 1 of the two phase lags (discrete values), revealing color-coded attraction basins of several fixed points (solid dots of the same colors) corresponding to the phase-locked rhythms generated by the 3-cell motif. A flattened torus is shown in Figure 3A.
Other attractors in the return maps on a 2D torus can be a stable invariant curve (IC) corresponding to rhythmic patterns with periodically varying phase lags. Such a curve can enclose a focal fixed point on the torus or wrap around the 2D torus [27,28] (see Figure 2B and Section 5.4 below). If the map has a single attractor, then the corresponding network is monostable; otherwise, it is a multifunctional or multistable network capable of producing several rhythmic outcomes robustly. Such as 2D return map, where Π : M n → M n+1 , for the phase lags can be represented as follows: with small μ i being associated with weak coupling; f i being some undetermined coupling function such that their zeros: f 1 f 2 0 correspond to fixed points: j1 of the map. These functions, similar to phase-resetting curves, can be numerically evaluated from all simulated trajectories {Δϕ (n) 21 , Δϕ (n) 31 } (see Figure 4C). By treating f i as partials zF/zϕ ij , one may try to restore a "phase potential"-some surface F(ϕ 21 , ϕ 31 ) C (see Figure 4). The f i quantities can be evaluated as the distance between two consecutive points (iterates) M n and M n+1 on every trajectory of the map (see Figure 3, for example). The shape of such a surface determines the location of critical points associated with FPs-attractors, repellers, and saddles of the map. With this approach, one can try to predict bifurcations due to landscape transformations and thus to interpret possible dynamics of the network as a whole. Figures 4A,B are meant to give an idea of how the potential surface may look like in the case of the 3-cell circuit with only two stable traveling wave patterns and in the case of three coexisting pacemakers only, respectively. Figure 4C shows a numerical reconstruction of the pseudopotential with the use of the distances between all pairs of successive iterates of the map with five stable FPs as depicted in Figure 3A.

MINIMALISTIC 2θ-BURSTER
The key feature of the 2θ-neuron given by which is the occurrence of two saddle-node bifurcations giving rise to the two slow transient phases in its dynamics alternating with fast transitions in between. Likewise endogenous bursters with two such slow states (Figure 1), the duration of the active Frontiers in Applied Mathematics and Statistics | www.frontiersin.org November 2020 | Volume 6 | Article 588904 tonic spiking, and the quiescent phases can be controlled independently in the 2θ-burster: the active "on" state and the inactive "off" state are due to the same bottleneck effects caused by the saddle-node bifurcations. This regulates the duty cycle of bursting, which is the fraction of the active-state duration compared with the burst period. As seen from Figure 5, the θ-model was meant to replicate phenomenologically fast-spiking cells, whereas the "spikeless" 2θ-neuron mimics burster dynamics instead. Next, we show that the network dynamics of a 3-cell motif of inhibitory coupled 2θ-bursters demonstrate that the key properties observed in such motifs are composed of Hodgkin-Huxley-type bursters. First, let us observe from Eq. 3 that the dynamics of the individual 2θ-burster is determined mainly by the driving term ω − cos2θ in symmetric α 0 and asymmetric cases due to small α-values. So, whenever the frequency 0 < ω ≤ 1, there exist two FIGURE 4 | Critical points of the sketched "pseudopotentials" with periodic boundary conditions reveal the location of potential wells-attractors, as well as saddles (including with six separatrices in (B)) and repellers (a single local maximum at the origin) in the (ϕ 21 , ϕ 31 ) phase space. These configurations correspond to the network with two traveling waves only (A) and with three pacemakers only (B), respectively. (C) A computational reconstruction of a pseudopotential/coupling function corresponding to the return map in Figure 3A.
FIGURE 5 | Comparison of the oscillatory dynamics generated by the spiking θ-neuron and the 2θ-burster. Panels (A) and (C) present snapshots of typical trajectories generated by both models on a unit circle S 1 (parametrized using Cartesian coordinates: x(t) sin(θ(t)) and y(t) −cos(θ(t)) with the origin 0 at 6 pm. (A) Clustering of purple spheres near the origin is due to a post-effect posteffect caused by a saddle-node bifurcation (SNIC) in the θ-model, whereas the 2θ-burster in (C) features two such bottleneck post-effect due to two heteroclinic saddle-node connections causing the stagnation of gray spheres near the top, "on" state and the inactive "off" state of the symmetric 2θ-burster and fast transitions in between. (B) Spiking trace (purple) of the θ-neuron, overlapped with 2-plateau traces of the 2θ-neuron with three values of the bursting duty cycle: x50%, 30%, and 70% (solid, short-, and long-dashed gray curves, respectively) controlled by small variations of the α-parameter in Eqs. 3 and 4.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org November 2020 | Volume 6 | Article 588904 pairs of stable and unstable equilibriums: one pair is near the bottom θx0 and the other is at the top around θxπ. The stable equilibria are associated with the hyperpolarized active and depolarized quiescent states of neurons, respectively.
Increasing ω > 1 makes the 2θ-neuron oscillatory through two simultaneous (if α 0) saddle-node bifurcations (SNIC) on a unit circle S 1 , which is its phase space. Moreover, as long as ω 1 + Δω, where 0 < Δω ≪ 1, the 2θ-burster possesses two slow phases: the active "on" state near θ π and the inactive "off" state near 0 on S 1 . These slow phases are alternated with fast counterclockwise transitions, which will be referred to as an upstroke and a downstroke, respectively. For greater values of ω, the active and inactive phases are defined more broadly they are: π/2 < θ ≤ 3π/2 and 3π/2 < θ ≤ π/2, respectively. This is convenient because the inactive phase remains below the synaptic threshold, which is set at θ th π/2 so that cosθ th 0 for the sake of simplicity, thus equally dividing the unit circle (see Figure 6A). The duty cycle of the 2θ-burster is controlled by the axillary term αcosθ in Eq. 3, provided that it remains oscillatory as long as ω − |α| > 1. Note that when α 0, the duty cycle of bursting is 50%, and the corresponding traces have two even plateaus (see Figure 5B). The active or inactive phases can be extended or shortened with small variations of the α-parameter so that α < 0 increases the duty cycle and α > 0 decreases the duty cycle of the individual 2θ-burster when it is shifted, respectively, closer to the top or to the bottom saddle-node "phantom," because the bottleneck effects become more profound, see Figure 5B,C.

THREE EQUATIONS FOR 3-CELL NETWORK
A 3-cell circuit of the 2θ-bursters coupled with chemical synapses is given by the following system: The 2θ-burster is coupled in the network using the fastinhibitory synapses driven by the fast-threshold modulation [33]. It is due to the sigmoidal term 1 1+e kcosθ i that, rapidly (here k 10) varying between 0 and 1, triggers an influx of inhibition flowing from the presynaptic neuron to the postsynaptic neuron, as soon as the former enters the active "on" phase above the synaptic threshold cosθ th 0, that is, π/2 < θ i < 3π/2. Note that the negative sign of this term makes the synapse inhibitory so that replacing it with "positive sign" makes the synapse excitatory because it would increase the rate of θ ′ during the upstroke, in contrast to slowing down the upstroke in the inhibitory case. The strength of the coupling is determined by the maximal conductance value β ij . Depending on the magnitude of βij values, the active cell in the "on" state can either slow down the inactive postsynaptic one due to the bottleneck effect (weak coupling) or shut it down completely with the saddle-node bifurcation in its perturbed state (strong coupling), which can be referred to as soft vs. hard locking, respectively. If the postsynaptic cell happens to be at the active phase, then the inhibition FIGURE 6 | (A) Sampling the moments in phase traces, y i (t) −cos(θ i (t)), plotted against time, when they reach a synaptic threshold θ syn 0, defines a sequence of the phase lags (τ (n) 21 , τ (n) 31 ) between upstrokes in the reference (blue) neuron 1 and the other two 2θ-neurons coupled in the 3-cell network. (B) Parametric representation of the 1D phase space of the coupled 2θ-bursters traversing counterclockwise (long gray arrows indicating rapid transition between "on and off" states) on a unit circle S 1 . Small downward blue and red arrows illustrate the inhibitory perturbations projected from the active green cell above the synaptic threshold that delays the forthcoming upstroke of the blue cell and speeds up the red cell toward the inactive phase. The gray arrows indicate the direction and its speed on a unit circle.
will shorten its duration significantly provided that βij is large enough. We deemphasize that it is the closeness to the saddle-node bifurcations in the postsynaptic cells that determines whether the coupling is weak or strong. Our coupling strategy is to ensure that θ ′ i > 0 in all three equations in Eq. 4, that is, the cells maintain endogenous bursting in isolation and on the network and converge to the phase-locked states exponentially (Figure 2A). This does not have been the case case. By increasing increasing βij and α or by decreasing bursting frequency ω or by manipulating all three parameters, one can speed up the convergence substantially ( Figure 7D) or even make the network rapidly reach any steady state in one or two steps [21].
The last term 1 − 2 1+e ksinθ , breaking the symmetry of coupling on upstrokes and downstrokes, converts the synaptic input into qualitative inhibition. Namely, its sign is switched from + to − upon crossing the values θ 0 and θ π. During the fast upstroke, when 0 < θ < π, this term is positive, thereby ensuring that inhibition slows down or delays the transition into bursting. When π < θ i < 2π, during the fast downstroke, this term 1 − 2 1+e ksinθ < 0 to ensure that the inhibition speeds up the transition from the active (tonic-spiking) phase of bursting into the inactive (quiescence) phase. This is phenomenologically consistent with neurophysiological recordings because inhibition projected onto the postsynaptic burster typically shortens the burst duration and extends the interburst intervals. Alternatively, this term can be replaced with 1 − 1 1+e ksinθ as it not only breaks the symmetry but also only acts during the upstroke of bursting.
The electrical coupling or the gap junction between the neurons is handled by the other term −C elec sin(θ pre − θ post ). It slows down the rate θ ′ post when θ post > θ pre and speeds it up if θ post < θ pre . The conductivity coefficient C elec is to be set around two orders of magnitude smaller than β values to maintain a balanced effect in the network. When C elec and β are of the same magnitude, the dynamics of network are solely dictated by the electrical coupling with the inhibitory synapses insignificantly affecting it.
Let us note that unlike the bidirectional electrical synapse, a chemical one is unidirectional and hence asymmetric because it has a synaptic threshold: the chemical synapse becomes functional when the membrane voltage in the presynaptic cell rises above the only synaptic threshold in the active phase; otherwise the synapse is silent. This is the reason why a network composed of identical cell, and identical chemical synapses can only be called symmetric loosely, in some permutation sense. This is the reason why a permutationsymmetric 3-cell network always possesses a pair of traveling wave patterns (stable or not) where the cells burst sequentially and/or may generate three pacemaker patterns where one cell bursts in antiphase with the two others. Note that the last rhythms cannot be produced by a properly symmetric network by default. Figure 6A shows how phase lags between upstrokes are introduced (here, cell 1 (blue) is the reference) between the three-networked 2θ-bursters turning counterclockwise on the unit circle S 1 ( Figure 6B). It is observed that inhibition generated by the green cell 2 in the active slow phase near θ π above the synaptic threshold (given by cosθ th 0) brings the other two cells closer to the bottom quiescent state bear θ 0 by accelerating the red burster 3 on the downstroke and by slowing down the blue burster 1 on the upstroke.

POINCARÉ RETURN MAPS FOR THE PHASE LAGS. RESULTS
Following the same approach used in the weakly coupled HHtype models above, we first create a uniform distribution of initial phases on S 1 , and therefore, the phase lags between the three 2θ-bursters. Next, we integrate the network (4) over a large number of cycles and record burst initiations (see Figure 5A) to determine the phase lags between the reference cell 1 and two other cells and determine the kind of phase-locked states they can converge. This approach is illustrated in Figure 2A for the symmetric 3-cell motif composed of identical 2θ-bursters and equal inhibitory synapses. The corresponding 2D Poincaré return map, with the coexisting stable fixed points and saddles is shown in Figure 3. By stitching together the opposite sides of this map, we wrap it around a 2D torus as shown in Figure 2B.
The fixed points and their attraction basins are coded with different colors in the map. For example, the Poincaré return map for the (Δϕ 21 , Δϕ 31 ) phase lags represented in Figure 3A has five stable fixed points representing robust three pacemaker FPs: red at 0, 1 2 , green at 1 2 , 0 , and blue at 1 2 , 1 2 and two traveling wave rhythmic patterns: yellow clockwise at 1 3 , 2 3 and teal counterclockwise at 2 3 , 1 3 . The borders of the attraction basins of these five FPs are determined by positions of stable sets (separatrices) of six saddles (gray dots). The origin is a repelling FP. In total, there are eight hyperbolic FPs in this Poincaré return map.
Let us underline another handy feature of the 2θ-burster paradigm. We can easily detect and explore repelling FPs or invariant circles, if any, existing in the 2D Poincaré map by reversing the integration direction of system (4), that is, multiplying the right-hand sides by −1, simulating the network in backward time. This reverses the direction to spin trajectories clockwise on S 1 , whereas the backward time integration will make dissipative systems run to infinity.

Homogeneous Motif With Identical Cells and Synapses
It will be shown below that 2θ-bursters weakly coupled in 3-cell networks, whether they are homogenous/symmetric or nonhomogeneous/asymmetric, can generate the same stable rhythms as the networks of biologically plausible HH-type models. We also discuss the bifurcations occurring in the networks and corresponding maps as synaptic connectivity and intrinsic temporal characteristics of the 2θ-bursters are changed. Bifurcations in the system are identified and classified by a change in the stable phase rhythms, which can be due to the stability loss of a particular FP or when it merges with a close saddle so that both disappear through a saddle-node bifurcation. Let us first consider a symmetric network with two bifurcation parameters: the coupling strength β β ii (i 1, 2, 3) and the α-parameter in Eq. 3 that controls the duty cycle (DC) of the 2θ-bursters. We use five different DC values as α is varied from −0.11 to 0.11l while synaptic strength is increased through four steps from β 0.0001 through β 0.1. The results are presented in Figure 7. The Panels A2/3 represent the most balanced, weakly coupled network that can produce all five bursting rhythms with 50% DC. One can see that by increasing the β-value, the saddles separating 2 TWs and 3 PMs move toward the latter ones, and after some critical values, three pairs, a saddle and the nearest stable PM, merge and vanish simultaneously. After that, the symmetric network can produce two only rhythms: counterclockwise and clockwise TWs corresponding to the teal and yellow stable FPs at 1 3 , 2 3 and 2 3 , 1 3 , respectively. This would correspond to the case of the "pseudopotential" depicted in Figure 4A.
The stable PMs promote or dominate the dynamics of the symmetric case at extreme α-values corresponding to the bursting rhythms with short or long burst durations. Once we compare panels, say A1 and D4 reveal this time, the separating saddles group around the stable TWs to minimize their attraction basins, and hence the likelihood of the occurrence of these rhythms in the network. These cases correspond to the "pseudopotential" depicted in Figure 4B.

"Winner Takes all" Motif
The first asymmetric case considered is a motif termed the "winner takes it all." In this modeling scenario, both outgoing inhibitory synapses from the given cell, here the reference blue burster 1, are evenly increased in the strength, see Figure 8A. It is observed that such a configuration breaks down the circular (and permutation) symmetries supporting traveling waves in the network. Let us start with Figure 8B: no surprise that with an initial increase in β 1,2/3 , two saddles shift away from the blue PM at (0.5, 0.5) toward 2 TWs, then merge with them to disappear pairwise. Next, as β 1,2/3 is FIGURE 9 | Mono-biased network motif (F) with one different synapse due to increasing β 21 . (A) The first of five (Δϕ 21 , Δϕ 31 ) return maps, an increase in β 21 value breaks down a counterclockwise symmetry so that the attraction basin (teal) of the corresponding TW at 2 3 , 1 3 shrinks as a nearby saddle moves closer to it and away from the green PM at 1 2 , 0 (A and B). (C) With further increase in β 23 , the counterclockwise TW at 2 3 , 1 3 vanishes through a saddle-node bifurcation after merging with the nearest saddle, followed by another saddle-node bifurcation eliminating the red PM at (0, 0.5) (D). At greater β 23 values, the green PM at 1 2 , 0 encompasses the majority of the network phase space, along with the blue PM at 1 2 , 1 2 , preserving the size of its attraction basin. The parameters are ω 1.15, α 0.07, and β′s 0.003 except β 21 0.00042, 0.0045, 0.01, and 0.02 for panels (A-D).
increased further, two other saddles annihilate the green and red PMs through similar saddle-node bifurcations ( Figure 8C). In the aftermath, the 3-cell network with a single burster generating the repulsive inhibition much stronger than the other two cells becomes a monostable one producing a single pacemaking rhythm with the phase lags locked at (0.5, 0.5).

Mono-Biased Motif
We refer as a mono-biased motif to another asymmetric network with a single different synapse. In this case, the strength β 21 of the outgoing synapse from cell 2 to cell 1 is increased, which violates the circular symmetry supporting the counterclockwise traveling wave in the network, see Figure 9F. So, as β 21 is increased, the counterclockwise stable FP at 2 3 , 1 3 first disappears through a saddle-node bifurcation, as seen in Figures 9A,B. Because this was the saddle between the TW and the green PM, the attraction basin of the latter increases after the first bifurcation in the sequence. The next saddle-node bifurcation eliminates the red stable FP at (0, 0.5). The reasoning is as follows: for this rhythm to persist, the red PM should evenly inhibit both green and blue PMs. However, a growing inhibition imbalance between them is no longer reciprocal. As we pointed out earlier, the stronger inhibition from cell 2 shortens the active phase of the blue burster. As so, they cannot longer line up by the burster 3, which causes the disappearance of this PM-rhythm and the FP itself ( Figure 9C). Similar arguments can be used to justify the disappearance of the green PM as cell 2 cannot inhibit cells 1 and 2 evenly to hold them together as β 21 is increased further (not shown). This is in this case is in good agreement with the 3-cell networks of the HH-type bursters. Frontiers in Applied Mathematics and Statistics | www.frontiersin.org November 2020 | Volume 6 | Article 588904

Dedicated HCO
The abbreviation HCO stands for a half-center oscillator, where a pair of neurons coupled reciprocally by inhibitory synapses to produce alternating bursting. Such a dedicated HCO is formed by cells 2 and 3 with stronger synapses due to β 23 β 32 in the configuration shown in Figure 10C. Again with start off with the symmetric case depicted in Figure 10A, one can observe at once that having the dedicated HCO should break down the circular symmetries of the network. So, the stable TWs become eliminated first as β 23 β 32 starts increasing. As these synapses become stronger, the attraction basin of the blue PM at (0.5 0.5) shrinks substantially, but the FP itself persists. Meanwhile, increasing β 23 β 32 further creates the inhibitory imbalance that makes the further existence of the green and red PMs impossible due to the factors that we outlined above for the mono-biased motif. Both vanish at the same time due to saddle-node bifurcations. However, at the bifurcation, both double FPs are connected by a heteroclinic orbit that transforms into a stable invariant curve wrapping around the torus (see Figure 10F). This stable invariant curve is associated with a phase-slipping rhythm that recurrently passes slowly through the "ghosts" of all four vanished FPs except for the coexisting blue PM, see the fragments of the corresponding traces presented in Figure 10G.

Clockwise-Biased Motif
The clockwise-biased motif in this case represents the 3-cell network with counterclockwise connections stronger than ones in the clockwise direction, see Figure 11E. This configuration does not break the circular symmetries of the network but implies that either TW should win over the opposite one, which should result in their attraction basins changing correspondingly. Figure 11 presents four transformation stages of the map as β 13 , β 32 , and β 21 sequentially increased. With a small increase, the shape of the map becomes slightly twisted with the three saddles shifting away from the stable PMs toward the teal TW at 2 3 , 1 3 . A further increase brings the saddle close to the teal one, thereby shrinking its attraction basin and substantially widening the basin of the clockwise TW at 1 3 , 2 3 . Finally, as some bifurcation threshold is reached, the saddles collapse at the stable FP that becomes a complex saddle with three outgoing and three incoming separatrices. This means that the counterclockwise TW becomes an unstable rhythm in such biased 3-cell motif that is fully dominated by the clockwise TW rhythm.

Gap Junction
In our last example, we consider the symmetric motif with a gap junction or electric synapse added between cells 1 and 2 as shown in Figure 12C. Recall that a gap junction is bidirectional unlike unidirectional chemical synapses with synaptic thresholds. Recall that it is modeled by this term −C elec sin(θ pre − θ post ) that slows down the rate θ ′ post when θ post > θ pre and speeds it up if θ post < θ pre . Due to this property, the electrical synapses, like excitatory ones, FIGURE 11 | (E) Clockwise-biased motif with three synaptic strengths, β 13 , β 32 , and β 21 sequentially increased. (A) As all three counterclockwise synapses are slightly strengthen, saddles shift away from the three stable PMs, blue at 1 2 , 1 2 , green at 1 2 , 0 , and red at 0, 1 2 , toward the teal clockwise TW at 2 3 , 1 3 (B) thus shrinking its basin and widening the attraction basin of the dominant counterclockwise TW (yellow) at 1 typically promote synchrony between such coupled oscillatory cells, as in our case between cells 1 and 2.
It is observed that introducing an electrical synapse between only two of the cells of the motif breaks down both circular symmetries in the system. This is documented in Figures 12A,B depicting the maps for the networks with C elec being increased from zero to 0.0003. One can see that both TWs vanished from the repertoire of the network. Further increase in C elec makes the stable green and blue stable PMs disintegrate as both cells become synchronous and burst in alternation with the red cell 3. This completes the consideration of the monostable network with a relatively strong gap junction between cells 1 and 2 that can only produce the only one pacemaker rhythm.

CONCLUSION
Our ultimate goal is to use the top-down approach to single out the fundamental principles of the rhythm formation in small networks that can be systematically generalized and applied for understanding larger network architectures. Due to the rhythmic nature of the bursting patterns, we employed Poincaré return maps defined on phases and phase lags between burst initiations in the interneurons. These maps allow us to study quantitative and qualitative properties of the stable rhythms and their corresponding attractor basins. The specific goal of this study is to demonstrate the simplicity and usability of the 2θ-bursters to construct multistable, polyrhythmic neural networks that have the same dynamical and bifurcation properties as ones composed of biologically plausible models of Hodgkin-Huxley-type bursters and chemical synapses [34,35].We argued that the maps derived from the HH-type bursters and ones from the 2θ-bursters have the same structure. As such, these maps serve as a detailed blueprint containing all necessary information about the network in question, including its rhythmic repertoire, stability of generated patterns, and the like. In addition, with such maps, one can predict possible transformations before they occur in the network. Furthermore, we showed that depending on strengths of inhibition, the maps and hence the corresponding networks may have different distributions of phase-locked states. As such, the proposed approach reveals the capacity of the network and the dependence of its outcomes on coupling strength, wiring circuitry, and synapses, thereby allowing one to identify necessary and sufficient conditions for rhythmic outcomes to occur. Our study is a further step toward the foundation of the bifurcation theory of multifunctional rhythmic circuits including network with a modular organization of subcircuits [25].
In this article, we did not discuss the rhythm-generating motifs composed of 2θ-cells that are quiescent in isolation. While motifs, made of coupled 2θ-cells initially placed at the upper "on" state, functions well using the escape mechanism for the rhythmogenesis, the other basic mechanism based on the postinhibitory rebound [24] is not (fully) applicable to 2θ-cells because it requires at least two dynamical variables, slow and the fast, to warrant the occurrence of specific transient dynamics in the system.
Our computational approach based on the reduction to the evident Poincaré return maps for phase lags extracted from voltage traces were inspired by neurophysiological recordings from biological CPGs, such as the 3-cell pyloric one and swim CPGs of sea slugs. The predictive power of the map approach is that it allows constructing a desired neural circuit with some preset properties. With such maps, one also gains generalizable insights helpful for the better understanding of the fundamental and universal rules of the pattern formation in various models of central pattern generators. Our findings can be employed for identifying or implementing the conditions for normal and pathological functioning of basic CPGs of animals and artificially intelligent prosthetics that can regulate various movements.