✬ ✫ ✩ ✪ Models of Multifunctional Central Pattern Generators: Polyrhythmic Bursting Jeremy Wojcik, Robert Clewley and Andrey Shilnikov Neuroscience Institute and Department of Mathematics and Statistics, Georgia State University, Atlanta, GA, USA ✎ ✍ ☞ ✌ Abstract We demonstrate a motif of three reciprocally inhibitory cells that is able to produce multiple patterns of bursting rhythms. Through the examination of the qualitative geometric structure of two-dimensional maps for phase lag between the cells we reveal the organizing centers of emergent polyrhythmic patterns and their bifurcations, as the asymmetry of the synaptic coupling is varied. The presence of multistability and the types of attractors in the network are shown to be determined by the duty cycle of bursting. This analysis does not require knowledge of the equations that model the system, and so provides a powerful new approach to studying regulatory networks. Thus, the approach is applicable to a variety of biological phenomena beyond motor control. ✎ ✍ ☞ ✌ 3 cell CPG motif A central pattern generator (CPG) is a neural microcircuit of cells whose synergetic interactions, without sen- sory input, rhythmically produce patterns of bursting that determine motor behaviors such as heartbeat, respi- ration, and locomotion in animals and humans [1]. While a dedicated CPG generates a single pattern robustly, a multifunctional CPG can flexibly produce distinct rhythms, such as swimming versus crawling, and alternation of blood circulation patterns, in leeches [2,3] Switching between such rhythms is attributed to input-dependent switching between attractors of the CPG. Here we analyze multistability and transformations of rhythmic patterns in a 9D model of a CPG motif (Fig. 1) comprised of 3 endogenously bursting leech heart interneurons (details in [4,8,9]. Our use of fast threshold modulation [5] implies that the post-synaptic current is zero (resp., maximized) when the voltage of a driving cell is below (resp., above) the synaptic threshold. This is an inherently bi-directionally asymmetric form of coupling. Many anatomically and physiologically diverse CPG circuits involve a 3-cell motif [6] , including the spiny lobster pyloric network, the Tritonia swim circuit, and the Lymnaea respiratory CPGs [7]. We show how rhythms of the multistable motif are selected by changing the relative timing of bursts by physiologically plausible perturbations. We also demonstrate how the set of possible rhythmic outcomes can be controlled by varying the duty cycle of bursts, and by varying the rotationally-asymmetric coupling around the ring. V 1 20 mV Θ th -50 mV V 2 20 mV -50 mV Δφ 21 (n+1) Δφ 21 (n) V 3 20 mV -50 mV Δφ 31 (n+1) Δφ 31 (n) (B) (A) 2s 2 3 1 Fig. 1. (A) 3-cell motif with asymmetric clockwise versus counter-clockwise connection strengths. (B) Voltage traces: the phase (mod 1) of reference cell 1 (blue) is reset when V 1 reaches threshold Θ th = −40 mV. The time between burst onsets in cell 2 (green) and 3 (red) determine phase lags {Δφ (n) 21 , Δφ (n) 31 }. AS the period of network oscillations varies over time, we define delays between cell 2 (green) and cell 3 (red) relative to the the reference cell 1 (blue) at the instances the voltages V i increase through a threshold of Θ th = −40 mV, thus indicating the onset of new bursting cycles (Fig. 1). Initial delays are controlled by the release of cells 2 and 3 from inhibition. The subsequent delays normalized over the period of the cell 1 define a sequence or forward trajectory { Δφ (n) 21 , Δφ (n) 31 } of phase-lag return maps on a torus [0, 1) × [0, 1) with Δφ i1 mod 1. The maps are tabulated on a 40 × 40 (or more) grid of initial points. The iterates are connected in figures below to preserve order. Now we can study the dynamical properties of the maps, locate fixed points and evaluate their stability, detect periodic and heteroclinic orbits, and identify the underlying bifurcations as the parameters, DC and g as , are varied. ✎ ✍ ☞ ✌ Leech heart interneuron model A reduced model of the leach heart interneuron is given by the following set of three nonlinear coupled differ- ential equations [8]: C dV dt = −I Na − I K2 − I L + I app + I syn , I L =¯ g L (V − E L ), I K2 =¯ g K2 m 2 K2 (V − E K ), I Na =¯ g Na m 3 Na h Na (V − E Na ), τ Na dh Na dt = h ∞ Na (V ) − h, m Na = m ∞ Na (V ), τ K2 dm K2 dt = m ∞ K2 (V ) − m K2 . (1) Here, C =0.5nF is the membrane capacitance; V is the membrane potential in mV; I Na is the sodium current with slow inactivation h Na and fast activation m Na ; I K2 is the slow persistent potassium current with activation m K2 ; I L is the leak current and I app =0.006 is an applied current. The values of maximal conductances are set as ¯ g K2 = 30nS, ¯ g Na = 200nS and g L =8nS. The reversal potentials are E Na =0.045 V, E K = −0.07 V and E L = −0.046 V. The time constants of gating variables are τ K2 =0.9 sec and τ Na =0.0405sec. The steady state values of gating variables, h ∞ Na (V ), m ∞ Na (V ), m ∞ K2 (V ), are given by the following Boltzmann equations: h ∞ Na (V ) = [1 + exp(0.5(V +0.325))] −1 m ∞ Na (V ) = [1 + exp(−0.15(V +0.0305))] −1 m ∞ K2 (V ) = [1 + exp (−0.083(V +0.018 + V shift K2 ))] −1 . (2) The synaptic current is modeled through the fast threshold modulation paradigm as follows [8,9]: I syn = − n ∑ j =1 g syn (1 ∓ g as )(E inh syn − V i )]Γ(V j − Θ syn ). (3) The reversal potential E inh syn = −0.0625 is set to be smaller than V i (t), i.e. the synapse is inhibitory. The synaptic coupling function is modeled by the sigmoidal function Γ(V j )=1/[1 + exp{−1000(V j − Θ syn )}]. The threshold Θ syn = −0.03 V is chosen so that every spike within a burst of the neuron burst can reach it. This implies that the synaptic current from the j -th neuron is initiated as soon as this neuron becomes active after its membrane potential exceeds the synaptic threshold. The intrinsic bifurcation parameter V shift K2 of the model is a deviation from V 1/2 =0.018 mV corresponding to the half-activated potassium channel at m ∞ K2 =1/2. In the model (1), decreasing V shift K2 elevates the slow nullcline dm K2 dt =0 in the V -direction, thereby delaying the activation of m K2 . In this study, we V shift K2 = −0.01895 mV, V shift K2 = −0.021 mV and V shift K2 = −0.0225 mV, for short (∼20%), medium (∼50%), and long (∼80%) bursting, respectively. ✓ ✒ ✏ ✑ Poincar ´ e return maps for phase lags Fig. 2A shows the (Δφ 31 , Δφ 21 ) phase-lag map for the S 3 -symmetric, medium bursting motif at V shift K2 = −21.0 mV. The map has five stable FPs (color-coded dots) corresponding to the coexisting phase-locked rhythms: red at (Δφ 31 ≈ 1 2 , Δφ 21 ≈ 0), green (0, 1 2 ), blue ( 1 2 , 1 2 ), black ( 2 3 , 1 3 ) and purple ( 1 3 , 2 3 ). The attraction basins of the FPs are divided by separatrices (incoming and outgoing sets) of six saddles (grey dots). (A) Fig. 2. (A) 3-cell motif with asymmetric clockwise versus counter-clockwise connection strengths. (B) Voltage traces: the phase (mod 1) of reference cell 1 (blue) is reset when V 1 reaches threshold Θ th = −40 mV. The time between burst onsets in cell 2 (green) and 3 (red) determine phase lags {Δφ (n) 21 , Δφ (n) 31 }. (B) (Δφ 21 , Δφ 31 ) phase-lag return maps on a torus [0, 1) × [0, 1). The coordinates of the FPs determine the phase locked states within bursting rhythms of the motif that are denoted symbolically as follows: (1 ≺ 2 ≺ 3) and (1 ≺ 3 ≺ 2) for clockwise and counter-clockwise traveling waves of bursting (resp.) around the ring (episodes (ii) and (iii) of the voltage trace in Fig. 2) corresponding to the FPs located near ( 1 3 , 2 3 ) and ( 2 3 , 1 3 ), respectively. Besides these rhythms (which stability is unknown a priori) in a symmetric motif, the three other FPs correspond to the rhythms where one cell fires in anti-phase (Δφ ∼ 1 2 , or ⊥) against two cells bursting simultaneously in-phase (Δφ =0, or ∥). The stable FP at ( 1 2 , 1 2 ) corresponds to the (1 ⊥{2 ∥ 3})-pattern (episode (i) of Fig. 2); FP (0, 1 2 ) corresponds to (2 ⊥{1 ∥ 3})-pattern (episode (iv)); and FP ( 1 2 , 0) corresponds to (3 ⊥{1 ∥ 2})-pattern (episode (v)). V 1 V 2 7 s V 3 (i) (ii) (iii) (iv) (v) Fig.3. Five polyrhythms in the medium bursting motif at g syn =5 × 10 -3 . Inhibitory pulses (horizontal bars) suppress the targeted cells, thus causing switching between co-existing rhythms: (1 ⊥{2 ∥ 3}) in episode (i), traveling waves (1 ≺ 2 ≺ 3) in (ii) and (1 ≺ 3 ≺ 2) in (iii), followed by (2 ⊥{1 ∥ 3}) led by cell 2 in (iv). Having released cells 1 and 2 simultaneously, cell 3 leads the motif in the (3 ⊥{1 ∥ 2}) rhythm in the last episode. Colors match those of the stable FPs in Figs.3–7. ✓ ✒ ✏ ✑ Polyrhythmic bursting toward inhibitory ring Next we examine bifurcations of FPs in asymmetric ring by increasing g as , which weakens counterclockwise and enforces clockwise-directed synapses. As g as increases to 0.154, the three saddles move away from the FP at ( 2 3 , 1 3 ), thereby increasing its attraction basin, and approach the stable FPs corresponding to anti-phase bursting rhythms. Meanwhile, the other three saddles move towards the FP at ( 1 3 , 2 3 ) corresponding to the (1 ≺ 3 ≺ 2) rhythm, narrowing its basin. As g as is increased further, the attractors and saddles in the bottom right of the unit square merge and vanish through saddle-node (SN) bifurcations. Increasing g as makes the FP ( 2 3 , 1 3 ) for the (1 ≺ 2 ≺ 3) rhythm globally dominant (Fig. 5B) after the three nearby saddles collapse onto the FP at ( 1 3 , 2 3 ). (A) (B) Fig. 4. (A) Asymmetric motif at g as =0.154 near the saddle-node bifurcations annihilating three stable FPs for anti-phase bursting rhythms. (B) Map for the asymmetric (g as =0.3) medium bursting motif depicting two persistent attractors: (1 ≺ 2 ≺ 3) prevails over (1 ≺ 3 ≺ 2). Further increase in g as leads to the only observable (1 ≺ 2 ≺ 3) rhythm, after FP ( 1 3 , 2 3 ) merges with nearby saddles. ✓ ✒ ✏ ✑ Duty cycle as Order parameter Comparison of the maps for the symmetric motifs demonstrates that Duty Cycle is an order parameter that determines the observable rhythms. I.e., short bursting makes has no clockwise (1 ≺ 2 ≺ 3) and counter-clockwise (1 ≺ 3 ≺ 2) traveling wave patterns. In contrast, the symmetric long bursting motif is unlikely to generate anti-phase rhythms, as the corresponding FPs have narrow attraction basins, divided equally between the soaringly dominating FPs corresponding to the (1 ≺ 2 ≺ 3) and (1 ≺ 3 ≺ 2) rhythms. SN bifurcations are same in the all motifs. The short motif undergoes another bifurcation that makes the (1 ≺ 2 ≺ 3) pattern the global attractor: this FP at ( 2 3 , 1 3 ) becomes stable through a supercritical Andronov-Hopf bifurcation. (A) (B) Fig. 5. (A) Phase-lag map for the symmetric, short bursting motif at V shift K2 = −18.95 mV, showing only 3 stable FPs (blue, red, and green) corresponding to anti-phase rhythms where one cell bursts followed by synchronized bursts in the other two cells. Unstable FPs at ( 2 3 , 1 3 ) and ( 1 3 , 2 3 ) exclude the clockwise (1 ≺ 2 ≺ 3) and counter-clockwise (1 ≺ 3 ≺ 2) traveling waves from the repertoire of the short bursting motif. (B) Map for the asymmetric motif (g as =0.2) depicting a stable invariant curve near a three-saddle connection around the FP at ( 2 3 , 1 3 ) The invariant curve is associated with the appearance of slow phase jitters within the (1 ≺ 2 ≺ 3) rhythm in voltage traces. Once it collapses onto the FP at ( 2 3 , 1 3 ) the asymmetric motif gains a new robust (1 ≺ 2 ≺ 3) rhythm, making four possible bursting outcomes. (A) Fig. 6. (A) Phase-lag map for the symmetric (A) long bursting motif at V shift K2 = −22.5 mV, revealing two equally dominant attractors: ( 2 3 , 1 3 ) and ( 1 3 , 2 3 ) for (1 ≺ 2 ≺ 3) and (1 ≺ 3 ≺ 2) traveling rhythms. (B) Asymmetric motif showing the dominance of the (1 ≺ 2 ≺ 3) rhythm over the (1 ≺ 3 ≺ 2)-rhythm whose basin of attraction is bounded by saddles – unstable phase locked states. ✓ ✒ ✏ ✑ Breaking symmetry: making dedicated pacemakers Fig.7. (A) Medium bursting motif without the (1 ≺ 2 ≺ 3) rhythm at one 5% stronger g 21 -connection. (B) Motif with g 21 = g 21 connections 50% stronger than all other connections, displaying a phase-slipping, which is due to cycling by-passing (1 ⊥{2 ∥ 3}) → (3 ⊥{1 ∥ 2}) through the ghosts of the vanished rhythms. The remaining two saddles define a narrow basin of attraction for the only (2 ⊥{1 ∥ 3}) phase locked rhythm. ✓ ✒ ✏ ✑ Acknowledgment We acknowledge support from NSF Grants CISE/CCF-0829742 (to R.C.), DMS-1009591 and RFFI Grant No. 08-01-00083 (to A.S.), GSU Brains & Behavior program, and Brain Corporation. ✎ ✍ ☞ ✌ References 1. A. Selverston (ed), Model Neural Networks and Behavior, (Springer, 1985); F. Skinner, N. Kopell, E. Marder, Comput. Neurosci. 1, 69 (1994); E. Marder, R. L. Calabrese, Physiol. Rev. 76, 687 (1996); N. Kopell, G. B. Ermentrout, PNAS USA 101, 15482 (2004). 2. K. L. Briggman, W. B. Kristan, Annu. Rev. Neurosci. 31 (2008). 3. W. B. Kristan, R. L. Calabrese, W. O. Friesen, Prog. Neurobiol. 76, 279 (2005). 4. A. Shilnikov, R. Gordon, I. Belykh, Chaos 18, 037120 (2008). 5. N. Kopell, D. Somers, Biol. Cybern. 68, 5 (1993). 6. R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, U. Alon, Science 298 (2002); M. I. Rabinovich, P. Varona, A. Selverston, and H. D. Abarbanel, Rev. Modern Phys. 78, 4 (2006). 7. W.N. Frost and PS Katz, PNAS 93, 422 (1996); A. G. M. Bulloch and N. I. Sayed, Trends Neurosci. 15, 11 (1992). 8. A. Shilnikov and G. Cymbalyuk, Phys. Rev. Lett. 94, 048101 (2005). 9. I. Belykh and A. Shilnikov, Phys. Rev. Lett. 101, 078102 (2008). Nature Precedings : doi:10.1038/npre.2011.5800.1 : Posted 17 Mar 2011