2 θ -burster for rhythm-generating circuits
22 θ -burster for rhythm-generating circuits Aaron Kelley
Neuroscience Institute, Georgia State University, Atlanta, GA, USA
Andrey Shilnikov
Neuroscience Institute and Department of Mathematics and Statistics,Georgia State University, Atlanta, USA (Dated: August 11, 2020)We propose and demonstrate the use of a minimal 2 θ model for endogenous bursters coupled in 3-cell neuralcircuits. This 2 θ model offers the benefit of simplicity of designing larger neural networks along with an acutereduction on the computation cost. I. INTRODUCTION
Neural networks called Central Pattern Generators (CPGs)[1–8]. produce and control a great variety of rhythmic mo-tor behaviors, including heartbeat, respiration, chewing, andlocomotion. Many anatomically and physiologically diverseCPG circuits involve a three-cell motifs including the spinylobster pyloric network [6], the
Tritonia swim circuit [4], andthe
Lymnaea respiratory CPGs [3]. Pairing experimental stud-ies and modeling studies have proven to be key to unlockinsights into operational and dynamical principles of CPGs[9 ? –14]. Although various circuits and models of specificCPGs have been developed, it still remains unclear how CPGsachieve the level of robustness and flexibility observed in na-ture. It is not clear either what mechanisms a single motor sys-tem can use to generate multiple rhythms, i.e., whether CPGsuse dedicated circuitry for each function, or whether the samecircuitry is multi-functional and can govern several behaviors[15–17].This paper capitalizes on our original work and the well-established principles in the characterization of 3-cell cir-cuits made of HH-type neurons [15, 18–20] and the Fitzhugh-Nagumo-like neurons [21]. We use the bottom-up approachto showcase the universality of rhythm-generation principlesin 3-cell circuits regardless of the model selected, which canbe a Hodgkin-Huxley (HH) type model of the leech heartinterneuron [22, 23], the the generalized Fitzhugh-Nagumo(gFN) model of neurons [24], and the minimal 2 θ burstingneuron, provided of course that all three models meet somesimple and generic criteria. II. RETURN MAPS FOR PHASE LAGS
We developed a computational toolkit for oscillatory net-works that reduces the problem of the occurrence of burstingand spiking rhythms generated by a CPG network to the bi-furcation analysis of attractors in the corresponding Poincaréreturn maps for the phase lags between oscillatory neurons.The structure of the phase space of the map is an individ-ual signature of the CPG as it discloses all characteristics ofthe functional space of the network. Recurrence of rhythmsgenerated by the CPG (represented by a system of coupledHodgkin-Huxley type neurons [23]) lets us employ Poincaré return maps defined for phase lags between spike/burst ini-tiations in the constituent neurons [25, 26] as illustrated inFig. 1,2 and 4. With such return maps, we can predict andidentify the set of robust outcomes in a CPG with mixed, in-hibitory/excretory and electrical synapses, which are differ-entiated by phase-locked or periodically varying lags corre-sponding, respectively, to stable fixed points and invariant cir-cles of the return map.Let us introduce a 3-cell network (Fig. 1A) made of weaklycoupled HH-like bursters; see the equations in the Appendixbelow. Here, “weakly” means that coupling cannot quitedisturb the shape of the stable bursting orbit in the 3Dphase space of the individual HH-model (Fig. 1A). Weak in-teractions, inhibitory (mainly repulsing) and excitatory/gap-junction (manly attracting) can only affect the phases of theperiodically varying states of the neurons, represented by thecolor-coded spheres, blue/green/red for cells 1/2/3, on thebursting orbit in the 3D phase space of the given interneu-ron model. As such weak-coupling can only gently alter thephase-differences or phase-lags between the coupled neurons(Fig. 2A). Being inspired by neuro-physiological recordingsperformed on various rhythmic CPGs, we employ only volt-age traces generated by such networks to examine the timedelays, τ and τ between the burst upstrokes on each cyclein the reference/blue cell 1 and in cells 2 (green) and 3 (red).In what follows, we will show that like the biologically plausi-ble HH-type networks, 3-cell circuits of coupled 2 θ -bursterscan stable produce similar phase-locked rhythms. They in-clude, but not limited, peristaltic patterns or traveling waves,in which the cells burst sequentially one after the other (seeFigs. 1 and 3C/E), as well as the so-called pacemaker rhythms,in which one cell effectively inhibits and bursts in anti-phasewith the other two bursting synchronously (Fig. 3B/D). Thesymmetric connectivity implies such 3-cell networks can pro-duce multiple rhythms due to cyclic permutations of the con-stituent cells (see Fig. 3 below). To analyze the existence andthe stability of various recurrent rhythms produced by suchnetworks, we employ our previously developed approach us-ing Poincaré return maps for phase-legs between constituentneurons. We introduce phase-lags defined at specific events intime when the voltage in cells reaches some threshold valuethis signaling the burst initiation (see Fig. 1B). The phase lag ∆ φ ( n ) j is then defined by a delay between n -th burst initiationsin in the given cell and the reference cell 1, normalized over a r X i v : . [ q - b i o . N C ] J u l FIG. 1. (A) Snapshots of the current states (represented by the blue, green and red spheres) of three weakly-coupled endogenous bursters at t = t =
10, on the periodic orbit (grey) in the 3D phase space of the Hodgkin-Huxley type model of the leechheart interneuron [22, 23]. The plane V = Θ syn represents the threshold for the chemical synapses, that divides the active “on” phase (aboveit) and the inactive “off” phase; here, the active red cell inhibits the quiescent green and blue ones. (B) Burst initiations in successive voltagetraces generated by the 3-cell neural network allow us to define the relative delays τ i ’s and hence, the phase lags (given by Eqs. (1)) betweenits constituent bursters; see further details in [25, 26]. the bursting period: ∆ φ ( n ) = t ( n ) − t ( n ) t ( n + ) − t ( n ) , ∆ φ ( n ) = t ( n ) − t ( n ) t ( n + ) − t ( n ) . mod 1 , (1)The sequence of phase lags { ∆ φ ( n ) , ∆ φ ( n ) } , defined forvalues between 0 and 1, gives the forward phase trajectoryon a 2D torus (Fig. 2B). The specific phase-lag values 0(or 1) and 0 . ( ∆ φ , ∆ φ ) -phase space of the 2D Poincaré return maps(such as one shown in Fig. 3A) of the 3-cell networks by ini-tiating multiple trajectories with a dense distribution of initialphase-lags (50 ×
50 grid), and by following their progressionsover large numbers of cycles. On long runs these trajecto-ries can eventually converge to some attractors, one or several.
FIG. 2. (A) Slow exponential convergence of initial states of ∆ φ (yellow curves) and ∆ φ (purple curves) to four phase-locked states: { ≡ , , , } , in the inhibitory 3-cell motif (4) with weak coupling β = 0.003. (B) Poincaré return map defined on a unit 2D torus, T = S ⊗ S of two phase-lags, showing color-coded attraction basins ofseveral fixed points (solid dots of same colors) corresponding to thephase-locked rhythms by the 3-cell motif. A flatten torus is shown inFig. 3A. Such an attractor can be a fixed point (FP) with constant val-ues ∆ φ ∗ and ∆ φ ∗ in (1)), which correspond to a stable rhyth-mic pattern with phase-lags locked (Fig. 2A). All phase tra-jectories converging to the same fixed point are marked by thesame color to reveal the attraction basins of the correspondingrhythms. This reduces the analysis of rhythmic activity gener-ated by a 3-cell network to the examination of the correspond-ing 2D Poincaré map for the phase-legs. For example, the mapshown Fig. 3A. reveals the existence of penta-stability in thesymmetric circuit generating three pacemakers (blue, greenand red) and two, clockwise and counter-clockwise, travel-ing waves (Fig. 3B). These three PM rhythms correspond tothe blue, green and red FPs located at ( . , . ) , ( . , ) and ( , . ) , respectively, while two traveling wave pattern are as-sociated with stable FPs located at ( / , / ) and ( / , / ) ,respectively, in the 2D return map. Other type of attractorscan be a stable invariant curve corresponding to rhythmic pat-tern wit (a)periodically varying phase-lags. Such a curve canbe a circle on and wrap around the 2D torus (see Figs. 2A and3A). If the map has a single attractor, then the correspond-ing network is mono-stable, otherwise it is a multifunctionalor multistable network capable of producing several rhythmicoutcomes robustly. The 2D return map: M n → M n + , for thephase-lags can be represented as follows: ∆ φ ( n + ) = ∆ φ ( n ) + µ f ( ∆ φ ( n ) , ∆ φ ( n ) ) , ∆ φ ( n + ) = ∆ φ ( n ) + µ f ( ∆ φ ( n ) , ∆ φ ( n ) ) (2)with small µ i being associated with weak coupling; f i aresome undetermined coupling functions such that their zeros: f = f = ∆ φ ∗ j = ∆ φ ( n + ) j = ∆ φ ( n ) j of the map. These functions, similar to phase-resettingcurves, can be numerically evaluated from the simulated dataon all trajectories { ∆ φ ( n ) , ∆ φ ( n ) } (see Fig. 4C). By treating f i as partials ∂ F / ∂ φ i j , one may try to restore a “phase poten-tial” – some surface F ( φ , φ ) = C (see Fig. 4). The shapeof such a surface defines the location of critical points asso- FIG. 3. Multistable outputs of the 3-cell homogeneous network withsix equal synaptic connections ( β = . ( ∆ φ , ∆ φ ) -phase lags with five stable fixed points rep-resenting robust three pacemaker (PM) patterns: red at ( , ) , greenat ( , ) and blue at ( , ) , and two traveling wave (TW) rhythmicpatterns: yellow clockwise at ( , ) and teal counter-clockwise at ( , ) . The color-coded attraction basins of these five FPs are deter-mined by positions of stable sets (separatrices) of six saddles (graydots). The origin is a repelling FP of the map with the even number –total eight of hyperbolic FPs in the map. Panels B-E depict the traceswith phases locked to the specific values (indicated by color-codeddots at top-left corners), corresponding to the selected FPs. ciated with FPs – attractors, repellers and saddles of the map.With this approach one can try to predict bifurcations due tolandscape transformations and therefore to interpret possibledynamics of the network as a whole. Figure 4A and B aremeant to give an idea how the potential surface may look likein the case of the 3-cell circuit with only two stable travelingwave patterns and in the case of three co-existing pacemakersonly, respectively. III. MINIMALISTIC 2 θ -BURSTER The concept of the 2 θ -burster is inspired by the dynamicsof endogenous bursters with two characteristic slow phases:depolarized tonic-spiking and quiescent, like one shown inFig. 1. These phases are often referred to as “on” or activeand “off” or inactive depending on whether the membranevoltage is above or below the synaptic threshold. During theactive phase the pre-synaptic cell releases neurotransmittersto inhibit or excite other cells on the network, while duringthe inactive phase, the cell does not “communicate” to any-one. This is a feature of chimerical synapses unlike the elec-tric synapses that let cells interact all the times regardless ofthe voltage values. The predecessor of the 2 θ -burster is theso-called “spiking” θ -neuron [34]. Mathematically, it is anormal form for the plain saddle-node bifurcation on a cir-cle through which two equilibrium state, stable and repelling,merge and disappear. After the phase point keeps traverse thecircle. That is why this bifurcation is referred to as a ho- FIG. 4. Critical points of the sketched “pseudo-potentials” with pe-riodic boundary conditions reveal the location of potential dwells– attractors, as well saddles (including one with six separatrices in(B)) and repellers in the ( φ , φ )-phase surface. These configu-rations correspond the network with only two traveling waves andwith only 3 pacemakers. (C) A computational reconstruction of apseudo-potential/coupling function corresponding to the return mapin Fig. 3A. moclinic Saddle-Node bifurcation on an Invariant Circle, orSNIC for short. The θ -neuron capitalizes on the pivotal prop-erty of the saddle-node bifurcation – the phantom bottle-neckeffect that underlies slow and fast time scales in the dynamicsof such systems, see Fig. 5A. Recall that a similar saddle-node bifurcation, which controls the tonic-spiking phase andthe number of spikes in bursters, is associated with the blue-sky catastrophe [23, 27–29].The key feature of the 2 θ -neuron given by θ ′ = ω − cos 2 θ + α cos θ , mod 1 (3)is the occurrence of two saddle-node bifurcations, thatintroduce two slow phases into its dynamics, with two fasttransition in between, see Fig. 5B. Similar to endogenousbursters with two slow transient states – the active tonic-spiking and the quiescent phases that can be controlledindependently, we can manage the durations of the twoanalogous states in the 2 θ -neuron: “on” at π and “off” at 0,using the same bottleneck post-effects of the two saddle-nodebifurcations. This allows us to regulate its duty cycle, whichis the fraction of the active-state duration compared to theburst period, see Fig. 5B. As seen from Fig. 5, the θ -modelwas meant to replicate phenomenologically fast spiking cells,while the 2 θ -neuron can be treated as a “spike-less” burster.Below, we demonstrate that the network dynamics producedby a 3-cell motif composed of inhibitory 2 θ -bursters, pre-serve all the key features seen in a motif composed of thethree Hodgkin-Huxley-type bursters (see Fig. 1).First, let us observe from Eq. 3 that the phase dynamics FIG. 5. Comparison of the oscillatory dynamics generated by thespiking θ -neuron and the 2 θ -burster. Panels A and C present snap-shots of typical trajectories generated by both models on a unit circle S (parametrized using Cartesian coordinates: x ( t ) = sin ( θ ( t )) and y ( t ) = − cos θ ( t ) ) with the origin 0 at 6pm. (A) Clustering of pur-ple spheres near the origin is due to a bottleneck post-effect causedby a saddle-node bifurcation (SNIC) in the θ model, while the 2 θ -burster in (C) features two such bottleneck post-effects due to twoheteroclinic saddle-node connections causing the stagnation of grayspheres near the top, “on” state and the inactive “off” state of the2 θ -burster and fast transitions in between. (B) Spiking trace (purple)of the θ -neuron, being overlapped with 2-plateau traces of the 2 θ -neuron with three values of the duty cycles ≃ of an individual 2 θ -burster is mainly governed by the terms ω − cos 2 θ . As long as the frequency 0 < ω ≤
1, there are twopair of stable and unstable equilibria: one pair is at the bot-tom around θ ≃
0, while the other is at the top near θ ≃ π .The stable states are associated with the hyperpolarized andthe depolarized quiescent states of the neuron. When ω > θ -burster becomes oscillatory through two simultaneous(provided α =
0) saddle-node bifurcations (SNIC) on a unitcircle S , where θ is defined on modulo 1. Moreover, aslonger as ω = + ∆ ω , where 0 < ∆ ≪
1, this new bursterpossesses two slow phases: the active “on” state near θ = π , and the inactive “off” state near 0 on S . These slow statedare alternated with fast counter-clockwise transitions, whichare referred to as an upstroke and a downstroke, respectively.For greater values of ω , the active and inactive phases are de-fined bore broadly: π / < θ ≤ π / π / < θ ≤ π / θ th = π / θ th = θ -burster iscontrolled by the term α cos θ , provided that it remains oscil-latory as long as ω − ∣ α ∣ >
1. Note that when α =
0, the dutycycle is 50% and the burster oscillations have even plagues
FIG. 6. (A) Sampling the moments in phase traces, y i ( t ) = − cos ( θ i ( t )) , plotted against time, when they reach a synaptic thresh-old θ syn =
0, defines a sequence of the phase lags ( τ ( n ) , τ ( n ) ) be-tween upstrokes in the reference, blue neuron and other 2 θ -neuronscoupled in the 3-cell network. (B) Parametric representation of the1D phase space of coupled 2 θ -bursters traversing counter-clockwise(long gray arrows indicating rapid transition between on-off states)on a unit circle S . Small-downward blue and red arrows illustrat-ing the inhibition perturbations from the active green cell above thesynaptic threshold that delays the forthcoming upstroke of the bluecell, and speeds up the red cell toward the inactive phase. (see Fig. 5B). The active or inactive phases can be extendedor shortened, respectively, with α <
0, making the duty cyclegreater, or vice versa – the duty cycle of individual burster isdecreased with α > IV. 3 EQUATIONS FOR 3-CELL NETWORK
A 3-cell circuit of the 2 θ -bursters coupled with chemicalsynapses is given by the following system: ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ θ ′ = ω − cos 2 θ + α cos θ − [ β + e k cos θ + β + e k cos θ ] ⋅ [ − + e k sin θ ] , θ ′ = ω − cos 2 θ + α cos θ − [ β + e k cos θ + β + e k cos θ ] ⋅ [ − + e k sin θ ] , θ ′ = ω − cos 2 θ + α cos θ − [ β + e k cos θ + β + e k cos θ ] ⋅ [ − + e k sin θ ] , mod 1. (4)The 2 θ -burster are coupled in the network using the fast inhibitory synapses driven by the fast-threshold modulation[33]. It is due to the positive “sigmoidal” term [ + e k cos θ i ] that, rapidly ((here k =
10) varying between 0 and 1, triggersan influx of inhibition flowing from the pre-synaptic neuroninto the post-synaptic neuron, as soon as the former enters theactive on-phase above the synaptic threshold cos θ th =
0, i.e., π / < θ i < π /
2. Note that the negative sign of this termmakes the synapse inhibitory; replacing it with “+” makes thesynapse excitatory because it would increase the rate of θ ′ during the upstroke, contrarily to slowing the upstroke downas in the inhibitory case. The strength of the coupling is de-termined by the maximal conductance values β i j .The last term [ − + e k sin θ ] , breaking the symmetry, con-verts the synaptic input into qualitative inhibition. Namely, itssign is switched from + to - upon crossing the values θ = θ = π . During the fast upstroke, when 0 < θ < π , the thisterm is positive, thereby ensuring that inhibition does slowdown or delay the transition into bursting. When π < θ i < π during the fast downstroke, this terms [ − + e k sin θ ] < [ − + e k sin θ ] that breaks the symmetry waswell as it acts during the upstroke only.The electrical coupling or the gap junction between the neu-rons is handled by the another term − C elec sin ( θ pre − θ post ) .It slows down the rate θ ′ post when θ post > θ pre and speeds itup if θ post < θ pre . The conductivity coefficient C elec is to beset around two orders of magnitude smaller than β -values tomaintain a balanced effect in the network. When C elec and β are of the same magnitude, the dynamics of network aresolely dictated by the electrical coupling with the inhibitorysynapses insignificantly affecting it. V. POINCARÉ RETURN MAPS FOR THE PHASE-LAGS.RESULTS
Figure 6 shows the snapshots of the phase progressions ofthe three coupled 2 θ -bursters on the unit circle S and depictshow phase-lags between the are introduced (here, we refer tocell 1 (blue) as the reference one). One can see from thisfigure that the active green neuron 2 in the active phase closeto θ = π , above the synaptic threshold, inhibits and pushes theother two closer to each other, near the bottom quiescent stateat θ =
0, by accelerating the red burster 3 on the downstroke,and by slowing down the blue burster 1 on the upstroke.Following the same approach used in the weakly coupledHH-type models above, we first create a uniform distribu-tion of initial phases on S , and therefore the phase-lags be-tween the three 2 θ -bursters. Next we integrate the network (4)over a large number of cycles, and record burst initiations(see Fig. 5A) to determine the phase-lags between the ref-erence cell 1 and two other cells and to what phase locked states they can converge with increasing number of the cy-cles. This approach is illustrated in Fig. 2A for the symmet-ric 3-cell motif composed of identical 2 θ -bursters and equalinhibitory synapses. The corresponding 2D Poincaré returnmap, with the co-existing stable fixed points and saddles isshown in Figs. 3. By stitching together the opposite sides ofthis map, we wrap it around a 2D torus as shown in Fig. 2B.The fixed points and their attraction basins are coded withdifferent colors in the map. For example, the Poincaré returnmap for the ( ∆ φ , ∆ φ ) -phase lags represented in Fig. 3Ahas five stable fixed points representing robust three pace-maker FPs located at: red ( , ) , green at ( , ) and blueat ( , ) , and two traveling-wave ones: yellow clockwise at ( , ) and teal counter-clockwise at ( , ) . The borders ofthe attraction basins of these five FPs are determined by po-sitions of stable sets (separatrices) of six saddles (gray dots).The origin is a repelling FP of the map. It totals up to eighthyperbolic FPs in the map.To conclude this section, we would like to point out anotherhelpful feature of the 2 θ -neuron paradigm, namely, the abil-ity to find repelling FPs, if any, in the 2D Poincaré map, byreversing the direction of integration of the system (4), i.e.,integrating it in the backward time after multiplying its right-hand sides by -1. Unlike HH-type dissipative systems wherethe backward integration will make solutions run to infinity, itis not the case for the 2 θ -bursters, as just reverses the directionand spin trajectories clockwise on S . A. Symmetric Motif
It will be shown below that the 2 θ -bursters weakly coupledin the 3-cell networks, symmetric, asymmetric and with mixedsynapses, can generate the same stable rhythms as the net-works of biologically plausible HH-type models. We will alsodiscuss the bifurcations occurring in the networks and corre-sponding maps as the synaptic connectivity or intrinsic tem-poral characteristics of the 2 θ -bursters are changed. Bifurca-tions in the system are identified and classified by a change ofthe stable phase rhythms which can be due to the stability lossof a particular FP, or when it merges with a close saddle soboth disappear through a saddle-node bifurcation.Let us first consider a symmetric network with two bifur-cation parameters: the coupling strength β = β ii ( i = , , ) and the α -parameter in Eq. (3) that controls the duty cycle(DC) of the 2 θ -bursters. We use five different DC-values as α is varied from -0.11 to 0.11l while synaptic strength is in-creased through four steps from β = 0.0001 through 0.1. Theresults are presented in Fig. 7. The Panels A2/3 represent themost balanced, weakly coupled network that can produce allfive bursting rhythms with the DC 50%. One can see thatincreasing the β -value, the saddles separating 2 TWs and 3PMs move toward the latter ones, and over some critical val-ues, 3 pairs: a saddle and the nearest stable PM merger andvanish simultaneously. After that, the symmetric network canproduce two only rhythms: counter- and clockwise TWs cor-responding to the teal and yellow stable FPs at ( , ) and FIG. 7. Bifurcations of FPs in the ( ∆ φ , ∆ φ ) -return map for the symmetric motif as the coupling β -parameter and the duty cycle (viavariations of α ) are changed; parameters: β -values are [0.001 ,0.003, 0.01, 0.03] from top to bottom labeled A to D, resp., while α -values are[-0.11,-0.05,0.0,0.11] from left to right labeled, 1 through 4, respectively, with 50% DC at α = . β -values, therate of convergence to the FPs increases. The TW-rhythms dominate the network dynamics when the DC is about 50%, as seen in the middlecolumns. The PM-rhythms become dominant at small and large DC-values, as depicted in the outer panels. ( , ) , respectively. This would correspond to the case of the“pseudo-potential” depicted in Fig. 4A.The stable PMs are promoted or dominate the dynamics ofthe symmetric at the extreme α -values corresponding to thebursting rhythms with short or long burst durations. Once cancompare panels, say A1 and D4 reveal that this time, the sepa-rating saddles group around the stable TWs to minimize theirattraction basins, and hence the likelihood of the occurrenceof these rhythms in the network. These case would correspondto the “pseudo-potential” depicted in Fig. 4B. B. “King of the mountain” motif
The first asymmetric case considered is a motif termed theKing of the Mountain. In this modeling scenario both out- going inhibitory synapses from the given cell, here the refer-ence blue burster 1 one, are evenly increased in the strength,see Fig. 8A. Observe that such a configuration breaks downboth circular symmetries supporting traveling waves in thenetwork. Let us start with Fig. 8B: no surprise that with initialincrease in β , / , two saddles shift away from the blue PMat (0.5,0.5) toward two TWs, then merge with them to dis-appear pair-wisely. Next, as β , / is increased further, twoother saddles annihilate the green and red PMs through sim-ilar saddle-node bifurcations (Fig. 8C). At the aftermath, the3-cell network with a single burster generating the repulsiveinhibition much stronger than the other two cells becomes amonostable one producing a single pacemaking rhythm withthe phase-lags locked at (0.5,0.5). FIG. 8. (A) “King of the mountain” network motif with twosynapse strengths, β and β , increased (indicated by darker con-nections), relative to the other synapse strengths. (B) The first ofthree ( ∆ φ , ∆ φ ) return maps, with β and β synaptic strengthsslightly greater than the other β s, the (blue) attraction area extendsso that the two saddles nearest the blue PM at ( , ) , move awayfrom the blue PM, closer towards the yellow and teal TWs at ( , ) and ( , ) , respectively. (C) With further increase of β and β ,these saddles and TWs merge with and annihilate each other throughsaddle-node bifurcations, and the blue PM basin grows. (D) Atgreater β and β values, the network becomes a winner-take-all, blue PM winning, after the red and green PMs, at ( , ) and ( , ) , respectively, vanish through subsequent saddle-node bifurca-tions. The parameters are: ω = 1.15, α = 0.07, and β = 0.003 except β and β = 0.0038, 0.004, 0.015 for panels B-D. C. Mono-biased motif
We refer as a mono-biased motif to another asymmetricthe network with a single different synapse: in this case thestrength β of the outgoing synapse from cell 2 to cell 1 is in-creased, which violates the circular symmetry supporting thecounter-clockwise traveling wave in the network, see Fig. 9F.So, as β is increased the counter-clockwise stable FP at ( , ) first disappears through a saddle-node bifurcation, asseen in Fig. 9A/B. Because this was the saddle between thisTW and he green PM, then the attraction basin of the lat-ter increases after the first bifurcation in the sequence. Thenext saddle-node bifurcation eliminates the red stable FP at(0, 0.5). The reasoning is the following: for this rhythm to per-sist the red PM is to evenly inhibit both green and blue PMs.However, a growing inhibition misbalance between them is nolonger reciprocal. As we pointed out earlier, the stronger inhi-bition from cell 2 shortens the active phase of the blue burster.As so they cannot be longer lined up by the burster 3, whichcauses the disappearance of this PM-rhythm and the FP itself(Fig. 9C). Same arguments can be just to justify the the dis-appearance of the green PM as cell 2 cannot not even inhibitcells 1 and 2 to hold them together as β is increased further(not shown). This is in the his case is in good agreement with FIG. 9. Mono-biased network motif (F) with one different synapsedue to increasing β . (A) The first of five ( ∆ φ , ∆ φ ) return maps,an increase in β value breaks down a counter-clockwise symmetryso that the attraction basin (teal) of the corresponding TW at ( , ) shrinks as a nearby saddle moves closer to it and away from the greenPM at ( , ) (A and B). (C) With further increase of β , the counter-clockwise TW at ( , ) vanishes through a saddle-node bifurcationafter merging with the nearby saddle, followed by another saddle-node bifurcation eliminating the red PM at (0, 0.5) (D). At greater β values the green PM ( , ) encompasses the majority of the networkphase space, along with the blue PM at ( , ) preserving the size ofits attraction basin. The parameters are: ω = 1.15, α = 0.07, and β ’s= 0.003 except β = 0.00042, 0.0045, 0.01, 0.02 for panels A-D. the 3-cell networks of the HH-type bursters. D. Dedicated HCO
The abbreviation HCP stands for a half-center oscillator,which a pair of neurons coupled reciprocally by inhibitorysynapses to produce alternating bursting. Such a dedicatedHCO is formed by cells 2 and 3 with stronger synapses due to β = β in the configuration shown in Fig. 10C. Again withstart off with the symmetric case depicted in Fig. 10A. Onecan observe at once, that having the dedicated HCO shouldbreaks down the circular symmetries of the network. So, thestable TWs become eliminated first as β = beta starts in-creasing. As these synapses become stronger the attractionbasin of the blue PM at (0.5 0.5) shrinks substantially, but theFP itself persists. Meanwhile increasing β = beta furthercreates the inhibitory misbalance that males the further exis-tence of the green and red PMs impossible due to the factorsthat we outlines above for the mono-biased motif. Both van-ish at the same time due to saddle-node bifurcations. How-ever, at the bifurcation both double FPs are connected by aheteroclinic orbit that transforms into a stable invariant curvewrapping around the phase torus (Figs 10F and 2B). This sta-ble invariant curve is associated with a phase-slipping rhythmsthat recurrently passes slowly through the “ghosts” of all fourvanished FPs except for the coexisting blue PM, see the frag- FIG. 10. (C) “Pairwise-biased” network motif with two recipro-cal synapse strengths β and β , increased. (A) The first of five ( ∆ φ , ∆ φ ) return maps, with β and β slightly greater thanother synaptic connections the network possesses all five attractingFPs. (B) Evenly increasing β and β values breaks down the ro-tational symmetry of the network so that both TWs at ( , ) and ( , ) vanish through saddle-node bifurcations while that the redand green PM basins equally expand and the blue basin shrinks.Here, two areas of the map, due to slow transitions throughout thesaddle-node ghosts, are color-coded in black because of uncertaintyin ultimate convergence/destination. (D-E) With further increasesof β , β values, the blue basin continues to shrink until red andgreen basins encompass almost all of the areas of the map. One cansee from Panel E that that the red and green PMs at ( , ) and ( , ) are also about to merge with nearby saddles and disappear throughtwo homoclinic saddle-node bifurcations (SNIC). (F) At greater val-ues of β , β , the blue PM at ( , ) has only a very narrow at-traction basin, corresponding to the only phase-locked rhythm, co-exists with a dominant phase-slipping repetitive pattern. The phaseslipping (its trace shown in Panel F1) corresponds to a stable invari-ant curve, passing throughout ( , ) and wrapping abound the 2Dtoroidal phase space to re-emerge near ( , ) and so forth. (F1) Fiveexemplary episodes of the traces vs. time showing periodically vary-ing phase lags between three cells (slipping). The parameters are: ω = 1.15, α = 0.07, and β = 0.003, except β and β are 0.005, 0.006,0.009, 0.035, in panels A, B, D-F and F1 is equal to F. ments of traces in Fig. 10F. E. Clockwise-biased motif
The clockwise-biased motif in this case represents the 3-cell network canter-clockwise connections stronger than onesin the opposite direction, see Fig. 11E. This configuration doesnot break circular symmetries of the network but infers that ei-ther TW should gain over the opposite one, which should re-sult in that their attraction basins should change correspond-ingly. Figure 11 presents four transformation stages of the
FIG. 11. (E) Clockwise-biased motif with three synaptic strengths, β , β and β sequentially increased. (A) As all three counter-clockwise synapses are slightly strengthen, saddles shift away fromthe three stable PMs, blue at ( , ) , green ( , ) and red ( , ) , to-wards the teal clockwise TW at ( , ) (B) thus shrinking its basinand widening the attraction basin of the dominant counter-clockwiseTW (yellow) at ( , ) (C). (D) With the stronger synaptic values, thethree saddles collapse into the CC TW, which becomes a complexsaddle with three incoming and three outgoing separatrices. The pa-rameters are ω = 1.15, α = 0.07, β = 0.003 except β , β and β = 0.0033, 0.025, 0.035, 0.055 for panels A-D. map as β , β and β sequentially increased. With a smallincrease, the shape of the map becomes a bit twisted with thethree saddles shifting away from the stable PMs toward theteal TW at ( , ) . The further increasing brings the saddleclose to the latter one thereby shrinking its attraction basinand substantially widening the basin of the clockwise TW at ( , ) . Finally, as some bifurcation threshold is reached, thesaddles collapse at the stable FP that becomes a complex sad-dle with three outgoing and three incoming separatrices. Thismeans that the counter-clockwise TW becomes an unstablerhythm in such biased 3-cell motif that is fully dominated bythe clockwise TW rhythm. F. Gap junction
In out last example we consider the symmetric motif witha gap junction or an electric synapses added between cells1 and 2 as shown in Fig. 12C. Recall that a gap junction isbi-directional unlike uni-directional chemical synapses withsynaptic thresholds. Recall that it is modeled by this term − C elec sin ( θ pre − θ post ) that slows down the rate θ ′ post when θ post > θ pre and speeds it up if θ post < θ pre . Due to this prop-erty, the electrical like excitatory synapse promote synchronybetween such coupled oscillatory cells, which in our case be-tween cells 1 and 2.Observe that introducing an electrical synapse betweenonly two of the cells of the motif ruins both circular sym-metries in the system. This is documented in Fig. 12A/B de- FIG. 12. Gap junction in the symmetric 3-cell network (C) is rep-resented by a resistor symbol placed between cells 1 and 2. (A) At C elec = . C elec breaks the circular symmetries of the network thatmakes both TWs at ( , ) and ( , ) vanish through saddle-nodebifurcations while the basin of the red PM at ( , ) widens. (D)With an even greater electrical coupling the red PM becomes thewinner-takes-all after the electrical connection ensures the in-phasesynchrony between cells 1 and 2 (C) that eliminates the blue andgreen PMs in the map after subsequent saddle-node bifurcation. Theparameters are: ω = 1.15, α = 0.07, β = 0.003, and C elec = . picting the maps for the networks with C elec being increasedfrom zero to 0 . C elec makes the stable green and blue stable PMs disintegrateas both cells become synchronous to burst in alternation withthe red cell 3. This completes the consideration of the mono-stable network with a relatively strong gap junction betweencells 1 and 2 that can only produce the only one pacemakerrhythm. VI. DISCUSSION
The goal of this paper is to demonstrate the simplic-ity and usability of the 2 θ -bursters to construct multistable,polyrhythmic neural networks that have the same dynami-cal and bifurcation properties as ones composed of biologi-cally plausible models of Hodgkin-Huxley type bursters andsynapses. Our de-facto approach is based on the computa-tional reduction to the visually evident Poincaré return mapsfor phase lags derived from multiple voltage traces. Thesemaps serve as a detailed blueprint containing all necessary in-formation about the network in questions, including its rhyth-mic repertoire, stability of generated patterns, etc, and in ad-dition to ability to predict possible transformations before thatoccur in the system. Our greater goal is to gain insight into thefundamental and universal rules governing pattern formation in complex networks of neurons. We believe that one shouldfirst investigate the rules underlying the emergence of coop-erative rhythms in basic neural motifs, as well as the role ofcoupling and in generating a multiplicity of coexisting rhyth-mic outcomes [35]. VII. ACKNOWLEDGEMENTS
We thank the Brains and Behavior initiative of GeorgiaState University for the pilot grant support. The authorsthank the past and current lab-mates of the Shilnikov NeurDSlab (Neuro–Dynamical Systems), specifically, K. Wojcik, K.Pusuluri, J. Collens, and J.Schwabedal, The NeurDS lab isgrateful to NVIDIA Corporation for donating the Tesla K40GPUs that were actively used in this study. This paper wasfunded in part by the NSF grant IOS-1455527.
VIII. APPENDIX
The time evolution of the membrane potential, V , of eachneuron is modeled using the framework of the Hodgkin-Huxley formalism, based on a reduction of a leech heart in-terneuron model, see [23] and the references therein: CV ′ = − I Na − I K2 − I L − I app − I syn , τ Na h ′ Na = h ∞ Na ( V ) − h , τ K2 m ′ K2 = m ∞ K2 ( V ) − m K2 . (5)The dynamics of the above model involve a fast sodium cur-rent, I Na with the activation described by the voltage depen-dent gating variables, m Na and h Na , a slow potassium current I K with the inactivation from m K2 , and an ohmic leak current, I leak : I Na = ¯ g Na m h Na ( V − E Na ) , I K2 = ¯ g K2 m ( V − E K ) , I L = ¯ g L ( V − E L ) . (6) C = . I app = . K2 = g Na = L = Na = K = − E L = − τ K2 = .
9s and τ Na = . h ∞ Na ( V ) , m ∞ Na ( V ) , m ∞ K2 ( V ) , of the of gating variables are determined by the fol-lowing Boltzmann equations: h ∞ Na ( V ) = [ + exp ( ( V + . ))] − m ∞ Na ( V ) = [ + exp ( − ( V + . ))] − m ∞ K2 ( V ) = [ + exp ( − ( V + . + V shiftK2 ))] − . (7)Fast, non-delayed synaptic currents in this study are modeledusing the fast threshold modulation (FTM) paradigm as fol-lows [33]: I syn = g syn ( V post − E syn ) Γ ( V pre − Θ syn ) , Γ ( V pre − Θ syn ) = /[ + exp { − ( V pre − Θ syn )}] ; (8)0here V post and V pre are voltages of the post- and the pre-synaptic cells; the synaptic threshold Θ syn = − .
03V is cho-sen so that every spike within a burst in the pre-synaptic cellcrosses Θ syn , see Fig. 1. This implies that the synaptic cur-rent, I syn , is initiated as soon as V pre exceeds the synapticthreshold. The type, inhibitory or excitatory, of the FTMsynapse is determined by the level of the reversal potential, E syn , in the post-synaptic cell. In the inhibitory case, it is setas E syn = − . V post ( t ) > E syn . In the excitatorycase the level of E syn is raised to zero to guarantee that theaverage of V post ( t ) over the burst period remains below the re-versal potential. We point out that alternative synapse models,such as the alpha and other detailed dynamical representation,do not essentially change the dynamical interactions betweenthese cells [19]. [1] Marder, E. and R.L. Calabrese (1996). Principles of rhythmicmotor pattern generation. Physiol Rev.
Prog. Neurobiol. , 76:279, 2005.[3] E. Marder. Invertebrate neurobiology: Polymorphic neural net-works.
Current Biology , 4(8):752 – 754, 1994.[4] RJ Calin-Jageman, MJ Tunstall, BD Mensh, PS Katz, and FrostWN. Parameter space analysis suggests multi-site plasticitycontributes to motor pattern initiation in tritonia.
J. Neurophys-iol. , 98:2382, 2007.[5] JM Newcomb, A Sakurai, JL Lillvis, CA Gunaratne, andPS Katz. Homology and homoplasy of swimming behaviorsand neural circuits in the nudipleura (mollusca, gastropoda,opistho-branchia).
Proc. Natl. Acad. Sci. , 109(1):10669–76,2012.[6] A Selverston, editor.
Model Neural Networks and Behavior .Springer, Berlin, 1985.[7] T. Bal, F. Nagy, and M. Moulins. The pyloric central patterngenerator in crustacea: a set of conditional neural oscillators.
Journal of Comparative Physiology A , 163(6):715–727, 1988.[8] PS Katz and SL Hooper. Invertebrate central pattern generators.In G North and R R. Greenspan, editors,
Invertebrate Neurobi-ology . Cold Spring Harbor Laboratory Press, NY, New York,2007.[9] N. Kopell and B. Ermentrout. Chemical and electrical synapsesperform complementary roles in the synchronization of in-terneuronal networks.
Proc. Natl. Acad, Sci. , 101(43):15482–15487, 2004.[10] K Matsuoka. Mechanisms of frequency and pattern control inthe neural rhythms generators.
Biol. Cybernetics , 1:1, 1987.[11] N. Kopell. Toward a theory of modeling central pattern gen-erators. In A.H. Cohen, S. Rossingol, and S. Grillner, editors,
Neural Control of Rhythmic Movements in Vertebrates . Wiley,New York, 1988.[12] C. C. Canavier, D. A. Baxter, J. W. Clark, and J. H. Byrne.Multiple modes of activity in a model neuron suggest a novelmechanism for the effects of neuromodulators.
J Neurophysiol ,72(2):872–882, Aug 1994.[13] F Skinner, N Kopell, and E Marder. Mechanisms for oscilla-tion and frequency control in networks of mutually inhibitoryrelaxation oscillators.
Comput. Neurosci. , 1:69, 1994.[14] R. O. Dror, C. C. Canavier, R. J. Butera, J. W. Clark, andJ. H. Byrne. A mathematical criterion based on phase responsecurves for stability in a ring of coupled oscillators.
Biol Cybern ,80(1):11–23, Jan 1999.[15] A.L. Shilnikov, R. Gordon, and I. Belykh. Polyrhyth-mic synchronization in bursting networking motifs.
Chaos ,18(3):037120, Sep 2008.[16] W.B. Kristan. Neuronal decision-making circuits.
Curr Biol ,18(19):R928–R932, Oct 2008. [17] K. L. Briggman and W. B. Kristan. Multifunctional pattern-generating circuits.
Annu Rev Neurosci , 31:271–294, 2008.[18] I.V. Belykh and A.L. Shilnikov. When weak inhibition synchro-nizes strongly desynchronizing networks of bursting neurons.
Phys Rev Lett , 101(7):078102, Aug 2008.[19] S. Jalil, I. Belykh, and AL. Shilnikov. Spikes matter inphase-locking of inhibitory bursting networks.
Phys. Rev. E. ,85:36214, 2012.[20] S Jalil, D. Allen, J. Youker, and AL. Shilnikov. Toward robustphase-locking in melibe swim central pattern generator models.
J. Chaos , 23(4), 046105–16, 2013.[21] D. Knapper, J Schwabedal and AL Shilnikov. Qualitative andquantitative stability analysis of penta-rhythmic circuits.
Non-linearity
Regular and Chaotic Dynamics
9, 281–297, 2004.[23] AL Shilnikov. Complete dynamical analysis of an interneuronmodel.
J. Nonlinear Dynamics , 68(3):305–328, 2012.[24] J. Collens, K. Pusuluri, A. Kelly, D. Knapper, T. Xing, S. Ba-sodi, D. Alacam and A.L. Shilnikov. Dynamics and bifurca-tions in multistable 3-cell neural networks
J. Chaos , J. Chaos30,072101, 2020 https://doi.org/10.1063/5.00113742020.[25] J. Wojcik, R. Clewley, and A.L. Shilnikov. Order parameter forbursting polyrhythms in multifunctional central pattern genera-tors.
Phys Rev E , 83:056209–6, 2011.[26] J. Wojcik, R. Clewley, J. Schwabedal, and A.L. Shilnikov. Keybifurcations of bursting polyrhythms in 3-cell central patterngenerators.
PLoS ONE , 9(4), 2014.[27] Shilnikov, A.L. and G. Cymbalyuk (2005). Transition betweentonic spiking and bursting in a neuron model via the blue-skycatastrophe.
Phys.Rev.Lett.
Scholarpe-dia
Inter. J. Bifurcation and Chaos
24, 1440003,2014, https://doi.org/10.1142/S0218127414400033.[30] A.L. Shilnikov, L.P. Shilnikov and D.V. Turaev, Moscow MathJ. (1), 205 (2005).[31] Channell, P., G. Cymbalyuk, and A. Shilnikov (2007a). Originof bursting through homoclinic spike adding in a neuron model. Phys.Rev.Lett.
Biol. Cybern. , 68:5, 1993.[34] B. Ermentrout and N. Kopell, Parabolic bursting in an excitablesystem coupled with a slow oscillation,