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θ-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.

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.

FIG. 1.

(a) Snapshot depicting the current states (represented by the blue, green, and red spheres) of three weakly coupled cells at t=0 and their further progressions at t=10, on the bursting orbit (gray) in the 3D phase space of the Hodgkin–Huxley type model of the leech heart interneuron.31,32 The plane Θsyn represents the threshold for the chemical synapses, which divides the active “on” phase (above it) and the inactive “off” phase; here, the active red cell inhibits the quiescent green and blue ones. (b) Burst initiations in successive voltage traces generated by the 3-cell neural network allow us to define the relative delays τ’s and hence the phase lags [given by Eq. (3)] between its constituent bursters; see Refs. 1 and 33 and Sec. II for further details.

FIG. 1.

(a) Snapshot depicting the current states (represented by the blue, green, and red spheres) of three weakly coupled cells at t=0 and their further progressions at t=10, on the bursting orbit (gray) in the 3D phase space of the Hodgkin–Huxley type model of the leech heart interneuron.31,32 The plane Θsyn represents the threshold for the chemical synapses, which divides the active “on” phase (above it) and the inactive “off” phase; here, the active red cell inhibits the quiescent green and blue ones. (b) Burst initiations in successive voltage traces generated by the 3-cell neural network allow us to define the relative delays τ’s and hence the phase lags [given by Eq. (3)] between its constituent bursters; see Refs. 1 and 33 and Sec. II for further details.

Close modal

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θ-neurons to build even simpler multifunctional oscillatory motifs.

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 

Vi˙=ViVi3hi+Iapp+jiGji(Vj,Vi),hi˙=ε[11+ek(ViV0)hi],i,j=1,2,3.
(1)

Here, the state of the ith 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; ε is the reciprocal of some time constant and regulates slow–fast dynamics (0<ε<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˙=0 and h˙=0, respectively [see Fig. 2(b)]. The default values for the parameters are k=10, ε=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 

Gji(Vj,Vi)=gji(VrevVi)Γ(Vj),Γ(Vj)=11+e100(VjVth).
(2)
FIG. 2.

(a) Symmetric 3-cell circuit with inhibitory synapses. (b) The (h,V)-phase portrait of the three coupled cells governed by Eq. (1), depicting two stable periodic orbits: relaxation-like and round-shaped (shown as gray solid and dotted curves), at ε=0.05 and ε=0.3, respectively, which are superimposed with the fast cubic nullcline (dark purple curves—solid unperturbed and dashed in the inhibited case) labeled as V˙=0 and the slow sigmoidal nullcline (orange curve), h˙=0. Blue, green, and red dots on the clockwise periodic orbit represent the time-evolution of the phases of the corresponding cells, 1, 2, and 3, coupled in the network. (c) and (d) Voltage traces generated by the network at ε=0.05 and ε=0.3; see the corresponding limit cycles in panel (b).

FIG. 2.

(a) Symmetric 3-cell circuit with inhibitory synapses. (b) The (h,V)-phase portrait of the three coupled cells governed by Eq. (1), depicting two stable periodic orbits: relaxation-like and round-shaped (shown as gray solid and dotted curves), at ε=0.05 and ε=0.3, respectively, which are superimposed with the fast cubic nullcline (dark purple curves—solid unperturbed and dashed in the inhibited case) labeled as V˙=0 and the slow sigmoidal nullcline (orange curve), h˙=0. Blue, green, and red dots on the clockwise periodic orbit represent the time-evolution of the phases of the corresponding cells, 1, 2, and 3, coupled in the network. (c) and (d) Voltage traces generated by the network at ε=0.05 and ε=0.3; see the corresponding limit cycles in panel (b).

Close modal

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, Γ(Vj)=1, from the inactive “off” state with Γ(Vj)=0, when Vj<Vth, provided that Γ 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=1.5 that makes (VrevVi)<0 defines the inhibitory synapse projected from the active neuron j, with Γ(Vj)=1, to neuron i, which slows down the rate V˙i in Eq. (1). Geometrically [see Fig. 2(b)], the term gji(VrevVi) shifts and skews the fast nullcline V˙i=0 closer to the slow nullcline h˙i=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˙=0) and the slow sigmoidal h-nullcline (h˙=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˙]2+[h˙]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 nth phase lags Δ12(n) and Δ13(n) of cells 2 and 3, respectively, as follows:

Δ12(n)=τ2(n)τ1(n)τ1(n+1)τ1(n),Δ13(n)=τ3(n)τ1(n)τ1(n+1)τ1(n),mod\,1,
(3)

where τin represents the time at which the ith cell reaches the threshold voltage, Vth=0, for the nth time [see Figs. 3(b) and 3(c)]. The sequence of phase lags {Δ12(n),Δ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 Δ12 and Δ13, by initiating multiple trajectories with different initial phase lags (on a grid of size 50×50) 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 Δ12 and Δ13 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(c2), green, and red] and two traveling waves [clockwise—Fig. 3(c1) 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 (Δ12, Δ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.

FIG. 3.

(a1) 2D torus wrapped around a Poincaré return map for the phase-lags between the three cells. (a2) Flattened map on a unit square, revealing five stable fixed points (FPs) (), with their color-coded attraction basins that correspond to the phase-locked states to which the phase lags Δ13(n) and Δ12(n) converge in panels b1 and b2. (c) Time-delays τ21(n) and τ31(n) between the upstrokes in the reference blue cell 1 and in cells 2 (green) and 3 (red), normalized over the network period, defining the phase-lags. These traces exponentially converge to multiple phase-locked rhythms such as the clockwise traveling wave or the pacemaker, corresponding to the black and blue FPs in panel a2.

FIG. 3.

(a1) 2D torus wrapped around a Poincaré return map for the phase-lags between the three cells. (a2) Flattened map on a unit square, revealing five stable fixed points (FPs) (), with their color-coded attraction basins that correspond to the phase-locked states to which the phase lags Δ13(n) and Δ12(n) converge in panels b1 and b2. (c) Time-delays τ21(n) and τ31(n) between the upstrokes in the reference blue cell 1 and in cells 2 (green) and 3 (red), normalized over the network period, defining the phase-lags. These traces exponentially converge to multiple phase-locked rhythms such as the clockwise traveling wave or the pacemaker, corresponding to the black and blue FPs in panel a2.

Close modal

Alternatively, it is also possible that the forward phase-lag trajectories {Δ12(n),Δ13(n)} could converge to attractors other than FPs, as demonstrated in Figs. 15(d), 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.

FIG. 4.

The (g,Iapp)-bifurcation diagram in (e) for the fully symmetric 3-cell motif [in Fig. 2(a)] has several (color-coded) regions where the network can produce only pacemaker (PM) rhythms, or only traveling wave (TW) rhythms, or both PMs/TWs. Panels a1–a4 give the snapshots of the Poincaré return maps for the network due to the escape mechanism as they are sampled at the parameter values (white dots) along the top horizontal line in panel (e). As the coupling strength g increases, the map with 3 stable FPs: blue PM at (0.5, 0.5), green PM at (0.5, 0), and red PM at (0, 0.5) in (A1) and gains two more stable FPs: black TW at (0.33, 0.67) and purple TW (0.67, 0.33) in (A). The growing TW and shrinking PM attraction basins are separated by the separatrices of the saddles (gray ) in (A3). Each stable PM becomes a saddle through a pitch-fork bifurcation at larger g-values, after merging with two nearby saddles in (A4). Panels (a1)–(d1) snapshot several maps as the network transitions from the escape to the release mechanism as Iapp is decreased [along the left vertical dotted lines in (e)]: from 3 stable PMs in (A1), to two stable TWs in (B1), to mixed dynamics with both PM/TWs in (C1), and back to 3 stable PMs in (D1). Parameters for (A1)–(A4): Iapp=0.5886, g=(0.0015,0.006,0.019,0.0225); parameters for (B1)–(D1): g=0.0015, Iapp=(0.493,0.419,0.393). The return maps sampled along the horizontal dashed lines A–D and A-D at the bottom of panel (E) are presented in Figs. 7 and 8 for the networks obeying the release mechanism.

FIG. 4.

The (g,Iapp)-bifurcation diagram in (e) for the fully symmetric 3-cell motif [in Fig. 2(a)] has several (color-coded) regions where the network can produce only pacemaker (PM) rhythms, or only traveling wave (TW) rhythms, or both PMs/TWs. Panels a1–a4 give the snapshots of the Poincaré return maps for the network due to the escape mechanism as they are sampled at the parameter values (white dots) along the top horizontal line in panel (e). As the coupling strength g increases, the map with 3 stable FPs: blue PM at (0.5, 0.5), green PM at (0.5, 0), and red PM at (0, 0.5) in (A1) and gains two more stable FPs: black TW at (0.33, 0.67) and purple TW (0.67, 0.33) in (A). The growing TW and shrinking PM attraction basins are separated by the separatrices of the saddles (gray ) in (A3). Each stable PM becomes a saddle through a pitch-fork bifurcation at larger g-values, after merging with two nearby saddles in (A4). Panels (a1)–(d1) snapshot several maps as the network transitions from the escape to the release mechanism as Iapp is decreased [along the left vertical dotted lines in (e)]: from 3 stable PMs in (A1), to two stable TWs in (B1), to mixed dynamics with both PM/TWs in (C1), and back to 3 stable PMs in (D1). Parameters for (A1)–(A4): Iapp=0.5886, g=(0.0015,0.006,0.019,0.0225); parameters for (B1)–(D1): g=0.0015, Iapp=(0.493,0.419,0.393). The return maps sampled along the horizontal dashed lines A–D and A-D at the bottom of panel (E) are presented in Figs. 7 and 8 for the networks obeying the release mechanism.

Close modal
FIG. 5.

Key asymmetric configurations analyzed in this study: (a) Mono-biased motif where only the synapse g31 is manipulated, (b) Double-biased motif, in which the reciprocal connections between cells 1 and 3 (g31 and g13) are manipulated equally, (c) Driver-biased motif, in which both the outgoing connections from cell 3 (g31 and g32) are affected equally, and (d) Clockwise-biased motif, in which all the clockwise connections (g12, g23, and g31) are affected equally.

FIG. 5.

Key asymmetric configurations analyzed in this study: (a) Mono-biased motif where only the synapse g31 is manipulated, (b) Double-biased motif, in which the reciprocal connections between cells 1 and 3 (g31 and g13) are manipulated equally, (c) Driver-biased motif, in which both the outgoing connections from cell 3 (g31 and g32) are affected equally, and (d) Clockwise-biased motif, in which all the clockwise connections (g12, g23, and g31) are affected equally.

Close modal

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˙=0, relative to that of the slow h-nullcline given by h˙=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.

FIG. 6.

Phase portraits of the gFN-neuron given by Eq. (1), demonstrating the release (a) and escape (b) mechanisms of rhythmogenesis in the 3-cell circuit. The clusters of the colored dots representing the phases of the coupled cells on the periodic orbit, at the lower-left and upper-right knees near the tangencies of the fast cubic and slow sigmoidal nullclines V˙=0 and h˙=0, respectively, are indicative of the stagnation due to the saddle-node bottleneck effect underlying the release and escape mechanisms. The solid (unperturbed) and dashed (perturbed) purple curves show the shifts in the V-nullclines for the pre- and post-synaptic cells, respectively. The intersection of the fast dashed V-nullcline in (a) with the slow h-nullcline near the left knee corresponds to a newly formed stable fixed point that makes the post-synaptic cell “hard-locked” by synaptic inhibition. In case (b), the post-synaptic cell escapes from the trap (hard or soft) of depolarization near the upper knee, after an inhibitory perturbation shifts V˙=0 leftward to create a gap between the nullclines, allowing the cell(s) to fall down onto the hyperpolarized branch, to start a new cycle of revolution.

FIG. 6.

Phase portraits of the gFN-neuron given by Eq. (1), demonstrating the release (a) and escape (b) mechanisms of rhythmogenesis in the 3-cell circuit. The clusters of the colored dots representing the phases of the coupled cells on the periodic orbit, at the lower-left and upper-right knees near the tangencies of the fast cubic and slow sigmoidal nullclines V˙=0 and h˙=0, respectively, are indicative of the stagnation due to the saddle-node bottleneck effect underlying the release and escape mechanisms. The solid (unperturbed) and dashed (perturbed) purple curves show the shifts in the V-nullclines for the pre- and post-synaptic cells, respectively. The intersection of the fast dashed V-nullcline in (a) with the slow h-nullcline near the left knee corresponds to a newly formed stable fixed point that makes the post-synaptic cell “hard-locked” by synaptic inhibition. In case (b), the post-synaptic cell escapes from the trap (hard or soft) of depolarization near the upper knee, after an inhibitory perturbation shifts V˙=0 leftward to create a gap between the nullclines, allowing the cell(s) to fall down onto the hyperpolarized branch, to start a new cycle of revolution.

Close modal

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 AIII 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.

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 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)–7(c). TWs finally become unstable through a secondary Andronov–Hopf/torus bifurcation in Fig. 7(d) and hence not observable in the network anymore.

FIG. 7.

Transformations in the 2D return maps corresponding to the fully symmetric networks with the release mechanism as the inhibitory coupling g is increased along the dashed white line (a)–(d) in the bifurcation diagram in Fig. 4(e). The TW attraction basins decrease in size, while the PM basins increase from (a) to (c). The TWs lose stability in (d) through a secondary Andronov–Hopf/torus bifurcation. Less smooth trajectories are indicative of fast convergence to the attractors in the network, at greater coupling strengths. Parameters: Iapp=0.4155, g=(0.0005,0.006,0.015,0.018).

FIG. 7.

Transformations in the 2D return maps corresponding to the fully symmetric networks with the release mechanism as the inhibitory coupling g is increased along the dashed white line (a)–(d) in the bifurcation diagram in Fig. 4(e). The TW attraction basins decrease in size, while the PM basins increase from (a) to (c). The TWs lose stability in (d) through a secondary Andronov–Hopf/torus bifurcation. Less smooth trajectories are indicative of fast convergence to the attractors in the network, at greater coupling strengths. Parameters: Iapp=0.4155, g=(0.0005,0.006,0.015,0.018).

Close modal
FIG. 8.

Poincaré return maps corresponding to the fully symmetric 3-cell motif with another example of the release mechanism for the parameter values sampled along the line (a)–(d) in the bifurcation diagram [Fig. 4(e)]. As g is increased, the unstable TWs in (a)–(b) become stable in panel c through a torus bifurcation. The TWs again lose stability with a further increase in g, leading to hard-locking in the system that results in the quick and jagged convergence to the three PMs in (d). Parameters: Iapp=0.3956, g=(0.0005,0.005,0.007,0.015).

FIG. 8.

Poincaré return maps corresponding to the fully symmetric 3-cell motif with another example of the release mechanism for the parameter values sampled along the line (a)–(d) in the bifurcation diagram [Fig. 4(e)]. As g is increased, the unstable TWs in (a)–(b) become stable in panel c through a torus bifurcation. The TWs again lose stability with a further increase in g, leading to hard-locking in the system that results in the quick and jagged convergence to the three PMs in (d). Parameters: Iapp=0.3956, g=(0.0005,0.005,0.007,0.015).

Close modal

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) and 8(b) to a configuration that supports both PMs and TWs in Fig. 8(c), following torus bifurcations. A further increase in the value of g leads to hard-locking in Fig. 8(d), 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.

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–F 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.

FIG. 9.

(g31,Iapp)-bifurcation diagram of the mono-biased motif presented in Fig. 5(a) shows five distinct regions corresponding to PM, TW, and PS rhythms along with the combinations PM/PS and PM/TW. Transitions between these regions are governed by saddle-node (SN) bifurcations that eliminate or restore FPs in the return map. The points A–F and A–F highlighted near the top dotted and bottom dashed lines, respectively, indicate the parameter values used for the return maps with the escape and release mechanisms elaborated in Figs. 10 and 11, respectively. The vertical line given by g31=gall=0.001 corresponds to the fully symmetric network where all g-values are identical.

FIG. 9.

(g31,Iapp)-bifurcation diagram of the mono-biased motif presented in Fig. 5(a) shows five distinct regions corresponding to PM, TW, and PS rhythms along with the combinations PM/PS and PM/TW. Transitions between these regions are governed by saddle-node (SN) bifurcations that eliminate or restore FPs in the return map. The points A–F and A–F highlighted near the top dotted and bottom dashed lines, respectively, indicate the parameter values used for the return maps with the escape and release mechanisms elaborated in Figs. 10 and 11, respectively. The vertical line given by g31=gall=0.001 corresponds to the fully symmetric network where all g-values are identical.

Close modal
FIG. 10.

Poincaré return maps corresponding to the mono-biased motif with the release mechanism and their evolution as g31 gradually increases, while the remaining synaptic strengths (gall) are held constant. When g31=0 in Panel (a), the green PM [stable FP at (0.5,0)] dominates the dynamics, coexisting with the blue PM [FP at (0.5,0.5)] having a smaller attraction basin. Increasing g31 leads to the formation of the purple TW pattern [stable FP at (0.6,0.3)] through a saddle-node bifurcation, while the blue basin increases in Panel (b), followed by the appearance of the red PM in (c) due to a saddle-node bifurcation, after the motif partially restores the anti-clockwise symmetry. Increasing g31 further leads to a single dominant anti-clockwise purple TW rhythm in the network, after all other attracting FPs vanish via a series of saddle-node bifurcations in Panels (d)–(f). Parameters: Iapp=0.412, gall=0.001, g31=(0, 0.000 81, 0.001 08, 0.001 35, 0.004, 0.008).

FIG. 10.

Poincaré return maps corresponding to the mono-biased motif with the release mechanism and their evolution as g31 gradually increases, while the remaining synaptic strengths (gall) are held constant. When g31=0 in Panel (a), the green PM [stable FP at (0.5,0)] dominates the dynamics, coexisting with the blue PM [FP at (0.5,0.5)] having a smaller attraction basin. Increasing g31 leads to the formation of the purple TW pattern [stable FP at (0.6,0.3)] through a saddle-node bifurcation, while the blue basin increases in Panel (b), followed by the appearance of the red PM in (c) due to a saddle-node bifurcation, after the motif partially restores the anti-clockwise symmetry. Increasing g31 further leads to a single dominant anti-clockwise purple TW rhythm in the network, after all other attracting FPs vanish via a series of saddle-node bifurcations in Panels (d)–(f). Parameters: Iapp=0.412, gall=0.001, g31=(0, 0.000 81, 0.001 08, 0.001 35, 0.004, 0.008).

Close modal
FIG. 11.

Poincaré return maps corresponding to the mono-biased motif with the escape mechanism as g31 gradually increases. At g31=0, there is a single dominant clockwise (black) TW FP in panel (a); this is elaborated further in Figs. 12 and 13. (b) As g31 increases, the unstable IC (the blue line in Fig. 12 ) first undergoes a reverse homoclinic saddle-node bifurcation, giving rise to a repelling fixed point (white dot) and a saddle, which then undergoes a pitch-fork bifurcation that makes it stable – the blue PM, with two additional saddles. The red and green PMs then emerge following additional saddle-node bifurcations, see panels (c) and (d). The green PM disappears through a saddle-node bifurcation, while the black TW at (0.3, 0.6) becomes repelling via a torus bifurcation, giving rise to a stable invariant circle (IC) in panel (e), and finally gets annihilated after merging with a nearby saddle. Parameters: Iapp=0.5825, gall=0.001, g31=(0, 0.000 27, 0.000 676, 0.000 81, 0.001 49, 0.004 05).

FIG. 11.

Poincaré return maps corresponding to the mono-biased motif with the escape mechanism as g31 gradually increases. At g31=0, there is a single dominant clockwise (black) TW FP in panel (a); this is elaborated further in Figs. 12 and 13. (b) As g31 increases, the unstable IC (the blue line in Fig. 12 ) first undergoes a reverse homoclinic saddle-node bifurcation, giving rise to a repelling fixed point (white dot) and a saddle, which then undergoes a pitch-fork bifurcation that makes it stable – the blue PM, with two additional saddles. The red and green PMs then emerge following additional saddle-node bifurcations, see panels (c) and (d). The green PM disappears through a saddle-node bifurcation, while the black TW at (0.3, 0.6) becomes repelling via a torus bifurcation, giving rise to a stable invariant circle (IC) in panel (e), and finally gets annihilated after merging with a nearby saddle. Parameters: Iapp=0.5825, gall=0.001, g31=(0, 0.000 27, 0.000 676, 0.000 81, 0.001 49, 0.004 05).

Close modal

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) 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)]. With further increases in the strength of g31 through 0.00135,0.004,and0.008 in Figs. 10(d)–10(f), 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 CIII H.

FIG. 12.

Four identical panels from Fig. 11(a) stitched together for a better understanding of the trajectories wrapping around the phase torus. The network is dominated by a single black TW around (0.3, 0.6). The saddle, located in a close proximity, causes two trajectories from close initial conditions to traverse different paths (red lines) to converge to the same FP. Shown in blue is the repelling invariant circle (IC), see also Fig. 13. Parameters: Iapp=0.5825, gall=0.001 except g31=0.

FIG. 12.

Four identical panels from Fig. 11(a) stitched together for a better understanding of the trajectories wrapping around the phase torus. The network is dominated by a single black TW around (0.3, 0.6). The saddle, located in a close proximity, causes two trajectories from close initial conditions to traverse different paths (red lines) to converge to the same FP. Shown in blue is the repelling invariant circle (IC), see also Fig. 13. Parameters: Iapp=0.5825, gall=0.001 except g31=0.

Close modal
FIG. 13.

Four identical panels stitched together to better visualize the trajectories wrapping around the phase torus, as a stable PS pattern emerges following the disappearance of the black TW and the nearby saddle of Fig. 12, through a homoclinic saddle-node bifurcation. The network has no phase-locked rhythms and all trajectories converge on to the stable invariant curve (red). The blue line marks an unstable invariant curve. Parameters: Iapp=0.5875, gall=0.001 except g31=0.

FIG. 13.

Four identical panels stitched together to better visualize the trajectories wrapping around the phase torus, as a stable PS pattern emerges following the disappearance of the black TW and the nearby saddle of Fig. 12, through a homoclinic saddle-node bifurcation. The network has no phase-locked rhythms and all trajectories converge on to the stable invariant curve (red). The blue line marks an unstable invariant curve. Parameters: Iapp=0.5875, gall=0.001 except g31=0.

Close modal

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).

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.

FIG. 14.

The (g31/13,Iapp)-bifurcation diagram of the double-biased motif. 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, and TW/PS in the color-mapped regions. The PM and PS behaviors dominate at weak and strong coupling, respectively. Points A–D and A–F indicate the sampled parameter values for the release (Fig. 15) and the escape (Fig. 16) mechanisms, respectively. The vertical line is where the network is symmetric at g31=gall=0.001.

FIG. 14.

The (g31/13,Iapp)-bifurcation diagram of the double-biased motif. 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, and TW/PS in the color-mapped regions. The PM and PS behaviors dominate at weak and strong coupling, respectively. Points A–D and A–F indicate the sampled parameter values for the release (Fig. 15) and the escape (Fig. 16) mechanisms, respectively. The vertical line is where the network is symmetric at g31=gall=0.001.

Close modal

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), 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), 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) 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) 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.

FIG. 15.

Poincaré return maps corresponding to the double-biased motif with the release mechanism. When the synapses g31 and g13 are weak, there is only the green PM in the map shown in panel (a). As their strength is increased, the blue and red PMs emerge through a series of saddle-node and pitch-fork bifurcations, along with multiple saddles and two repelling FPs corresponding to the unstable TWs in panel (b). Further increase in the synaptic strengths causes the unstable TW FPs to disappear through a heteroclinic saddle-node bifurcation in panel (c). The blue and red PMs also disappear through a heteroclinic saddle-node bifurcation, thus giving rise to a stable “phase-slipping” invariant circle (shown gray), coexisting with the green PM in panel (d). Their basins are partitioned by the incoming separatrices (black curves) of the saddles. Parameters: Iapp=0.399, gall=0.001 except for g31=g13= (0.0005, 0.001, 0.0012, 0.0015).

FIG. 15.

Poincaré return maps corresponding to the double-biased motif with the release mechanism. When the synapses g31 and g13 are weak, there is only the green PM in the map shown in panel (a). As their strength is increased, the blue and red PMs emerge through a series of saddle-node and pitch-fork bifurcations, along with multiple saddles and two repelling FPs corresponding to the unstable TWs in panel (b). Further increase in the synaptic strengths causes the unstable TW FPs to disappear through a heteroclinic saddle-node bifurcation in panel (c). The blue and red PMs also disappear through a heteroclinic saddle-node bifurcation, thus giving rise to a stable “phase-slipping” invariant circle (shown gray), coexisting with the green PM in panel (d). Their basins are partitioned by the incoming separatrices (black curves) of the saddles. Parameters: Iapp=0.399, gall=0.001 except for g31=g13= (0.0005, 0.001, 0.0012, 0.0015).

Close modal

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.

FIG. 16.

Poincaré return maps corresponding to the double-biased motif with the escape mechanism. When both synapses g31=g13=0, the network produces only the green PM rhythm in panel (a) (elaborated further in Fig. 17). The blue and red PMs emerge following pitch-fork bifurcations of the saddles in panel (b). The two TWs emerging in panel (c), disappear through saddle-node bifurcations in panel (d). Next, the green PM disappears via a saddle-node bifurcation, after merging with a nearby saddle in panel (e). The red and blue PMs finally disappear through a heteroclinic saddle-node bifurcation, which gives rise to the only stable invariant circle (gray) in panel (f), moving in the opposite direction compared with Fig. 15(d). Parameters: Iapp=0.5716, gall=0.001, except for g31=g13= (0, 0.000 676, 0.001, 0.001 28, 0.001 55, 0.003 31).

FIG. 16.

Poincaré return maps corresponding to the double-biased motif with the escape mechanism. When both synapses g31=g13=0, the network produces only the green PM rhythm in panel (a) (elaborated further in Fig. 17). The blue and red PMs emerge following pitch-fork bifurcations of the saddles in panel (b). The two TWs emerging in panel (c), disappear through saddle-node bifurcations in panel (d). Next, the green PM disappears via a saddle-node bifurcation, after merging with a nearby saddle in panel (e). The red and blue PMs finally disappear through a heteroclinic saddle-node bifurcation, which gives rise to the only stable invariant circle (gray) in panel (f), moving in the opposite direction compared with Fig. 15(d). Parameters: Iapp=0.5716, gall=0.001, except for g31=g13= (0, 0.000 676, 0.001, 0.001 28, 0.001 55, 0.003 31).

Close modal
FIG. 17.

Four identicals panels from Fig. 16(a) stitched together to continuously visualize trajectories wrapping around the 2D torus. This network is monostable with a single PM rhythm—the green FP at (0.5,0), to which some trajectories converge along quite a long path, as its attraction basin is shaped by complex interactions of the separatrices (black lines) of the two saddles (black ).

FIG. 17.

Four identicals panels from Fig. 16(a) stitched together to continuously visualize trajectories wrapping around the 2D torus. This network is monostable with a single PM rhythm—the green FP at (0.5,0), to which some trajectories converge along quite a long path, as its attraction basin is shaped by complex interactions of the separatrices (black lines) of the two saddles (black ).

Close modal

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.

FIG. 18.

(g31/32,Iapp)-bifurcation diagram for the driver-biased motif. The network is dominated by the red PM for strong synaptic coupling, as expected from the asymmetry. The network also produces TWs, PS, and a combination of PMs with TWs or PS at weaker coupling, along with the blue and green PMs. The vertical line represents the network symmetry where g31=g32=gall=0.001. Sampled parameter values for the release and escape mechanism are as described previously and elaborated in Figs. 19 and 20, respectively.

FIG. 18.

(g31/32,Iapp)-bifurcation diagram for the driver-biased motif. The network is dominated by the red PM for strong synaptic coupling, as expected from the asymmetry. The network also produces TWs, PS, and a combination of PMs with TWs or PS at weaker coupling, along with the blue and green PMs. The vertical line represents the network symmetry where g31=g32=gall=0.001. Sampled parameter values for the release and escape mechanism are as described previously and elaborated in Figs. 19 and 20, respectively.

Close modal

Figure 19 shows the return maps for the cells obeying the release mechanism. For small values of g31=g32 in Fig. 19(a), the blue and green PMs coexist. As the synaptic strengths are increased in Fig. 19(b), 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), as well as the disappearance of the two TWs in Fig. 19(d) through saddle-node bifurcations, resulting in a single red PM rhythm of the circuit, under the control of the biased-driver cell.

FIG. 19.

Poincaré return maps corresponding to the driver-biased motif with the release mechanism. For small g31=g32-values, the network produces both the blue and green PMs in panel (a). As their value is increased, the red PM and the two TWs then emerge via a series of saddle-node bifurcations in panel (b). Next, the green and the blue PMs disappear in (c), as well as the black and purple traveling waves in (d), via saddle-node bifurcations, resulting in the single dominant red PM. Parameters: Iapp=0.426, gall=0.001 except for g31=g32=(0.0001, 0.001, 0.001 15, 0.0015).

FIG. 19.

Poincaré return maps corresponding to the driver-biased motif with the release mechanism. For small g31=g32-values, the network produces both the blue and green PMs in panel (a). As their value is increased, the red PM and the two TWs then emerge via a series of saddle-node bifurcations in panel (b). Next, the green and the blue PMs disappear in (c), as well as the black and purple traveling waves in (d), via saddle-node bifurcations, resulting in the single dominant red PM. Parameters: Iapp=0.426, gall=0.001 except for g31=g32=(0.0001, 0.001, 0.001 15, 0.0015).

Close modal
FIG. 20.

Poincaré return maps corresponding to the driver-biased motif with the escape mechanism: with weak synapses g31=g32, the network is dominated by a single phase-slipping rhythm in panel (a). As they become stronger, through a series of saddle-node bifurcations, the red PM first emerges in (b), followed by the blue and green PMs and the two TWs in (c). Furthermore, the blue and green PMs then disappear in (d), followed by the loss of the two TWs in (e), via a heteroclinic saddle-node bifurcation, giving rise to a stable PS pattern—the stable IC (gray basin) that coexists with the red PM. The IC then disappears and the network becomes mono-stable with the single dominant PM rhythm corresponding to the red FP in (f). Parameters: Iapp=0.57, gall=0.001 except for g31=g32=(0.000 01, 0.000 65, 0.001, 0.0011, 0.001 36, 0.0025).

FIG. 20.

Poincaré return maps corresponding to the driver-biased motif with the escape mechanism: with weak synapses g31=g32, the network is dominated by a single phase-slipping rhythm in panel (a). As they become stronger, through a series of saddle-node bifurcations, the red PM first emerges in (b), followed by the blue and green PMs and the two TWs in (c). Furthermore, the blue and green PMs then disappear in (d), followed by the loss of the two TWs in (e), via a heteroclinic saddle-node bifurcation, giving rise to a stable PS pattern—the stable IC (gray basin) that coexists with the red PM. The IC then disappears and the network becomes mono-stable with the single dominant PM rhythm corresponding to the red FP in (f). Parameters: Iapp=0.57, gall=0.001 except for g31=g32=(0.000 01, 0.000 65, 0.001, 0.0011, 0.001 36, 0.0025).

Close modal

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.

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.

FIG. 21.

g,Iapp-bifurcation diagram of the clockwise-biased motif. The network is dominated primarily by traveling wave rhythms as expected from this asymmetry. A single TW is seen for both small and large g-values, while both the TWs are seen for moderate values. PMs and PM/TWs are also seen close to symmetry in the network, indicated by the vertical dotted line where g=gall=0.001. Parameter values sampled for the release and the escape mechanisms in Figs. 22 and 23 are also shown.

FIG. 21.

g,Iapp-bifurcation diagram of the clockwise-biased motif. The network is dominated primarily by traveling wave rhythms as expected from this asymmetry. A single TW is seen for both small and large g-values, while both the TWs are seen for moderate values. PMs and PM/TWs are also seen close to symmetry in the network, indicated by the vertical dotted line where g=gall=0.001. Parameter values sampled for the release and the escape mechanisms in Figs. 22 and 23 are also shown.

Close modal

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). Following a series of saddle-node bifurcations, the three PMs emerge in Fig. 22(b). This is followed by the emergence of the anti-clockwise TW (purple) in Fig. 22(c) and the disappearance of the clockwise TW (black) in Fig. 22(d), through their respective torus bifurcations. The three PMs then disappear via saddle-node bifurcations in Fig. 22(e), leading to a single anti-clockwise TW rhythm dominating the network. Further strengthening of the synapses in Fig. 22(f) leads to hard locking, with a single TW rhythm.

FIG. 22.

Poincaré return maps corresponding to the clockwise-biased motif with the release mechanism. The network is initially dominated by a single stable clockwise TW in panel (a), while the purple TW FP remains unstable. As the clockwise synapses are strengthened, the three PMs emerge in panel (b) via saddle-node bifurcations around the black TW FP. The purple TW then becomes stable in (c), while the black TW in (d) loses stability, via torus bifurcations. The three PMs disappear via saddle-node bifurcations in panel (e), resulting in the single dominant purple TW rhythm that finally becomes hard-locked in (f) as indicated by the rapid, non-smooth convergence of trajectories to it. Parameters: Iapp=0.4, gall=0.001 except g12=g23=g31= (0.0005, 0.000 65, 0.001, 0.001 35, 0.0016, 0.002).

FIG. 22.

Poincaré return maps corresponding to the clockwise-biased motif with the release mechanism. The network is initially dominated by a single stable clockwise TW in panel (a), while the purple TW FP remains unstable. As the clockwise synapses are strengthened, the three PMs emerge in panel (b) via saddle-node bifurcations around the black TW FP. The purple TW then becomes stable in (c), while the black TW in (d) loses stability, via torus bifurcations. The three PMs disappear via saddle-node bifurcations in panel (e), resulting in the single dominant purple TW rhythm that finally becomes hard-locked in (f) as indicated by the rapid, non-smooth convergence of trajectories to it. Parameters: Iapp=0.4, gall=0.001 except g12=g23=g31= (0.0005, 0.000 65, 0.001, 0.001 35, 0.0016, 0.002).

Close modal

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).

FIG. 23.

Poincaré return maps corresponding to the clockwise-biased motif with the escape mechanism look almost identical to those with the release mechanism shown in Fig. 22. The network transitions from a single dominant black TW in panel (a) to purple TW in panel (f). Along the way, the three PMs emerge via saddle-node bifurcations in (b), the black TW loses stability in (c), while the purple TW becomes stable in panels (d) and (e) via torus bifurcations (resulting in the transitory invariant circle shown), and the three PMs are lost via saddle-node bifurcations (f). Parameters: Iapp=0.5886, gall=0.001 except g12=g23=g31= (0.0005, 0.000 743, 0.001, 0.001 16, 0.001 32, 0.002).

FIG. 23.

Poincaré return maps corresponding to the clockwise-biased motif with the escape mechanism look almost identical to those with the release mechanism shown in Fig. 22. The network transitions from a single dominant black TW in panel (a) to purple TW in panel (f). Along the way, the three PMs emerge via saddle-node bifurcations in (b), the black TW loses stability in (c), while the purple TW becomes stable in panels (d) and (e) via torus bifurcations (resulting in the transitory invariant circle shown), and the three PMs are lost via saddle-node bifurcations (f). Parameters: Iapp=0.5886, gall=0.001 except g12=g23=g31= (0.0005, 0.000 743, 0.001, 0.001 16, 0.001 32, 0.002).

Close modal

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 breakdown49 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.

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 (0.5,0) and (0.5,0.5), respectively. In addition to the persistent repeller at the origin, there are three more saddle FPs (labelled by ’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 ω-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 Δ12=0 and comes back into the map from the right wall given Δ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 ω-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 ω-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.

FIG. 24.

Six maps demonstrating the stages in the emergence of a transitive torus in the given asymmetric motif (top). By increasing Iapp, the network migrates from the release to the escape mechanisms. Both forward (gray) and backward (red) trajectories elucidate the unstable FPs and ICs. Step 1: The network initially has a single stable phase-slipping pattern (vertical gray IC with the decreasing Δ13-direction) in panel (a), along with a repelling FP near the origin and a saddle. Step 2: both merge and vanish through a homoclinic saddle-node bifurcation to give rise to an additional repelling PS pattern (red IC with the increasing Δ13-direction) in panel (b). Step 3: the stable IC breaks down via a homoclinic saddle-node bifurcation as two FPs, a stable one and a saddle, emerge on it so that the outgoing separatrices of the saddle end up on the stable FP around (0.3, 0.6). The gap between both FPs increases, which makes them closer to each other after crossing the boundary Δ13=0=1. Step 4: The FPs disappear via a reverse homoclinic saddle-node bifurcation, and the stable PS pattern re-emerges with the opposite direction, matching the orientation of the unstable IC in panel (d). Step 5: the distance between the stable and the unstable ICs decreases, as they start to merge and vanish in panel (e). This completes the bifurcation sequence giving rise to a transitive torus in (f), without any FPs or ICs such that a single trajectory densely fills its surface. Parameters: g13=g31=0.0003, g23=g32=0.0005, g12=001, g21=0.0051, Iapp= (0.4, 0.46, 0.5, 0.572 982, 0.594, 0.61).

FIG. 24.

Six maps demonstrating the stages in the emergence of a transitive torus in the given asymmetric motif (top). By increasing Iapp, the network migrates from the release to the escape mechanisms. Both forward (gray) and backward (red) trajectories elucidate the unstable FPs and ICs. Step 1: The network initially has a single stable phase-slipping pattern (vertical gray IC with the decreasing Δ13-direction) in panel (a), along with a repelling FP near the origin and a saddle. Step 2: both merge and vanish through a homoclinic saddle-node bifurcation to give rise to an additional repelling PS pattern (red IC with the increasing Δ13-direction) in panel (b). Step 3: the stable IC breaks down via a homoclinic saddle-node bifurcation as two FPs, a stable one and a saddle, emerge on it so that the outgoing separatrices of the saddle end up on the stable FP around (0.3, 0.6). The gap between both FPs increases, which makes them closer to each other after crossing the boundary Δ13=0=1. Step 4: The FPs disappear via a reverse homoclinic saddle-node bifurcation, and the stable PS pattern re-emerges with the opposite direction, matching the orientation of the unstable IC in panel (d). Step 5: the distance between the stable and the unstable ICs decreases, as they start to merge and vanish in panel (e). This completes the bifurcation sequence giving rise to a transitive torus in (f), without any FPs or ICs such that a single trajectory densely fills its surface. Parameters: g13=g31=0.0003, g23=g32=0.0005, g12=001, g21=0.0051, Iapp= (0.4, 0.46, 0.5, 0.572 982, 0.594, 0.61).

Close modal
FIG. 25.

2D return maps showing how the outgoing and the incoming separatrices (gray lines wrapping around the torus) of the two saddles (labelled with gray ) shape the boundaries and the sizes of the attraction basins of the stable FPs, green and blue. The direction in which the heteroclinic connection between the two saddle FPs splits determines which stable FP has the largest attraction basin: blue in panel (a) at g31=0.0071 or green in panel B at g31=0.0072; other parameters: g12=g32=0.0038 and g21=g23=g13=0.0041 and ε=0.3.

FIG. 25.

2D return maps showing how the outgoing and the incoming separatrices (gray lines wrapping around the torus) of the two saddles (labelled with gray ) shape the boundaries and the sizes of the attraction basins of the stable FPs, green and blue. The direction in which the heteroclinic connection between the two saddle FPs splits determines which stable FP has the largest attraction basin: blue in panel (a) at g31=0.0071 or green in panel B at g31=0.0072; other parameters: g12=g32=0.0038 and g21=g23=g13=0.0041 and ε=0.3.

Close modal

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 previously50 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 Δ12=Δ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=0 and the slow h=0, this newly formed attractor co-exists with the stable PMs (green, red, and blue) in Fig. 26(a1) or with stable PMs and TWs (black and purple spirals) in Fig. 26(a2).

FIG. 26.

(a)–(c) Snapshot of MotifToolbox43 panels showing the Poincaré return maps with an additional stable “synchronous” FP at the origin (0,0) (yellow), along with three stable PMs in panel (a1) (parameters: gij=0.0035, I=0.409, V0=0.113, ε=0.3) or three stable PMs and two stable TWs in panel (a2) (parameters: gij=0.0021, I=0.406, V0=0.113, ε=0.3). Panels (b1) and (b2) reveal the attraction basins of the coexisting FPs in greater detail and allow us to determine the locations of the repelling FPs, at the junctions of three distinct color-coded regions. Traces in panel (c1) demonstrate an asymptotic convergence to the synchronous rhythm corresponding to the FP at the origin. (d) Magnification of a neighborhood of the stable synchronous state (yellow) in the return map presented in (a1), by stitching together four identical panels, to disclose its structure with six repellers (white dots) and six saddles (gray diamonds) surrounding the stable FP (yellow) at the origin (0,0).

FIG. 26.

(a)–(c) Snapshot of MotifToolbox43 panels showing the Poincaré return maps with an additional stable “synchronous” FP at the origin (0,0) (yellow), along with three stable PMs in panel (a1) (parameters: gij=0.0035, I=0.409, V0=0.113, ε=0.3) or three stable PMs and two stable TWs in panel (a2) (parameters: gij=0.0021, I=0.406, V0=0.113, ε=0.3). Panels (b1) and (b2) reveal the attraction basins of the coexisting FPs in greater detail and allow us to determine the locations of the repelling FPs, at the junctions of three distinct color-coded regions. Traces in panel (c1) demonstrate an asymptotic convergence to the synchronous rhythm corresponding to the FP at the origin. (d) Magnification of a neighborhood of the stable synchronous state (yellow) in the return map presented in (a1), by stitching together four identical panels, to disclose its structure with six repellers (white dots) and six saddles (gray diamonds) surrounding the stable FP (yellow) at the origin (0,0).

Close modal

This figure also introduces another concept incorporated into MotifToolbox43 to reveal the attraction basins of the co-existing stable FPs (labelled by color-matching dots) in greater details, as presented in Fig. 26(b1) and 26(b2). Using these diagrams, we can identify the locations of the repellers in the (Δ12,Δ13)-plane at the junctions of three distinctively colored regions. For example, in Figs. 26(a1) and 26(b1), 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(a2) and 26(b2) 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 Δ12=0, Δ13=0, and Δ12=Δ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.

The 2θ-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” θ-neuron54 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 θ-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=0 and h=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

FIG. 27.

Unit circle represents the phase space of a θ-neuron with a single “quiescent” saddle-node phantom at 0 (near the quadratic tangency at the bottom) (a) and a 2θ-neuron with two single saddle-node phantoms at 0 (quiescent phase) and π (tonic-spiking phase) in (c), where the counter-clockwise rotations slow down as illustrated by the spheres representing the phase points on S1. Varying the gaps or the distances from/to the saddle-nodes let one control the inter-spike/burst interval (b) or the duty cycle (d) of bursting traces.

FIG. 27.

Unit circle represents the phase space of a θ-neuron with a single “quiescent” saddle-node phantom at 0 (near the quadratic tangency at the bottom) (a) and a 2θ-neuron with two single saddle-node phantoms at 0 (quiescent phase) and π (tonic-spiking phase) in (c), where the counter-clockwise rotations slow down as illustrated by the spheres representing the phase points on S1. Varying the gaps or the distances from/to the saddle-nodes let one control the inter-spike/burst interval (b) or the duty cycle (d) of bursting traces.

Close modal

A key feature of the 2θ-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θ-neuron: “on” at π 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 θ-model phenomenologically depicts spiking cells, while the 2θ-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θ-bursters preserve all the key features seen in a motif composed of the three gFN-neurons as well.

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

{θ1=ωcos2θ1+αcosθ1(β211+ekcosθ2+β311+ekcosθ3)[121+eksinθ1],θ2=ωcos2θ2+αcosθ2(β121+ekcosθ1+β321+ekcosθ3)[121+eksinθ2],θ3=ωcos2θ3+αcosθ3(β131+ekcosθ1+β231+ekcosθ2)[121+eksinθ3],mod 1.
(4)

One can observe that the phase dynamics of an individual 2θ-neuron are governed by the terms ωcos(2θ). As long as the frequency 0<ω1, there are two stable and two unstable equilibria: at the bottom around θ0 and at the top near θπ; they are associated with the hyperpolarized and the depolarized quiescent states of the neuron. When ω>1, the 2θ-neuron becomes oscillatory through two simultaneous saddle-node bifurcations (SNIC) on a unit circle S1, where θ is defined on modulo one. Moreover, whenever ω=1+ε, where 0<ε1, this new “burster” possesses two slow phases: the active “on” state near θ=π 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 ω, the active and inactive phases should be defined by π/2<θ3π/2 and 3π/2<θπ/2, respectively. The latter phase is below the synaptic threshold, which is set by θth=π/2 so that cos(θth)=0, thus equally dividing the unit circle. The duty cycle of the 2θ-neuron is controlled by the term αcosθ, provided that it remains oscillatory as long as ω|α|>1. Note that when α=0, the duty cycle is 50% and the oscillations are even. The active or inactive phases can be extended or shortened, respectively, with α<0, making the duty cycle greater or vice versa—the duty cycle of individual neurons can be decreased with α>0.

The 2θ-neurons are coupled in the 3-cell network using the fast inhibitory FTM42 synapses. The “sigmoidal” term 11+ekcosθi, ranging between 1 and 0, rapidly triggers (here, k=10) an influx of inhibition flowing from the ith pre-synaptic neuron into the jth post-synaptic neuron, as soon as the former enters the active on-phase above the synaptic threshold cos(θth)=0, i.e., π/2<θi<3π/2. The strength of the inhibitory coupling is determined by the maximal conductance βij that slows down the rate of increase of θj in the jth 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 θ=0 and θ=π. This is achieved by multiplying all the coupling terms of each ODE by [121+eksinθ]. When 0<θ<π, the inhibition is a negative input to slow the transition into bursting. When π<θi<2π, 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 [111+eksinθ], 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θ-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 θ=π, above the synaptic threshold, inhibits and pushes the other two closer to each other, near the bottom quiescent state at θ=0, by accelerating the red neuron on the downstroke and by slowing down the blue neuron on the upstroke.

FIG. 28.

(a) Phase lags between the three 2θ-neurons are recorded when the phase/voltage reaches the synaptic threshold θsyn=0 from below. (b) Phase progressions of the coupled 2θ-neurons and their color-coded phase points on a unit circle S1.

FIG. 28.

(a) Phase lags between the three 2θ-neurons are recorded when the phase/voltage reaches the synaptic threshold θsyn=0 from below. (b) Phase progressions of the coupled 2θ-neurons and their color-coded phase points on a unit circle S1.

Close modal

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θ-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θ-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.

FIG. 29.

(a) Time progression of multiple initial conditions exponentially converging to several phase-locked states [four of them shown in (c1) and (c2)], generated by a symmetric 3-cell motif [see Fig. 2(a)]. These phase-locked states correspond to the co-existing stable FPs of the 2D Poincaré map depicted in panel (d), wrapping around the 2D torus in panel (b).

FIG. 29.

(a) Time progression of multiple initial conditions exponentially converging to several phase-locked states [four of them shown in (c1) and (c2)], generated by a symmetric 3-cell motif [see Fig. 2(a)]. These phase-locked states correspond to the co-existing stable FPs of the 2D Poincaré map depicted in panel (d), wrapping around the 2D torus in panel (b).

Close modal

One can see that the 2D return maps for the gFN-neurons and 2θ-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θ-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θ-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 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θ-bursters, as the phases on S1 will just reverse the direction and spin clockwise on the unit circle.

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θ-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 tongues59,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.

FIG. 30.

Diagram outlining several nested synchronization zones in the parameter space of a 3-cell neural network. Pairs of stable and unstable FPs in the 2D return map are eliminated sequentially (see Figs. 10, 11, 15, and 16) when the boundaries of these zones, corresponding to saddle-node (SN) bifurcations, are crossed outwardly as the coupling Δgij between neurons i and j is increased.

FIG. 30.

Diagram outlining several nested synchronization zones in the parameter space of a 3-cell neural network. Pairs of stable and unstable FPs in the 2D return map are eliminated sequentially (see Figs. 10, 11, 15, and 16) when the boundaries of these zones, corresponding to saddle-node (SN) bifurcations, are crossed outwardly as the coupling Δgij between neurons i and j is increased.

Close modal

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 circuits44 and their repertoire of multistable polyrhythms. Certain network topologies display a rich multiplicity (due to permutation-phase shift symmetries65) 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 

FIG. 31.

Larger network configurations: The techniques in this study may be extended to the analysis of larger network dynamics, including (a) the multistable fully connected, (b) the bistable one-way inhibitory loop, and (c) the robust monostable mixed 4-cell inhibitory circuits,44 as well as (d) a 6-cell network composed of two connected 3-cell motifs, reciprocally coupled with cross-inhibitory synapses (dashed gray).

FIG. 31.

Larger network configurations: The techniques in this study may be extended to the analysis of larger network dynamics, including (a) the multistable fully connected, (b) the bistable one-way inhibitory loop, and (c) the robust monostable mixed 4-cell inhibitory circuits,44 as well as (d) a 6-cell network composed of two connected 3-cell motifs, reciprocally coupled with cross-inhibitory synapses (dashed gray).

Close modal
FIG. 32.

Cons and pros of 3D visualization: unit-torus with a return map for three phase-lags between constituent cells in the 4-cell network depicted in Fig. 31(a). The visible FPs (blue and green dots) correspond to stable pacemaker rhythms.

FIG. 32.

Cons and pros of 3D visualization: unit-torus with a return map for three phase-lags between constituent cells in the 4-cell network depicted in Fig. 31(a). The visible FPs (blue and green dots) correspond to stable pacemaker rhythms.

Close modal
FIG. 33.

Voltage traces depicting some of the stable rhythms seen in larger inhibitory network configurations: 4-cell networks [see Figs. 31(a)31(c)] can produce four paired half-center rhythms in (a) or four traveling wave rhythms in (b), among many other outcomes.44 A modular 6-cell network [sketched in Fig. 31(d)], with 2 coupled 3-cell motifs, cells (1,2,3) and (4,5,6), can generate multiple pacemaker-like rhythms, depicted in panels (c) and (d), due to the reciprocal cross-inhibition.

FIG. 33.

Voltage traces depicting some of the stable rhythms seen in larger inhibitory network configurations: 4-cell networks [see Figs. 31(a)31(c)] can produce four paired half-center rhythms in (a) or four traveling wave rhythms in (b), among many other outcomes.44 A modular 6-cell network [sketched in Fig. 31(d)], with 2 coupled 3-cell motifs, cells (1,2,3) and (4,5,6), can generate multiple pacemaker-like rhythms, depicted in panels (c) and (d), due to the reciprocal cross-inhibition.

Close modal

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).

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

1.
J.
Wojcik
,
J.
Schwabedal
,
R.
Clewley
, and
A. L.
Shilnikov
, “
Key bifurcations of bursting polyrhythms in 3-cell central pattern generators
,”
PLoS One
9
,
e92918
(
2014
).
2.
J. P.
Miller
and
A. I.
Selverston
, “Neural mechanisms for the production of the lobster pyloric motor pattern,” in Model Neural Networks and Behavior (Springer, 1985), pp. 37–48.
3.
T.
Bal
,
F.
Nagy
, and
M.
Moulins
, “
The pyloric central pattern generator in Crustacea: A set of conditional neuronal oscillators
,”
J. Compar. Phys. A
163
,
715
727
(
1988
).
4.
E.
Marder
and
R. L.
Calabrese
, “
Principles of rhythmic motor pattern generation
,”
Physiol. Rev.
76
,
687
717
(
1996
).
5.
W. B.
Kristan, Jr.
,
R. L.
Calabrese
, and
W. O.
Friesen
, “
Neuronal control of leech behavior
,”
Prog. Neurobiol.
76
,
279
327
(
2005
).
6.
R. J.
Calin-Jageman
,
M. J.
Tunstall
,
B. D.
Mensh
,
P. S.
Katz
, and
W. N.
Frost
, “
Parameter space analysis suggests multi-site plasticity contributes to motor pattern initiation in Tritonia
,”
J. Neurophysiol.
98
,
2382
2398
(
2007
).
7.
W. E.
Sherwood
,
R.
Harris-Warrick
, and
J.
Guckenheimer
, “
Synaptic patterning of left-right alternation in a computational model of the rodent hindlimb central pattern generator
,”
J. Comput. Neurosci.
30
,
323
360
(
2011
).
8.
J. M.
Newcomb
,
A.
Sakurai
,
J. L.
Lillvis
,
C. A.
Gunaratne
, and
P. S.
Katz
, “
Homology and homoplasy of swimming behaviors and neural circuits in the Nudipleura (Mollusca, Gastropoda, Opisthobranchia)
,”
Proc. Natl. Acad. Sci. U.S.A.
109
,
10669
10676
(
2012
).
9.
R.
Milo
,
S.
Shen-Orr
,
S.
Itzkovitz
,
N.
Kashtan
,
D.
Chklovskii
, and
U.
Alon
, “
Network motifs: Simple building blocks of complex networks
,”
Science
298
,
824
827
(
2002
).
10.
O.
Sporns
and
R.
Kötter
, “
Motifs in brain networks
,”
PLoS Biol.
2
,
e369
(
2004
).
11.
M. I.
Rabinovich
,
P.
Varona
,
A. I.
Selverston
, and
H. D.
Abarbanel
, “
Dynamical principles in neuroscience
,”
Rev. Mod. Phys.
78
,
1213
(
2006
).
12.
A.
Bulloch
and
N.
Syed
, “
Reconstruction of neuronal networks in culture
,”
Trends Neurosci.
15
,
422
427
(
1992
).
13.
E.
Marder
, “
Invertebrate neurobiology: Polymorphic neural networks
,”
Curr. Biol.
4
,
752
754
(
1994
).
14.
W. N.
Frost
and
P. S.
Katz
, “
Single neuron control over a complex motor program
,”
Proc. Natl. Acad. Sci. U.S.A.
93
,
422
426
(
1996
).
15.
P. S.
Katz
, “
Evolution of central pattern generators and rhythmic behaviours
,”
Philos. Trans. R. Soc. B
371
,
20150057
(
2016
).
16.
D.
Alacam
and
A.
Shilnikov
, “
Making a swim central pattern generator out of latent parabolic bursters
,”
Int. J. Bifurcat. Chaos
25
,
1540003
(
2015
).
17.
T. G.
Brown
, “
The intrinsic factors in the act of progression in the mammal
,”
Proc. R. Soc. Lond. B
84
,
308
319
(
1911
).
18.
S.
Jalil
,
D.
Allen
,
J.
Youker
, and
A.
Shilnikov
, “
Toward robust phase-locking in Melibe swim central pattern generator models
,”
Chaos
23
,
046105
(
2013
).
19.
A.
Sakurai
and
P.
Katz
, “Distinct neural circuit architectures produce analogous rhythmic behaviors in related species,” in Society for Neuroscience Abstracts, 918.04 Neuroscience Meeting Planner, Washington, DC (Society for Neuroscience, 2011), Vol. 37
20.
A.
Sakurai
,
J. M.
Newcomb
,
J. L.
Lillvis
, and
P. S.
Katz
, “
Different roles for homologous interneurons in species exhibiting similar rhythmic behaviors
,”
Curr. Biol.
21
,
1036
1043
(
2011
).
21.
A.
Sakurai
,
C. A.
Gunaratne
, and
P. S.
Katz
, “
Two interconnected kernels of reciprocally inhibitory interneurons underlie alternating left-right swim motor pattern generation in the mollusk Melibe leonina
,”
J. Neurophysiol.
112
,
1317
1328
(
2014
).
22.
P. S.
Katz
, “
Comparison of extrinsic and intrinsic neuromodulation in two central pattern generator circuits in invertebrates
,”
Exp. Physiol.
83
,
281
292
(
1998
).
23.
A.
Sakurai
and
P. S.
Katz
, “
Artificial synaptic rewiring demonstrates that distinct neural circuit configurations underlie homologous behaviors
,”
Curr. Biol.
27
,
1721
1734
(
2017
).
24.
N.
Kopell
and
B.
Ermentrout
, “
Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks
,”
Proc. Natl. Acad. Sci. U.S.A.
101
,
15482
15487
(
2004
).
25.
K.
Matsuoka
, “
Mechanisms of frequency and pattern control in the neural rhythm generators
,”
Biol. Cybern.
56
,
345
353
(
1987
).
26.
N.
Kopell
, “Toward a theory of modelling generators,” in Neural Control of Rhythmic Movements in Vertebrates (Wiley, 1988).
27.
C.
Canavier
,
D.
Baxter
,
J.
Clark
, and
J.
Byrne
, “
Multiple modes of activity in a model neuron suggest a novel mechanism for the effects of neuromodulators
,”
J. Neurophysiol.
72
,
872
882
(
1994
).
28.
F.
Skinner
,
N.
Kopell
, and
E.
Marder
, “
Mechanisms for oscillation and frequency control in networks of mutually inhibitory relaxation oscillators
,”
J. Comput. Neurosci.
1
,
69
87
(
1994
).
29.
R.
Dror
,
C. C.
Canavier
,
R. J.
Butera
,
J. W.
Clark
, and
J. H.
Byrne
, “
A mathematical criterion based on phase response curves for stability in a ring of coupled oscillators
,”
Biol. Cybern.
80
,
11
23
(
1999
).
30.
A. A.
Prinz
,
C. P.
Billimoria
, and
E.
Marder
, “
Alternative to hand-tuning conductance-based models: Construction and analysis of databases of model neurons
,”
J. Neurophysiol.
90
,
3998
4015
(
2003
).
31.
A. L.
Shilnikov
and
G.
Cymbaluk
, “
Homoclinic bifurcations of periodic orbits en a route from tonic spiking to bursting in neuron models
,”
Regul. Chaotic Dyn.
9
,
281
297
(
2004
).
32.
A.
Shilnikov
, “
Complete dynamical analysis of a neuron model
,”
Nonlinear Dyn.
68
,
305
328
(
2012
).
33.
J.
Wojcik
,
R.
Clewley
, and
A.
Shilnikov
, “
Order parameter for bursting polyrhythms in multifunctional central pattern generators
,”
Phys. Rev. E
83
,
056209
(
2011
).
34.
A.
Shilnikov
,
R.
Gordon
, and
I.
Belykh
, “
Polyrhythmic synchronization in bursting networking motifs
,”
Chaos
18
,
037120
(
2008
).
35.
J. E.
Rubin
and
D.
Terman
, “
Explicit maps to predict activation order in multiphase rhythms of a coupled cell network
,”
J. Math. Neurosci.
2
,
4
(
2012
).
36.
R. L.
Calabrese
,
B. J.
Norris
,
A.
Wenning
, and
T. M.
Wright
, “
Coping with variability in small neuronal networks
,”
Integr. Comp. Biol.
51
,
845
(
2011
).
37.
W. B.
Kristan
, “
Neuronal decision-making circuits
,”
Curr. Biol.
18
,
R928
R932
(
2008
).
38.
K. L.
Briggman
and
W.
Kristan, Jr.
, “
Multifunctional pattern-generating circuits
,”
Annu. Rev. Neurosci.
31
,
271
294
(
2008
).
39.
J. T.
Schwabedal
,
D. E.
Knapper
, and
A. L.
Shilnikov
, “
Qualitative and quantitative stability analysis of penta-rhythmic circuits
,”
Nonlinearity
29
,
3647
(
2016
).
40.
T.
Xing
, “Computational study in chaotic dynamical systems and mechanisms for pattern generation in three-cell networks,” Dissertation (Georgia State University, 2015).
41.
X.-J.
Wang
and
J.
Rinzel
, “
Alternating and synchronous rhythms in reciprocally inhibitory model neurons
,”
Neural Comput.
4
,
84
97
(
1992
).
42.
D.
Somers
and
N.
Kopell
, “
Rapid synchronization through fast threshold modulation
,”
Biol. Cybern.
68
,
393
407
(
1993
).
43.
J.
Schwabedal
,
D.
Knapper
,
K.
Pusuluri
, and
D.
Alacam
(
2014
). “MotifToolbox,” GitHub.
https://github.com/jusjusjus/Motiftoolbox
44.