We disclose the generality of the intrinsic mechanisms underlying multistability in reciprocally inhibitory 3-cell circuits composed of simplified, low-dimensional models of oscillatory neurons, as opposed to those of a detailed Hodgkin–Huxley type [Wojcik *et al*., PLoS One **9**, e92918 (2014)]. The computational reduction to return maps for the phase-lags between neurons reveals a rich multiplicity of rhythmic patterns in such circuits. We perform a detailed bifurcation analysis to show how such rhythms can emerge, disappear, and gain or lose stability, as the parameters of the individual cells and the synapses are varied.

Complex multistable bursting rhythms can arise even in simple biological neural circuits. We employ a computational technique combining phase-lag return maps and fast (parallel) Graphics Processing Unit (GPU)-based sweeps of the phase and parameter spaces to identify multistable patterns, rhythm switching, and attractor robustness in reciprocally inhibitory 3-cell circuits composed of the proposed generalized Fitzhugh–Nagumo model of “spike-less” bursters. With such maps, we can thoroughly examine how internal and external factors such as the synaptic strengths, network asymmetry, and externally injected currents determine what stable rhythmic patterns can co-exist, emerge, or disappear, as well as their underlying bifurcation mechanisms. Depending on intrinsic mechanisms, such as release and escape in individual cells, these networks can produce a variety of multistable rhythmic states, ranging from penta-stability with phase-locked bursting pacemakers and traveling wave patterns, or stable chimeras admitting recurrent phase slipping (PS) of one cell with respect to the other two that remain phase-locked over time, to more exotic behaviors such as a robust stable synchronous state with all three cells oscillating together or a lack of phase locked rhythmic states altogether. We present detailed transition mechanisms between such rhythms including saddle-node, pitch-fork, and secondary Andronov–Hopf or torus bifurcations, as well as the emergence of a transitive torus. Last, we introduce the concept of $2\theta $-neurons to build even simpler neural circuits capable of desired dynamics. Our qualification promotes the use of simplified, low-dimensional modeling of multistable bursting patterns arising in oscillatory neural circuits in lieu of computationally intensive high-dimensional Hodgkin–Huxley type models.

## I. INTRODUCTION

A central pattern generator (CPG) is a small network of coupled neurons that determines and autonomously controls rhythmic oscillations underlying sensory, motor, or cognitive behaviors of an animal. CPGs are implicated in a variety of functions ranging from respiration, heartbeat, and circulation to sleep and locomotion.^{2–15} While many insights into the operational principles of CPGs have been obtained from experimental and computational studies, the basic principles of robustness and stability of many CPGs observed in nature are yet poorly understood and cannot be inferred *a priori.* The cooperative dynamics of coupled cells is an area of active ongoing research, with both biological and phenomenological approaches employed.^{8,16–23} The smallest building unit commonly tied to many CPGs is a pair of bilaterally symmetric neurons that mutually inhibit each other to produce anti-phase bursting, called a *half-center oscillator* (HCO). The current study is focused on the rhythmic dynamics, transitions, and bifurcations occurring in the context of a 3-cell neural network motif, made up of interconnected HCO circuits. Various 3-cell biological circuits that form constituent blocks or centers for larger networks have been reported previously.^{24–30}

Several modeling paradigms have been applied for studies of such circuits, including biologically relevant Hodgkin–Huxley (HH) type models, in which the individual parameters can be related to specific ionic currents or concentration gradients. The high dimensionality of the detailed HH-type models presents obvious difficulties in performing a thorough dynamical and bifurcation analysis to classify their generic properties. Such understanding is essential to reliably assemble and configure small neural networks with common oscillatory characteristics. The simpler integrate-and-fire models, which belong to the opposite end of the spectrum of mathematical models, are often inadequate to connect their parameters to biological mechanisms that may be directly manipulated or affected, failing to capture nuances in the dynamic behaviors that intermediate systems can. In this paper, we employ the so-called generalized FitzHugh–Nagumo (gFN) neurons to model 3-cell networks. This 2D gFN-neuron model presents a better description of some of the pivotal properties of typical HH-type square wave bursters. We use it to showcase the essential characteristics of the building blocks of rhythmic circuits to stably generate the desired dynamics, regardless of the specific models of neurons and synapses. The gFN-equations are simpler and hence more practical for computational studies, especially for intense GPU-based sweeping of parameters and initial conditions.

This paper capitalizes on our previous work and the well established principles in the characterization of 3-cell circuits made of HH-type neurons (as depicted in Fig. 1).^{1,33,34} The goal is to present novel findings using parameter sweeps that reveal how the multistability of a 3-gFN motif is aligned with its parameter space. It also serves as a tutorial blueprint, providing a complete framework for interested researchers to borrow and employ our methods for the analysis of similar oscillatory networks. We employ gFN-neurons to study a variety of polyrhythmic dynamics arising in 3-cell circuits with symmetric and asymmetric connectivity, as well as various fast, slow, and delayed effects. While maintaining generic behaviors, this reduction in complexity allows for more extensive exploration of key parameters in neural systems ranging from inherently quiescent or tonic spiking to intrinsic bursters. It also aids our search for biologically plausible circuitries that ensure the robustness of the rhythmic patterns observed in nature. One important aspect of the so-called multifunctional CPGs is the ability of the same circuit to produce more than one observable rhythmic outcomes, as well as to switch between its rhythms.^{29,33,35–40} We examine how changes in parameters can trigger transitions or bifurcations in the rhythms of an otherwise robust network. We summarize the use of phase lags, their corresponding Poincaré return maps, and fast GPU-based parallel simulations to characterize the state space of the 3-cell motif and to identify its stable and unstable rhythms.

Specific results are presented in terms of two primary mechanisms underlying network rhythmogenesis—synaptic release and escape.^{39,41} For each mechanism, we contrast the behaviors of symmetric networks with uniform all-to-all inhibitory connections, to those observed in asymmetric motifs with one or several distinct connections. We study how network circuitry influences the occurrence and stability of various rhythms with phase-locked states, periodic phase slipping, or chimera-like behaviors. We determine the ranges of network parameter values, specifically the synaptic coupling strengths and the external current drives, that give rise to these behaviors and the underlying bifurcations. We also connect overarching mechanisms and present the stages of transition from the release to the escape mechanisms, in terms of bifurcation sequences that the corresponding fixed points (FPs) in the Poincaré return maps for the phase lags between the constituent neurons may undergo due to network parameter variations. We elaborate on several specific examples chosen for both relevance and novelty.

This paper is organized as follows: first, we discuss the gFN models of neurons coupled in a network with fast inhibitory synapses. Next, we introduce phase lags between oscillatory neurons, followed by how 2D Poincaré return maps for such phase lags are defined and the correlations between the fixed points of these maps and the multiple phase-locked states of the corresponding voltage traces. We examine several representative networks to demonstrate the association between the repertoire of possible rhythmic outcomes for a given motif and its parameter space. We also discuss the role of non-local bifurcations and how these shape the borderlines separating the attraction basins in the 2D Poincaré map. This is followed by a discussion of some exotic motifs without any phase-locked rhythmic states. Also of our special interest is the case of the motif that opens up another possibility—the robust synchronous state when all three cells oscillate together. Finally, we introduce the concept of the reduced $2\theta $-neurons to build even simpler multifunctional oscillatory motifs.

## II. METHODS

Fitzhugh–Nagumo-like cells in biological sciences, or generic relaxation oscillators, present a mathematical generalization that captures some plain dynamical features often observed and reported in detailed HH-type models. The generalized Fitzhugh–Nagumo (gFN) model of neurons employed in this paper adds a set of extra dynamical and temporal features to more realistically replicate biological (endogenous) bursters in isolation and, what is more important, under perturbations. We consider a fully connected 3-cell circuit of such gFN-neurons,^{39}

Here, the state of the $i$th neuron is described by its activity variable $V$, representing the membrane voltage, and a recovery voltage-gated variable $h$, introduced in a way similar to the Hodgkin–Huxley formalism; $\epsilon $ is the reciprocal of some time constant and regulates slow–fast dynamics ($0<\epsilon <1$) in the gFN-neuron [see voltage traces in Figs. 2(c) and 2(d)], while an applied current $Iapp$ is used as a bifurcation parameter for individual cells; $V0$ and $k$ influence the position and shape of the cubic and sigmoidal nullclines given by $V\u02d9=0$ and $h\u02d9=0$, respectively [see Fig. 2(b)]. The default values for the parameters are $k=10$, $\epsilon =0.3$, and $Vth=0$. By construction, active ($Vi>0$) driving or pre-synaptic neurons slow down or repress the recovery dynamics of the driven or post-synaptic oscillators via a fast inhibitory coupling given by $Gji$, modeled using a sigmoidal coupling function employed via fast-threshold modulation (FTM),^{42}

The FTM-formalism can sharply differentiate the active “on” state of the neuron, when its voltage $Vj$ is above the synaptic threshold $Vth=0$, and hence, $\Gamma (Vj)=1$, from the inactive “off” state with $\Gamma (Vj)=0$, when $Vj<Vth$, provided that $\Gamma $ is stiff enough (due to the factor 100 in its equation). The strength of coupling is controlled by the maximal conductance $gji$; its default value is set to $0.001$, unless otherwise specified, to ensure weak coupling in the network. The choice of $Vth=0$ and $Vrev=\u22121.5$ that makes $(Vrev\u2212Vi)<0$ defines the inhibitory synapse projected from the active neuron $j$, with $\Gamma (Vj)=1$, to neuron $i$, which slows down the rate $V\u02d9i$ in Eq. (1). Geometrically [see Fig. 2(b)], the term $gji(Vrev\u2212Vi)$ shifts and skews the fast nullcline $V\u02d9i=0$ closer to the slow nullcline $h\u02d9i=0$ of the inhibited post-synaptic neuron in the phase space. As depicted in Fig. 2(b), an individual gFN-neuron has a single unstable equilibrium state at the intersection of the fast cubic $V$-nullcline ($V\u02d9=0$) and the slow sigmoidal $h$-nullcline ($h\u02d9=0$), surrounded by a stable limit cycle producing the relaxation oscillations shown. Note that strong inhibition or negative applied current $Iapp$, which, respectively, induces temporary or permanent shifts of the $V$-nullcline, can make it cross the $h$-nullcline near the lower knee to give rise to another stable equilibrium on the lower-hyperpolarized stable branch of the cubic nullcline, at which the inhibited neuron will rest (described further in Fig. 6). If inhibition is not strong enough, it slows down the recovery or, equivalently, extends the inactive “off” state of the post-synaptic neuron by bringing the lower knee of the fast $V$-nullcline closer to the sigmoidal $h$-nullcline and narrowing the gap between them. This is in accordance with the famous bottleneck effect of the saddle-node bifurcation, where the dwelling speed, estimated as $[V\u02d9]2+[h\u02d9]2$, decreases further closer the nullclines approach to tangency.

In what follows, we will show that 3-cell gFN networks can produce various stable phase-locked rhythms including the traveling waves, in which the cells fire sequentially one after the other [see Figs. 2(c), 2(d), and 3(b)], as well as the pacemakers, in which one cell effectively inhibits and fires in anti-phase with the remaining pair [see Fig. 3(c)]. The symmetric connectivity in this network implies the coexistence of multiple rhythms that result from cyclic permutations or relabeling of the cells. In order to analyze the stability of various recurrent rhythms produced in a network, we employ the approach of Poincaré return maps. First, we introduce the notion of phase lags between the constituent cells, defined at specific events in time when the cells cross the threshold voltage from below, signaling the burst initiation. The phase lag of a cell is then defined as the delay in its burst initiation with respect to that of the reference cell 1, normalized over the bursting period. We define the $n$th phase lags $\Delta 12(n)$ and $\Delta 13(n)$ of cells 2 and 3, respectively, as follows:

where $\tau in$ represents the time at which the $i$th cell reaches the threshold voltage, $Vth$=0, for the $n$th time [see Figs. 3(b) and 3(c)]. The sequence of phase lags ${\Delta 12(n),\Delta 13(n)}$, defined for values between $0$ and $1$, gives the forward phase trajectory on a 2D torus [Fig. 3(a)]. The specific phase lag values $0$ (or $1$) and $0.5$ represent in-phase and anti-phase relationships, respectively, with the reference cell 1. We examine the 2D phase space of the Poincaré return maps [Fig. 3(a)] of the 3-cell networks, given by the phase lags $\Delta 12$ and $\Delta 13$, by initiating multiple trajectories with different initial phase lags (on a grid of size $50\xd750$) and by following their evolutions for a large number of cycles. As we compute long traces of firing activity of the circuit and evaluate the corresponding phase lag trajectories, these phase lags eventually converge to some attractor, which can be a fixed point of the 2D map [with steady coordinates $\Delta 12\u2217$ and $\Delta 13\u2217$ in (3)], which implies the existence of a stable rhythmic pattern in the circuit with phase-locked bursting between the cells. All the phase trajectories that converge to the same fixed point are marked by the same color and reveal the basin of attraction for the corresponding stable rhythm. The rhythmic characteristics of a CPG can be identified by analyzing the phase space of the corresponding Poincaré map, such as Fig. 3(a), which reveals the existence of penta-stability in the circuit, composed of three pacemakers [blue—Fig. 3(c$2$), green, and red] and two traveling waves [clockwise—Fig. 3(c$1$) and anti-clockwise]. We depict stable fixed points with colored dots, saddles with gray diamonds, and unstable fixed points with white dots in the 2D maps. The phase lag values ($\Delta 12$, $\Delta 13$) at the fixed points (FPs) corresponding to the blue, green, and red pacemakers (PMs) are given by $(0.5,0.5)$, $(0.5,0)$, and $(0,0.5)$, respectively, while those for the clockwise (black) and anti-clockwise (purple) traveling-waves (TWs) are given by $(0.33,0.67)$ and $(0.67,0.33)$, respectively.

Alternatively, it is also possible that the forward phase-lag trajectories ${\Delta 12(n),\Delta 13(n)}$ could converge to attractors other than FPs, as demonstrated in Figs. 15(d$\u2032$), 16(f), and 20(a), 20(b), and 20(e), or the map has no attractors at all as in Fig. 24. One can also trace down the backward trajectories, as shown in Fig. 24, to effectively depict repelling fixed points and unstable phase-slipping orbits.

Numerical integration of the trajectories is performed using the fourth order Runge–Kutta method, with a fixed step size. Computation of voltage and phase lag trajectories across multiple initial conditions is parallelized on a Tesla K40 GPU using CUDA, while visualizations are performed in Python. The open source software toolkit developed is available at https://github.com/jusjusjus/Motiftoolbox.^{43} GPU parallelization allows the construction of typical phase sweeps such as one shown in Fig. 3(a) within just a few seconds.

With detailed simulations of the Poincaré return maps for phase lags, we can visualize and analyze their stable and unstable fixed points (FPs) and other limiting sets such as invariant circles (ICs). We can also detect various bifurcations including homoclinic and heteroclinic and, therefore, can identify multiple possible rhythmic patterns generated by the neural circuits. With changes in the bifurcation parameters of the network, the constructed phase-leg sequences also vary, thus allowing us to determine the basins of the coexisting attractors and to reveal bifurcations through which FPs can emerge, disappear, or lose their stability. By varying two bifurcation parameters of the circuit (typically, the external drive $Iapp$ and the coupling strength $g$) and identifying all the permissible rhythmic states, we construct bi-parametric sweeps such as one shown in Fig. 4(e). The given bifurcation diagram frames the neighborhoods of the release and the escape mechanisms (described later) in the parameter plane for the symmetric network [depicted in Fig. 2(a)] over a wide behavioral range of the individual neurons. We can characterize the parametric regions whose phase-lag return maps contain only the pacemakers (PM), only the traveling waves (TW), or a combination of both PM/TW. The return maps shown in panels (a1)–(a4) and (a1)–(d1) in Fig. 4 reveal the sequence of bifurcations that stable rhythms undergo near the borderlines. Such bifurcation diagrams have proven useful for studying the dynamics of small CPG circuits and other nonlinear systems.^{44–48} In addition to the fully symmetric motif in Fig. 2(a), we perform a detailed bifurcation analysis of some other key asymmetric motifs (see Fig. 5). They include (1) mono-biased motif, in which only a single synaptic connection is varied, while all others are held constant; (2) double-biased motif, in which the reciprocal connections between two cells are manipulated equally, while all others are held constant; (3) biased-driver motif, in which both the outgoing connections from one cell are varied identically, while the rest are held constant; and (4) clockwise-biased motif, in which all the clockwise connections are changed simultaneously, while the anti-clockwise connections remain fixed. This is followed by a brief description of another asymmetric network configuration (see Fig. 24) without any phase-locked rhythmic states but only quasiperiodic phase slipping or an ergodic flow with neither FPs nor ICs.

## III. RESULTS AND DISCUSSION

The terms “release” and “escape” mechanisms referring to anti-phase oscillations in 2-cell networks were first introduced by Wang and Rinzel.^{41} By construction, the release mechanism requires that isolated cells are intrinsically bursting, while the escape mechanism requires the cell to wait at the depolarized quiescent state, which can also be associated with the tonic-spiking state in the case of the “spike-less” gFN-burster. The initial state of the gFN-neuron depends on the value of the applied current stimulus $Iapp$ in Eq. (1), which horizontally shifts the position of the fast $V$-nullcline given by $V\u02d9=0$, relative to that of the slow $h$-nullcline given by $h\u02d9=0$; see Figs. 2 and 6. In this study, we examine how these mechanisms influence the dynamical behavior, rhythmic outputs, and bifurcations in the 3-cell motifs with several distinct circuitries presented in Fig. 2(a) and in Fig. 5, using the methodology outlined in Sec. II.

A weakly coupled network obeys the release mechanism when the slow $h$-nullcline lies close enough to the lower knee of the fast cubic $V$-nullcline in the phase space of an intrinsic HH burster in Fig. 1 or in the ($h,V$)-plane of the gFN-neuron in Fig. 6(a). As the value of the $Iapp$-parameter decreases, or alternatively, if a cell receives an inhibitory synaptic current from its pre-synaptic cell(s), its $V$-nullcline shifts horizontally leftward in the ($h,V$)-plane. In the release case, this can lead to the occurrence of tangency of both the nullclines near the lower knee of the $V$-nullcline. Further strengthening of inhibition causes the nullclines to cross locally twice, giving rise to a new stable equilibrium state on the lower hyperpolarized branch on the $V$-nullcline, which emerges through a saddle-node bifurcation. This ceases oscillations in the cell, as it remains effectively hard-locked in the hyperpolarized state. Otherwise, if the inhibition is not strong enough to cause the saddle-node bifurcation, it merely narrows the gap between the nullclines, which makes the cell slow down due to the bottleneck effect preceding the saddle-node bifurcation. After the presynaptic cell(s) traverses its active on-state and becomes inactive on its lower hyperpolarized branch of the fast $V$-nullcline below the synaptic threshold, the post-synaptic cell is released from inhibition and its $V$-nullcline shifts rightward back to its original position. This eliminates the stable hyperpolarized equilibrium state; therefore, the given cell can start a new oscillatory cycle and during its active phase can in turn play the opposite role of inhibiting its post-synaptic cells in the network.

The escape mechanism requires that (i) one of the unperturbed gFN-cells, say the red cell 3, is initially locked at the stable depolarized equilibrium state at the transverse intersection of the slow $h$-nullcline and the fast $V$-nullcline near its upper right knee, as shown Fig. 6(b). The other condition (ii) is that meanwhile the other cell(s) must transition along the perturbed hyperpolarized branch of the $V$-nullcline until either the post-synaptic cell, say the green cell 2, reaches the lower knee from where it jumps up, switching from the “off” to “on” states. As soon as its voltage goes over the synaptic threshold, cell 2 results in a fast inhibitory current that shifts the $V$-nullcline of cell 3 leftward, away from the $h$-nullcline, and eliminates the stable equilibrium state so that cell 3 jumps down onto the inhibited hyperpolarized branch of the $V$-nullcline and so forth.

Below, we will show how inhibitory 3-cell networks can employ the two generic mechanisms to synergistically orchestrate the shifts in the positions of the nullclines of post-synaptic neurons so that the constituent cells can alternate cyclically between their active and inactive states, resulting in several stable rhythmic outcomes. We reiterate that we only consider weakly coupled networks here to ensure the visual continuity of the Poincaré return maps for the phase lags. Our choice of the sigmoidal shape for the slow $h$-nullcline ensures that the system can exploit the bottleneck effect of saddle-node bifurcations to produce a variety of rhythmic outcomes. Increasing coupling strength causes faster non-smooth convergence and hard-locking to stable FPs in the maps, corresponding to phase-locked rhythms; see Fig. 8 illustrating the effect of increasing inhibition in the symmetric motif.

In Secs. III A–III H, using Poincaré return maps and bifurcation diagrams, we will demonstrate the onset of various rhythms and their transitions in the fully symmetrical system, as we vary control parameters that can also be manipulated in neurophysiological experiments. We will then show the behavioral ranges and transitions occurring in various asymmetrical network configurations, otherwise impossible in the fully symmetric system. This is followed by a brief discussion of the stable synchronous state with all three cells oscillating together, as well as a special network configuration without any phase-locked rhythms.

### A. Symmetric motif

Figure 4(e) represents the so-called bi-parametric sweep, or the bifurcation diagram, of the fully symmetric 3-cell motif depicted in Fig. 2(a) as two parameters: the coupling strength $g$ of all six inhibitory synapses and the applied current $Iapp$ are varied. One can see that the parameter space has three color-coded regions where the network produces either three stable pacemakers (PM), two traveling waves (TW), or five coexisting rhythms: 3 PMs and 2 TWs.

The corresponding 2D Poincaré maps for the phase lags, depicted in Figs. 4(A1)–4(A4), demonstrate the transitions and bifurcations in the network due to the escape mechanism as the inhibitory coupling is increased. The initial gap between the slow $h$-nullcline and the upper knee of the fast $V$-nullcline is small enough so that the modulating bottleneck effect makes either cell linger longer in the active on-phase near the upper knee, while the other two cells transition through the inactive off-phase, thereby promoting pacemaker rhythms corresponding to the three stable FPs: blue at $(0.5,0.5)$, green at $(0.5,0)$, and red $(0,0.5)$. As the synaptic coupling $g$ is increased in strength, the gap between the nullclines of the post-synaptic cells widens [see Fig. 6(b)], thus weakening the bottleneck effect so that the circular motion on the limit cycles in the phase plane becomes more uniform. As we increase $g$ further, the unstable TWs at $(0.33,0.67)$ and $(0.67,0.33)$ become stable through a secondary Andronov–Hopf or torus bifurcation [Fig. 4(a2)], with the TW attraction basins gradually increasing, while those of the PMs diminishing in size [Fig. 4(a3)]. Note that there are an even number of saddle FPs (labelled by gray $\u22c4$s) in the maps: their (separatrices’) role is to separate the attraction basins of the coexisting stable FPs. With a further increase in $g$, a pair of nearby saddles approach each stable PM fixed point and merge with it through a pitch-fork bifurcation. The three new saddles now equally partition the attraction basins of the two remaining TWs in Fig. 4(a4).

The maps depicted in Figs. 4(a1)–4(d1) illustrate the bifurcation stages as the symmetric 3-cell motif transitions from the escape to the release mechanism (Fig. 6). This occurs as $Iapp$ is decreased along the vertical dotted line in the bifurcation diagram in Fig. 4(e). Decreasing $Iapp$ shifts the fast $V$-nullcline leftward and, therefore, gradually increases the gap between its upper knee and the slow $h$-nullcline, while simultaneously decreasing the gap near the lower knee of the $V$-nullcline. As seen in the corresponding return maps in Fig. 4, the network transitions from stage A1: a motif producing only three stable PMs; to stage B1: a motif producing only two stable TWs, after the pitch-fork and torus bifurcations; to stage C1: a motif with co-existing TWs and PMs (as PMs re-emerge after reverse pitch-fork bifurcations). Finally, reverse torus bifurcations make the TWs unstable again and restore the motif with only three stable PMs at stage D1. This is due to the dominating bottleneck effect near the lower knee, which substantially amplifies the slow–fast separation in the rate of circulation along the stable limit cycles in the ($h,V$)-plane.

Figure 7 shows the transitions occurring based on the release mechanism, as the synaptic inhibition is increased in the 3-cell network. The initial gap between the slow $h$-nullcline and the lower knee of the fast $V$-nullcline is chosen so that both stable PM and TW rhythms can coexist for the given coupling strength. As the synaptic inhibition $g$ is gradually increased, this gap narrows, thereby slowing the progression of the post-synaptic cells near the lower fold (dwelling time here is inversely proportional to the square root of the size of the gap between the nullclines). This makes the pacemaker activity dominant in the network. Further increase in $g$ can lead to hard-locking, as the nullclines intersect at the lower branch, after crossing the knee. This results in the gradual shrinking of the attraction basins of TWs and a corresponding increase in those of PMs, as seen in Figs. 7(a$\u2032$)–7(c$\u2032$). TWs finally become unstable through a secondary Andronov–Hopf/torus bifurcation in Fig. 7(d$\u2032$) and hence not observable in the network anymore.

Let us now consider another bifurcation pathway in the bifurcation diagram in Fig. 7(e), for a lower value of $Iapp=0.3956$, and, therefore, displaying a more pronounced release mechanism in the network. The return maps with increasing synaptic inhibition are presented in Fig. 8. As $g$ is gradually increased, the network bifurcates from the dominant PM case in Figs. 8(a$\u2032\u2032$) and 8(b$\u2032\u2032$) to a configuration that supports both PMs and TWs in Fig. 8(c$\u2032\u2032$), following torus bifurcations. A further increase in the value of $g$ leads to hard-locking in Fig. 8(d$\u2032\u2032$), due to the tangency or the intersection of the nullclines near the lower knee, for the temporarily inhibited post-synaptic cells. This can be seen in the jagged phase-lag trajectories in the return map, which quickly converge to the pacemaker fixed points.

Finally, we can conclude that increasing coupling strength leads to pronounced pacemaker behaviors in the networks featuring the release mechanism. In contrary, traveling waves become more dominant at stronger coupling values in networks featuring the escape mechanism. This means that symmetry *per se* is insufficient to predict *a priori* what rhythmic outcomes a given network can produce, without knowing the qualitative mechanisms of rhythmogenesis, escape, or release and the quantitative strength of synaptic connectivity. This knowledge is vital to make testable predictions regarding possible dynamics in various biological systems of coupled oscillators and neurophysiological experiments with CPG circuits.

### B. Mono-biased motif

We will now investigate the effect of a single asymmetric connection within an otherwise fully symmetric, weakly coupled system [see Fig. 5(a)]. While we focus on asymmetric increase or decrease in the strength of a single connection ($g31$), the results can be extended to any of the other connections by symmetry and interpreted accordingly without loss of generality. As seen previously with the symmetric motif, the rhythmic behaviors of this network also vary depending upon the gap between the nullclines, the transitions between soft and hard locking, as well as the release and escape mechanisms of the cells. Here and in subsequent discussions of other asymmetric motifs, we will show a bifurcation diagram varying the parameters $g$ (or $g31$ in this case) and $Iapp$, which effectively depicts a range of rhythmic behaviors exhibited by the network. We then pick a few parameter values to demonstrate detailed phase lag return maps and their transitions for cells obeying the release mechanism, as the synaptic strength is gradually increased (Fig. 10). This is then followed by a similar demonstration for cells obeying the escape mechanism. For these asymmetric network configurations, we observe asymmetric bifurcations in which only one or two pacemakers may appear or disappear, rather than all pacemakers (or all traveling waves) simultaneously, as seen in the fully symmetric network for any given parameters.

Figure 9 represents the bifurcation diagram for the mono-biased motif [shown in Fig. 5(a)], as the parameters $Iapp$ and $g31$ are varied, while the remaining synaptic connections are held constant at $gall=0.001$. The network can produce several distinct multistable behaviors composed of just pacemakers (PM), just traveling waves (TW), just phase slipping (PS), or a combination of pacemakers with phase slipping (PM/PS) or traveling waves (PM/TW). Transitions between these regions are due to saddle-node (SN) bifurcations eliminating or restoring FPs to the map. The points A–F highlighted along the dotted line near the top of the bifurcation diagram indicate the parameter values used for the phase lag return maps for the cells obeying the escape mechanism (Fig. 11), while the points A$\u2032$–F$\u2032$ on the dashed line near the bottom highlight the parameter values for the cells obeying the release mechanism (Fig. 10). The vertical dotted line represents the parameter values where the network retains full symmetry, with $g31=gall$. As such, the rhythmic behaviors and the transitions as we move along this line are identical with the fully symmetric motif.

For cells obeying the release mechanism shown in Fig. 10, we initially disable the synapse $g31=0$, while all other synaptic connection strengths are held constant at $gall=0.001$. In this case, we observe that the network is dominated by the green PM rhythm with the largest basin of attraction, while the blue PM rhythm is also stable although with a smaller basin [Fig. 10(a)]. One may observe the presence of two saddle nodes (gray diamonds) and an interesting pattern of whorls in the phase space near the original location of the purple TW (0.66, 0.33), which is currently not seen in the network. Restoring the missing synapse $g31=0.00081$ in Fig. 10(b$\u2032$) increases the blue PM basin and also leads to the formation of the purple TW pattern through a saddle-node bifurcation, which gives rise to a third saddle around the purple TW FP. Strengthening of this synapse at $g31=0.00108$ leads to the appearance of another saddle node and a fixed point corresponding to the red pacemaker rhythm, both of which rapidly diverge [see Fig. 10(c$\u2032$)]. With further increases in the strength of $g31$ through $0.00135,0.004,and0.008$ in Figs. 10(d$\u2032$)–10(f$\u2032$), the purple TW becomes the dominant rhythm of the network via a series of saddle-node bifurcations where the red, green, and blue PM FPs, respectively, merge with the three saddles surrounding the purple TW and disappear one after the other.

Figure 11 shows the evolution of Poincaré return maps for cells obeying the escape mechanism as $g31$ gradually increases. In Fig. 11(a) with $g31=0$, the network is dominated by a single clockwise TW (black FP). Figure 12 shows four identical panels [the same as Fig. 11(a)] stitched together to aid the visual inspection of the trajectories and their convergence. A saddle (gray diamond) exists very close to the black TW FP in the phase space, which gives rise to interesting dynamics such that two trajectories (red lines) starting from very close initial conditions (on either side of the blue line) take entirely different paths but ultimately converge to the same fixed point (black TW). Figure 13 shows a similar alignment of four identical return maps for this network, but at a slightly higher value of $Iapp=0.5875$, chosen within the phase slipping (PS) region of the bifurcation diagram in Fig. 9, close to the parameter values corresponding to the map depicted in Fig. 11(a). Comparing this with Fig. 12 shows that the black TW FP and the saddle merge and disappear following a saddle-node bifurcation and give rise to a stable PS pattern in Fig. 13. The network has no phase-locked rhythms in this state, and all the trajectories from different initial conditions converge onto the red stable PS invariant circle, after an initial transient. Conversely, the blue line represents an unstable PS pattern. If a trajectory starts with initial conditions exactly on this blue line, it will continue to move along this unstable PS pattern, but slight perturbations would lead to diverging paths that then ultimately converge onto the stable red PS pattern. Note that within these PS patterns, cells 2 and 3 remain nearly phase-locked, while cell 1 undergoes phase-slipping and so runs at a different frequency compared to the other two cells. Due to the existence of such subpopulations of cells that run at distinct frequencies, PS has also been referred to as a chimera state. Stable and unstable PS patterns are elaborated further in Secs. III C–III H.

As the coupling strength $g31$ increases from $0$ to $0.00027$ in Fig. 11(b), the unstable PS pattern first undergoes a saddle-node bifurcation, giving rise to an unstable fixed point (white dot) and a saddle. The saddle then undergoes pitch-fork bifurcation to give rise to the blue pacemaker rhythm, as well as two additional saddles. Further increases in $g31$ give rise to the red and green pacemakers, following their respective saddle-node bifurcations in Figs. 11(c) and 11(d). The basin of attraction of the black TW continues to diminish and undergoes an additional torus bifurcation in Fig. 11(e), creating an invariant circle, while the green PM rhythm is lost following a saddle-node bifurcation, resulting in a larger red PM basin. In Fig. 11(f), the invariant circle and the black TW rhythm are finally lost as the black repellor converges with a saddle, and the system is dominated by just the red and blue PMs, with very rapid convergence for many initial conditions (exemplified by large white regions on the return map, where convergence is so rapid that traces are not apparent in these areas).

### C. Double-biased motif

In this section, we study the dynamics of the double-biased motif [Fig. 5(b)], where asymmetry in the network is achieved by simultaneously altering a pair of connections, $g31$ and $g13$ between cells 1 and 3. Figure 14 shows the bifurcation diagram for this motif, with the sampled parameter values for the release (Fig. 15) and the escape (Fig. 16) mechanisms highlighted along the bottom dashed and top dotted lines, respectively. The vertical dotted line represents symmetric network configuration as described previously. The bifurcation diagram reveals that the network can produce just pacemakers (PM), just traveling waves (TW), just phase slipping (PS), or any combination of a pair of such rhythms (PM/PS, PM/TW, TW/PS). Pacemaker behavior dominates at weak coupling, while phase slipping does so at strong coupling. Other rhythms exist primarily near the mid-ranges of values for $Iapp$ or close to full symmetry for the values of $g$.

Figure 15 shows the evolution of the Poincaré return maps for the release mechanism as the synaptic strengths of both $g31$ and $g13$ increase from $0.0$ to $0.0045$, while all other synaptic strengths remain constant at $0.001$. When both the synapses are absent or very weak ($g31=g13=0.0005$) in Fig. 15(a$\u2032$), the network produces a single stable rhythm of the green PM. This could be inferred from the fact that only cell 2 (green) has outgoing inhibitory connections to the other two cells. As we strengthen the coupling ($g31=g13=0.001$) in Fig. 15(b$\u2032$), following a series of saddle-node and pitch-fork bifurcations, the blue and red PMs emerge along with multiple saddles and two repelling FPs corresponding to unstable TWs. Further increase in the synaptic strength ($g31=g13=0.0012$) in Fig. 15(c$\u2032$) causes the unstable TW FPs to disappear through saddle-node bifurcations and the corresponding increase in the blue and red PM basins. In Fig. 15(d$\u2032$) with $g31=g13=0.0015$, the blue and red PMs disappear and give rise to a stable invariant circle (gray) through a heteroclinic saddle-node bifurcation. This invariant circle corresponds to a PS rhythmic pattern that wraps around the torus and coexists along with the green PM. Hand-drawn lines in black are sampled to illustrate the attraction basins bounded by the incoming separatrices of the saddles.

Figure 16 illustrates the dynamics and the corresponding return maps for the escape mechanism as the synaptic strengths $g31$ and $g13$ are gradually increased through (0, 0.000 676, 0.001, 0.001 28, 0.001 55, 0.003 31) at the points from A through F. When both the synapses are turned off in Fig. 16(a), the network produces a single stable rhythm with the green PM, as described previously. Figure 17 shows four identical panels [the same as in Fig. 16(a), with greater detail] placed next to each other to aid visual inspection of the trajectories and their convergence. Sample trajectories are shown in black, along with the green pacemaker, two saddles, and an unstable fixed point. In Fig. 16(b), the two saddles undergo pitch-fork bifurcations to give rise to the blue and red PMs, as well additional saddles. As the synapses are further strengthened in Figs. 16(c) and 16(d), the two TWs emerge and then disappear through saddle-node bifurcations. Further increases lead gradually to the disappearance of the green pacemaker in Fig. 16(e) via a saddle-node bifurcation and of the red and blue PMs in Fig. 16(f) through a heteroclinic saddle-node bifurcation that gives rise to a single invariant circle of the gray phase-slipping pattern that wraps around the torus.

### D. Driver-biased motif

Another type of asymmetry we investigate is the driver-biased motif, where the two outgoing synapses from cell 3 ($g31=g32$) are manipulated, while the remaining connection strengths are held constant. Figure 18 shows the bifurcation diagram and the sampled values in it, corresponding to the release and escape mechanisms, as described previously. As can be expected from this asymmetry, for sufficiently strong synaptic coupling for the outgoing connections from cell 3 (red) that acts as the driver, the network is dominated by the red PM rhythm. For weaker coupling strengths, one can observe TWs, PS, and a combination of PMs with TWs or PS. At weaker coupling, one may also see the blue and green pacemaker rhythms.

Figure 19 shows the return maps for the cells obeying the release mechanism. For small values of $g31=g32$ in Fig. 19(a$\u2032$), the blue and green PMs coexist. As the synaptic strengths are increased in Fig. 19(b$\u2032$), the red PM and the two TWs emerge following a series of saddle-node bifurcations. Further increases in the synaptic strengths lead to the disappearance of the original blue and green PMs in Fig. 19(c$\u2032$), as well as the disappearance of the two TWs in Fig. 19(d$\u2032$) through saddle-node bifurcations, resulting in a single red PM rhythm of the circuit, under the control of the biased-driver cell.

Following the return maps for the escape mechanism in Figs. 20(a)–20(f) as the values of $g31$ and $g32$ are simultaneously increased shows that the network is initially dominated by a single PS rhythm in Fig. 20(a). A red PM rhythm then emerges following a saddle-node bifurcation in Fig. 20(b). A series of saddle-node bifurcations then gives rise to the blue and green PMs, as well as the two TWs in Fig. 20(c). With further increase in synaptic strengths, the blue and green PMs first disappear in Fig. 20(d), followed by the loss of the two TWs in Fig. 20(e), through saddle-node bifurcations, resulting in a phase slipping pattern that coexists with the red PM. Finally, for very strong coupling in Fig. 20(f), the phase slipping pattern also disappears, giving rise to a single red PM rhythm for the circuit, driven by the dominant cell.

### E. Clockwise-biased motif

The next asymmetric network configuration examined is the clockwise-biased motif, in which all the clockwise connections: $g12$, $g23$, and $g31$ are manipulated simultaneously, while the anti-clockwise connections are held constant. Figure 21 shows the bifurcation diagram for this motif, which reveals that the TW rhythms dominate this network. A single TW is seen at either end of the coupling strength spectrum, as can be expected from such asymmetry, while both the TWs are seen in between these two regions. PMs and PM/TW combinations are also seen in parametric regions close to the fully symmetric network configuration, for both the escape and release mechanisms.

Figure 22 shows the return maps at parameter values sampled for the release mechanism, as the clockwise synapses are gradually strengthened. Initially, the network is dominated by a single clockwise TW (black) in Fig. 22(a$\u2032$). Following a series of saddle-node bifurcations, the three PMs emerge in Fig. 22(b$\u2032$). This is followed by the emergence of the anti-clockwise TW (purple) in Fig. 22(c$\u2032$) and the disappearance of the clockwise TW (black) in Fig. 22(d$\u2032$), through their respective torus bifurcations. The three PMs then disappear via saddle-node bifurcations in Fig. 22(e$\u2032$), leading to a single anti-clockwise TW rhythm dominating the network. Further strengthening of the synapses in Fig. 22(f$\u2032$) leads to hard locking, with a single TW rhythm.

Figure 23 depicts the dynamical transitions for the escape mechanism and looks almost identical to that of the release mechanism shown in Fig. 22. The transition from the black clockwise TW in Fig. 23(a) to the purple counter-clockwise TW in Fig. 23(f) occurs via the formation of the three PMs, the disappearance of the black TW, the appearance of the purple TW, followed by the disappearance of the three PMs. The escape mechanism is more conducive to permitting the transitory limit cycle behavior following torus bifurcations, and therefore, these can be observed in finer detail here in Fig. 23(c).

### F. Emergence of a transitive torus

In this section, we describe a route to the emergence of a transitive torus without any fixed points or invariant circles in the return map. For the asymmetric network under consideration (Fig. 24, top), this would correspond to a lack of any phase-locked or periodically varying rhythmic outcomes. Such a transitive torus has what is called the *everywhere dense covering*, and a trajectory starting from any initial condition fills in the entire phase portrait over time.

Figure 24 (top) shows the asymmetric network, where cells 1 and 3 have weak inhibitory coupling ($g13=g31=0.0003$), cells 2 and 3 have a slightly stronger coupling ($g23=g32=0.0005$), while cells 1 and 2 have stronger asymmetric coupling ($g12=0.001$, and $g21=0.0051$). Unlike the previous examples, we maintain the synaptic strengths constant in this case, while varying the external current drive $Iapp$ for all the cells simultaneously through 0.4, 0.46, 0.5, 0.572 982, 0.594, and 0.61 from panel A through panel F in Fig. 24, transforming individual cells gradually from the release to the escape mechanisms. In order to effectively demonstrate both the stable and unstable fixed points, invariant circles, as well as their transitions, we compute both the forward (gray) and backward (red) trajectories and plot them in the Poincaré return maps.

In Fig. 24(a), the network produces a single stable invariant circle or a phase slipping pattern, shown in gray. In addition, the red trajectories show the presence of an unstable fixed point as well as a saddle, both of which disappear after undergoing a saddle-node bifurcation in Fig. 24(b) to give rise to an unstable PS pattern, which coexists with the stable PS. Also, note that the stable and unstable PS patterns run in opposite orientations. In Fig. 24(c), the stable PS pattern is lost due to a saddle-node bifurcation, giving rise to a stable fixed point and a saddle. The network therefore has a single stable rhythm (the black TW rhythm, shown in previous examples) that coexists with an unstable PS pattern. The fixed point and the saddle then gradually move further apart from each other, wrap around the torus, and then merge and disappear through a saddle-node bifurcation in Fig. 24(d), again giving rise to a stable PS pattern, which runs along the same orientation as the unstable PS pattern. The PS patterns then start disappearing via torus breakdown^{49} as the stable and unstable invariant circles start merging in Fig. 24(e) and finally give rise to the transitive torus in Fig. 24(f). Only a single long trajectory is plotted in Fig. 24(f), and it retraces the entire phase space over time, without any fixed points or invariant circles.

### G. Invisible heteroclinic bifurcations

Here, we illustrate the “invisible” role of some heteroclinic bifurcations of the saddles and how they determine the structure of the 2D Poincaré return maps and specifically how they shape the attraction basins of the co-existing stable FPs. Let us examine the transformation of the attraction basins depicted in Fig. 25. In both the cases, the map has two stable FPs, green and blue, corresponding to the pacemaker rhythms with phase lags $\u223c(0.5,0)$ and $\u223c(0.5,0.5)$, respectively. In addition to the persistent repeller at the origin, there are three more saddle FPs (labelled by $\u22c4$’s) so that the total number of hyperbolic FPs on the torus is even. Of special interest here is the saddle to the left of the stable (blue) FP around $(0.5,0.5)$. More specifically, let us follow its left outgoing (unstable) separatrix (set) to find its destination, or the $\omega $-limit set, as the number of iterates increases. As the map on the torus is defined on mod 1, the separatrix disappears when it reaches the left wall given by $\Delta 12=0$ and comes back into the map from the right wall given $\Delta 12=1$. Next, it slides above the incoming separatrix of the other saddle (to the right of the blue FP) to converge to the blue FP—its $\omega $-limit set [Fig. 25(a)]. As $g31$ is slightly increased from 0.0071 to 0.0072, this separatrix first merges with the incoming separatrix of the other saddle (right) to form a one-way heteroclinic connection, which is followed by its shift further downwards below the saddle to switch to another $\omega $-limit set—the stable FP around $(0.5,0)$, which corresponds to the green PM [Fig. 25(b)]. This heteroclinic bifurcation drastically repartitions the sizes of the attraction basins of the co-existing FPs: the blue FP is no longer dominant as a majority of initial conditions would now converge to the green FP. This example highlights the pivotal role of homoclinic and heteroclinic bifurcations that underlie major reconstructions of the phase space, while preserving the existence and stability of FPs in these return maps and other systems, in general.

### H. Stable synchronous state

In this section, we discuss the unexpected case of a stable synchronous state, where all the cells oscillate together, which has not been previously observed in 3-cell circuits coupled by the *fast* FTM synapses. We note, however, that the discovery of the stable synchronous state, in addition to several other exponentially stable polyrhythms, was recently reported in 4-cell inhibitory circuits, made up of gFN-model neurons.^{44} So far, we have known that, as soon the fast inhibitory connections in 3-cell motifs are replaced with excitatory ones,^{1} all the attractors of the map become repellers, the incoming (stable) separatrices (sets) of the saddles become outgoing (unstable) and vice versa. This is also true for half-center oscillators made up of two reciprocally inhibitory cells. It was shown previously^{50} that the synchronous state in such a HCO made up of two plain FN-oscillators is unstable, in the case of FTM or similar inhibitory coupling. Note that we say a synapse is fast when the corresponding current is only slightly delayed compared to the timing of the spike or the burst that initiates it and decays quickly after the voltage of the pre-synaptic cell lowers below the synaptic threshold. A synapse is slow when the decaying current generated by the pre-synaptic cell lasts much longer; therefore, its duration can be compared with the inter-spike/burst period and not with the spike/burst duration as in the case of fast synapse. This property of long decay is key to understanding how two neurons coupled reciprocally with slow inhibitory synapses can suppress in-phase oscillations.^{51} Basically, if both are given a “window of opportunity” to start together, they will continue oscillating in phase. Otherwise, if the initial conditions are different and the neurons do not start within this window, either one may surpass the other if the coupling is strong enough and the decay is long enough. As such, there is no asymptotic convergence to the synchronous state but to the anti-phase rhythm. This is not the case when one examines HCOs made of endogenous bursters with weak inhibitory coupling, using the fast FTM-synapses, provided that the level of the synaptic threshold goes through the fast spikes within bursts; if the level is below the spikes, the burster effectively becomes a FN-neuron or a gFN-neuron. It was demonstrated in Refs. 52 and 53 that in-phase synchrony of two inhibitory coupled HH-like bursters can be stable and asymptotic for three different models of the fast inhibitory synapses. Moreover, such busters can also converge to other close synchrony-like states, with one or several spikes apart, due to the spike timing and interactions, which make inhibition act like excitation; the reader can find further details in the above references. Increasing the coupling strength breaks down the synchrony arising from weak spike interactions; therefore, the neurons start bursting in alternation.

One can see from Fig. 26 that the FP at the origin in the depicted 2D maps is no longer a repeller, unlike all the previous cases, where it corresponds to an unstable synchronous rhythm with $\Delta 12=\Delta 13=0$. For the given parameters values, the origin becomes an *asymptotic* attractor with a relatively large basin (shown in yellow), to which the nearby initial conditions converge. Depending on the gap between the nullclines, the fast $V\u2032=0$ and the slow $h\u2032=0$, this newly formed attractor co-exists with the stable PMs (green, red, and blue) in Fig. 26(a$1$) or with stable PMs and TWs (black and purple spirals) in Fig. 26(a$2$).

This figure also introduces another concept incorporated into MotifToolbox^{43} to reveal the attraction basins of the co-existing stable FPs (labelled by color-matching dots) in greater details, as presented in Fig. 26(b$1$) and 26(b$2$). Using these diagrams, we can identify the locations of the repellers in the $(\Delta 12,\Delta 13)$-plane at the junctions of three distinctively colored regions. For example, in Figs. 26(a$1$) and 26(b$1$), one can spot such repellers at the locations of the TW FPs at $(13,23)$ and $(23,13)$, as well as six more repellers (white dots) and six new saddles (gray diamonds). The number of repellers in Figs. 26(a$2$) and 26(b$2$) provides an explanation as to how the synchronous FP at the origin becomes stable. Due to its location and symmetry on the 2D torus, it undergoes a degenerate pitch-fork bifurcation simultaneously along the lines $\Delta 12=0$, $\Delta 13=0$, and $\Delta 12=\Delta 13$ and becomes stable in all three directions. This bifurcation gives rise to six new saddles and six new repellers nearby. Given that the return maps are defined on a phase torus in modulo 1, the particular visualization approach stitching together four identical panels can be employed to magnify the vicinity of the origin, as shown in Fig. 26(d). This panel reveals how the boundary of the attraction basin of the synchronous state at the origin is geometrically determined by the stable and the unstable separatrices of the saddles, which emerges through this bifurcation sequence.

## IV. BOTTOM UP APPROACH: A 2*θ*-BURSTER FOR 3-CELL MOTIFS

The $2\theta $-neuron model is motivated by the dynamics of endogenous bursters, with two characteristic slow phases: tonic-spiking and quiescent. Similar to the notion of the relaxation oscillator vs the FN-neuron, the equation describing the so-called “spiking” $\theta $-neuron^{54} in the context of neuroscience was known for a long time since it was introduced in the classical mathematical theory of synchronization.^{55} Its core is a homoclinic saddle-node bifurcation on a torus or on an invariant circle, or the SNIC bifurcation, which occurs on the V-shaped boundaries of synchronization zones, also known as Arnold tongues, in the parameter plane; see Fig. 30. The $\theta $-neuron capitalizes on the pivotal property of the saddle-node bifurcation—the phantom bottleneck effect that gives rise to slow and fast time scales in the dynamics of systems ranging from simple 1D to higher-order models. In the gFN-model, the saddle-node bifurcation occurs at the quadratic tangency of the nullclines, $V\u2032=0$ and $h\u2032=0$, in the phase plane (Fig. 6). Figure 1 illustrates the same principle in the 3D phase space of the reduced leech heart interneuron, where the quiescent phase of bursting can be controlled by varying the gap between the slow nullcline and the right hyper-polarized knee. Recall that a similar saddle-node bifurcation in this model, which controls the tonic spiking phase and the number of spikes per burst, is associated with the famous blue-sky catastrophe.^{32,56–58}

A key feature of the $2\theta $-neuron is the occurrence of two saddle-node bifurcations, which introduce two slow transitions into its dynamics, with two fast switches in between. Similar to endogenous bursters with two slow transient states—the active tonic-spiking and the quiescent phases that can be controlled independently, we can manage the durations of the two analogous states in the $2\theta $-neuron: “on” at $\pi $ and “off” at $0$, using the same bottleneck post-effects of the two saddle-node bifurcations. This allows us to regulate its duty cycle, which is the fraction of the active-state duration compared

to the burst period; see Fig. 27. As seen in this figure, the $\theta $-model phenomenologically depicts spiking cells, while the $2\theta $-neuron can be treated as a “spike-less” burster, similar to the gFN-neuron discussed previously. Below, we demonstrate that the network dynamics produced by a 3-cell motif composed of inhibitory $2\theta $-bursters preserve all the key features seen in a motif composed of the three gFN-neurons as well.

A 3-cell motif comprising $2\theta $-neurons, coupled with fast inhibitory FTM-synapses, is given by the following system:

One can observe that the phase dynamics of an individual $2\theta $-neuron are governed by the terms $\omega \u2212cos\u2061(2\theta )$. As long as the frequency $0<\omega \u22641$, there are two stable and two unstable equilibria: at the bottom around $\theta \u22430$ and at the top near $\theta \u2243\pi $; they are associated with the hyperpolarized and the depolarized quiescent states of the neuron. When $\omega >1$, the $2\theta $-neuron becomes oscillatory through two simultaneous saddle-node bifurcations (SNIC) on a unit circle $S1$, where $\theta $ is defined on modulo one. Moreover, whenever $\omega =1+\epsilon $, where $0<\epsilon \u226a1$, this new “burster” possesses two slow phases: the active “on” state near $\theta =\pi $ and the inactive “off” state near $0$ on $S1$, alternating with fast counter-clockwise transitions that are referred to as the upstroke and the downstroke, respectively. For greater values of $\omega $, the active and inactive phases should be defined by $\pi /2<\theta \u22643\pi /2$ and $3\pi /2<\theta \u2264\pi /2$, respectively. The latter phase is below the synaptic threshold, which is set by $\theta th=\pi /2$ so that $cos\u2061(\theta th)=0$, thus equally dividing the unit circle. The duty cycle of the $2\theta $-neuron is controlled by the term $\alpha cos\u2061\theta $, provided that it remains oscillatory as long as $\omega \u2212|\alpha |>1$. Note that when $\alpha =0$, the duty cycle is 50% and the oscillations are even. The active or inactive phases can be extended or shortened, respectively, with $\alpha <0$, making the duty cycle greater or vice versa—the duty cycle of individual neurons can be decreased with $\alpha >0$.

The $2\theta $-neurons are coupled in the 3-cell network using the fast inhibitory FTM^{42} synapses. The “sigmoidal” term $11+ekcos\u2061\theta i$, ranging between $1$ and $0$, rapidly triggers (here, $k=10$) an influx of inhibition flowing from the $i$th pre-synaptic neuron into the $j$th post-synaptic neuron, as soon as the former enters the active on-phase above the synaptic threshold $cos\u2061(\theta th)=0$, i.e., $\pi /2<\theta i<3\pi /2$. The strength of the inhibitory coupling is determined by the maximal conductance $\beta ij$ that slows down the rate of increase of $\theta j\u2032$ in the $j$th post-inhibitory neuron because of the negative sign of the coupling term. To translate the synaptic input into qualitative inhibition, the sign of the input is switched upon crossing the values $\theta =0$ and $\theta =\pi $. This is achieved by multiplying all the coupling terms of each ODE by $[1\u221221+eksin\u2061\theta ]$. When $0<\theta <\pi $, the inhibition is a negative input to slow the transition into bursting. When $\pi <\theta i<2\pi $, the inhibition is a positive input making the transition out of bursting toward quiescence faster. This is logically consistent as, in general, inhibition should shorten the duration taken for the post-synaptic neuron to leave its active phase. Optionally, one can replace this term with $[1\u221211+eksin\u2061\theta ]$, which breaks the symmetry well.

Figure 28 shows the snapshots of the phase progressions of the three calls on the unit circle $S1$ and depicts how phase lags between the three $2\theta $-neurons are introduced (here, the reference cell is cell 1, in blue), just like in the case of the 3-cell gFN motif in Fig. 3. One can see from Fig. 28(b) that the active green neuron in the active phase near $\theta =\pi $, above the synaptic threshold, inhibits and pushes the other two closer to each other, near the bottom quiescent state at $\theta =0$, by accelerating the red neuron on the downstroke and by slowing down the blue neuron on the upstroke.

Following the same approach used in the weakly coupled HH and gFN models, we use a uniform distribution of initial phase conditions, and hence the phase lags between the three $2\theta $-neurons, and determine the phase locked states that they can converge to, with increasing number of cycles. This approach is illustrated in Fig. 29(a) (compare with Fig. 3) for the symmetric 3-cell motif composed of identical $2\theta $-neurons and equal inhibitory synapses. The corresponding 2D Poincaré return map, with the co-existing fixed points and saddles, is shown in Fig. 29(d), defined on mod 1. By stitching together the opposite sides of this flat map, we wrap it around a 2D torus shown in Fig. 29(b), color-coded accordingly.

One can see that the 2D return maps for the gFN-neurons and $2\theta $-neurons are nearly identical. This implies that our descriptions and modeling approaches to reveal the intrinsic properties of individual neurons and their phenomenological interactions at the network-level are generic and universal. We would like to underline that the proposed concept of $2\theta $-bursters bares a great promise for studies of collective dynamics exhibited by larger and modular networks with a combination of inhibitory, excitatory, and electric synapses, as well as for modeling biologically plausible circuits such as central pattern generators. To conclude this section, we point out another helpful feature of the $2\theta $-neuron paradigm, namely, the ability to find repelling FPs, if any, in the 2D Poincaré map, by reversing the direction of integration of the system (4), i.e., integrating it in the backward time by multiplying the right-hand sides in Eq. (4) by $\u2212$1. Unlike the gFN and other HH-like dissipative neural systems where the backward integration will make solutions run to infinity, it is not the case for $2\theta $-bursters, as the phases on $S1$ will just reverse the direction and spin clockwise on the unit circle.

## V. CONCLUSIONS AND FUTURE DIRECTIONS

The goal of this paper is to *de facto* illustrate that 3-cell and larger neural networks can universally produce the same emergent behaviors in response to parameter variations, provided that the dynamical properties are properly chosen for the synapses and the constituent neurons, whether those are biologically plausible Hodgkin–Huxley type bursters, reduced generalized Fitzhugh–Nagumo neurons, or toy $2\theta $-bursters. In all these cases, we can employ the reduction to the visually evident Poncaré return maps for phase lags, solely derived from multiple voltage traces. This presents a potent computational approach for the thorough analysis of rhythmic behaviors arising in a range of symmetrical and asymmetrical neural networks. By taking advantage of the latest GPU computing paradigms, we perform fast parallel computations of numerous network trajectories and construct the return maps, such as Figs. 7 and 8, within just a few seconds. We demonstrate how the reduced 2D gFN neuron model, given by Eq. (1), in conjunction with these computational techniques, allows for comprehensive examinations of rhythmic behaviors arising in these networks, their underlying mechanisms of release and escape, as well as the construction of detailed bifurcation diagrams [Fig. 4(e)] to study rhythm transitions as the parameters such as the external current drive and the strength of one or more synapses are manipulated. The parameters are carefully chosen so they can also be controlled in neurophysiological experiments with a dynamic clamp to replicate these behaviors in real animal CPGs.^{21,23} Symmetric and asymmetric 3-cell configurations can produce a range of stable and unstable rhythmic behaviors, including phase-locked bursting with pacemakers or traveling waves, as well as the recurrent phase slipping chimeras. A rich set of bifurcations can be induced in these networks including saddle-node, pitch-fork, and secondary Andronov–Hopf/torus bifurcations through parametric changes, resulting in the emergence/disappearance of rhythmic states and the gain/loss of their stability. Finally, we also demonstrate how the 3-cell motif can lose all its stable/unstable fixed points and invariant circles, causing the emergence of a transitive torus, where the network can produce voltage traces whose phase lags vary in a weakly chaotic manner.

We emphasize that 3-cell circuits composed of the gFN-type neurons replicate the multiplicity of rhythmic behaviors and bifurcations seen in more detailed Hodgkin–Huxley type models,^{1} albeit at much lower computational costs. The reduction of the analysis of 3-cell network dynamics to 2D phase lag return maps allows for the simple visual inspection of fixed points, invariant circles, and their bifurcations. A stable FP representing a phase-locked bursting rhythm also remains structurally stable under variations of network parameters, as seen in biparametric scans such as Figs. 9 and 14. These scans show the regions of existence and stability for the FPs, also referred to as synchronization zones,^{55} in the 2D parameter diagrams. The boundaries of such zones correspond to homoclinic and heteroclinic saddle-node bifurcations of fixed points or periodic orbits on an invariant circle; see Ref. 49 and the references therein for more details about tori in neural models. Since there exist several controllable network parameters such as the current drive for each cell, the connection strengths, and other dynamical properties of individual synapses, the nested organization of synchronization zones in the higher dimensional parameter space is depicted in Fig. 30. Such zones have also been known as Arnold tongues^{59,60} for weakly coupled oscillators and other complex vibrating systems. Gradually changing a single parameter can cause a cascade of saddle-node bifurcations as the boundaries of the Arnold tongues are crossed inward/outward, and rhythmic behaviors emerge/disappear as seen in Figs. 22 and 23. Although there exist several controllable network parameters, we demonstrate that many of these manipulations produce qualitatively similar dynamics and transitions, highlighting the effectiveness of these reduction tools.

The extension of these techniques to study larger neural networks with more than 3 cells as demonstrated in Fig. 31 would require additional enhancements, including taking advantage of unsupervised machine learning techniques to analyze the phase-lag return maps in higher dimensions. The *cons and pros* of the visualization approach applied for a 4-cell network are elucidated in Fig. 32. One can observe that in contrast to 2D return maps, the trajectory density undermines the clarity of the representation in this case, and therefore, alternative ways should be used for thorough studies of larger networks. We used a clustering approach based on unsupervised machine-learning for the examination of 4-cell inhibitory neural circuits^{44} and their repertoire of multistable polyrhythms. Certain network topologies display a rich multiplicity (due to permutation-phase shift symmetries^{65}) of multistable [in the motif depicted Fig. 31(a)] or bistable [in the uni-directional motif shown in Fig. 31(b)] polyrhythms, while others can only demonstrate monostability [the specific configuration presented in Fig. 31(c)]. Figures 33(a) and 33(b) present two such stable rhythms produced by 4-cell inhibitory networks: the so-called paired half-center rhythm and the traveling wave. Future work could study the dynamic behaviors arising in modular networks composed of smaller motifs. An example is the 6-cell circuit shown in Fig. 31(d), which is made up of two 3-cell motifs that are also bond together via reciprocally inhibitory synapses. The rhythmic outcomes of the larger network depend on the dynamics of individual motifs as well as the connectivity between them. When each 3-cell motif in Fig. 31(d) is configured to produce pacemaker rhythms only, the cross-inhibition (between cells 1–6, 2–5, and 3–4) suppresses certain combinations of pacemaker patterns in the individual motifs, while promoting others, like the one shown in Figs. 33(c) and 33(d). Using the known principles and rhythmic outcomes of the smaller CPGs, one might simplify the analysis of such modular networks. Our analysis of phase-lags and return maps does not depend on the underlying mathematical equations governing the system. As such, the approach can be generalized to a variety of biological and non-biological complex systems spanning across engineering, economics, population dynamics, dynamic memory, and decision making in animals,^{61} as well as the development of efficient robot locomotion.^{62–64}

## ACKNOWLEDGMENTS

We thank the Brains and Behavior initiative of Georgia State University for the pilot grant support, as well as for the B&B fellowships awarded to J. Collens, D. Alacam, and K. Pusuluri. The authors—former Ph.D. students in the Shilnikov NeurDS lab (Neuro-Dynamical Systems)—thank the past and current lab-mates, specifically J. Schwabedal, J. Wojcik, J. Scully, J. Bourahman, and H. Ju, for inspiring discussions, and also acknowledge S. Kniazev who, using Motiftoolbox,^{43} found that the 3-cell motif could also exhibit a stable synchronous state (depicted in Fig. 26) as a project assignment for *Dynamical Foundation of Neuroscience* course offered by A.S. at GSU. The NeurDS lab is grateful to NVIDIA Corporation for donating the Tesla K40 GPUs that were actively used in this study. This so-called “lab” paper was funded in part by the National Science Foundation (NSF) (Grant No. IOS-1455527).

## DATA AVAILABILITY

The toolkit that supports the findings of this study is openly available as Motiftoolbox in GitHub at https://github.com/jusjusjus/Motiftoolbox, Ref. 43.

## REFERENCES

*Model Neural Networks and Behavior*(Springer, 1985), pp. 37–48.

*Society for Neuroscience Abstracts, 918.04 Neuroscience Meeting Planner, Washington, DC*(Society for Neuroscience, 2011), Vol. 37

*Neural Control of Rhythmic Movements in Vertebrates*(Wiley, 1988).