Generalized half-center oscillators with short-term synaptic plasticity
Valentina Baruzzi, Matteo Lodi, Marco Storace, Andrey Shilnikov
AAPS/123-QED
Generalized half-center oscillators with short-term synapticplasticity
V. Baruzzi, M. Lodi, M. Storace, and A. Shilnikov Department of Electrical, Electronics andTelecommunication Engineering and Naval Architecture,University of Genoa, 16145 Genoa, Italy Department of Mathematics and Statistics,Neuroscience Institute, Georgia State University, Atlanta, GA 30303 USA. (Dated: May 22, 2020)
Abstract
How can we develop simple yet realistic models of the small neural circuits known as centralpattern generators (CPGs), which contribute to generate complex multi-phase locomotion in livinganimals? In this paper we introduce a new model (with design criteria) of a generalized half-centeroscillator (gHCO), (pools of) neurons reciprocally coupled by fast/slow inhibitory and excitatorysynapses, to produce either alternating bursting or synchronous patterns depending on the sensoryor other external input. We also show how to calibrate its parameters, based on both physio-logical and functional criteria and on bifurcation analysis. This model accounts for short-termneuromodulation in a bio-physically plausible way and is a building block to develop more realisticand functionally accurate CPG models. Examples and counterexamples are used to point out thegenerality and effectiveness of our design approach. a r X i v : . [ q - b i o . N C ] M a y . INTRODUCTION Central pattern generators (CPGs) are small neural circuits that can autonomously (i.e.,in the absence of sensory feedback or higher motor planning centers inputs) produce var-ious rhythmic patterns of neural activity [17]. They bear a fundamental function in bothinvertebrate and vertebrate animals as they determine multi-phase locomotion – the innatemotor behavior that requires sequential activation of body muscles in a coordinated way[22]. Various approaches to the modeling of CPGs and CPG-inspired control systems havebeen explored in the last decades [7, 12, 18, 30, 37]. Recently, new methods have been pro-posed to reduce large models of detailed neural networks to smaller CPG circuits, tradingoff biological plausibility and complexity of the model [2, 7, 25, 27, 30].Although CPGs function autonomously, their activity is modulated through the influenceof hierarchically higher areas, which can, for example, prompt transitions between gaits[9, 15, 34]. A single gait in a typical CPG model is obtained by fixing the connectivity.By contrast, to generate multiple gaits the CPG connections between constituent neuronsare typically changed acting on the synaptic weights to model the control action of thebrainstem [12, 25–27]. The modulation from higher areas that controls the synchronizationbetween the CPG neurons, and thus triggers gait switches, is conveniently integrated inCPG models to directly affect the synaptic conductance strengths. However, in real CPGschanges in conductance values are the result of long-term synaptic plasticity, and thereforeit is hardly a cause for quick gait switches, which can instead be accounted for more re-alistically by short-term neuromodulation. Indeed, most natural CPGs exhibit patterns offunctional connectivity between neurons or synchronized clusters of neurons that can un-dergo spontaneous fluctuations and be highly responsive to perturbations, e.g., induced bysensory input or cognitive tasks, on a timescale of milliseconds or hundreds of milliseconds,respectively, thus ensuring robustness and stability. This short-term neuromodulation lacksin most CPG models.One of the pivotal building blocks of many CPGs is a half-center oscillator (HCO). TheHCO-concept is widely used to model two synchronous pools of neurons reciprocally in-hibiting each other to produce stable rhythmic alternation in animal locomotion [6, 10].This basic structure has been largely studied from both biological and nonlinear dynamicsstandpoints. For example, in [4] transitions between stable synchronous states in the HCO2ccur through direct manipulations with synaptic weights, whereas in [14] a large databaseof HCO models is swept using a brute-force approach, without a focus on gait transitions.While the importance of an interplay between inhibitory and excitatory coupling has alreadybeen outlined [4], the thorough understanding of its functional role for determining multiplestates or patterns in such neural networks and how transitions between them may stablyoccur remains yet insufficient. Moreover, there is the growing evidence that (i) post-synapticpotential (PSP) summation increasing with the spike frequency in the pre-synaptic cell is acrucial factor for stable functioning of some CPGs [11, 16, 29, 32], while other experimentsindicate that (ii) the activity of some synapses is barely affected by the spike frequency [12].In this paper, we propose a generalized half-center oscillator (gHCO) composed of two neu-rons or of two neural pools that are coupled reciprocally by excitatory synapses, in additionto the standard HCO’s reciprocally inhibitory synapses. We show that this circuitry war-rants a more biologically plausible mechanism of short-term plasticity to implicitly controlthe phase-lag between the gHCO cells by varying their spike frequency through sensory driveor external currents, rather then directly manipulating the synaptic conductance strengths.Moreover, we show how to calibrate the gHCO parameters in order to obtain the desiredbehaviors, also carrying out a numerical bifurcation analysis.
II. THE GHCO AND ITS DESIGN CONSTRAINTS
The proposed generalized half-center oscillator is shown in Fig. 1. It is made of twoneurons or two neural pools, coupled by both excitatory (marked by a black circle) andinhibitory (marked by a black triangle) synapses.There are a few simple constraints that neurons and synapses must meet for the circuit togenerate stably the desired rhythmic outcomes: (a) both neurons are endogenous bursterswith (b) the spiking voltage range above the hyperpolarized voltage (i.e., they do not un-dershoot [19]) within each burst, while (c) the mean spike frequency can be controlled. ThegHCO bursters are coupled by (d) slow synapses with PSP summation whose strength in-creases with the growing spike frequency in presynaptic cells, as well as by (e) fast synapseswithout PSP summation.In what follows, both gHCO cells are represented by the Hodgkin-Huxley (HH) type modelof the thalamic reticular neuron [13, 28] (see Appendix). This slow-fast model with seven3tate variables can exhibit endogenous bursting activity of alternating trains of fast actionpotentials with long quiescent intervals, as depicted in Fig. 2. The dynamics of the mem-brane potential V j and of the voltage-dependent state variables (the vector y j ) are governedby a generic set of HH-like equations ddt V j y j = − (cid:80) k I k + I synj f ( V j , y j ) , where j = 1 , . (1)where f ( V j , y j ) is a vector function describing y j -dynamics; in particular, each f componentfor the HH gating variables is a logistic function. In addition to intracellular currents, (cid:80) k I k includes a further external contribution, namely a control current I c acting essentially onthe spike frequency within bursts. For the given model, bursting activity occurs when I c ∈ [ − . , . (cid:2) µAcm (cid:3) , with the mean spike frequency decreasing from 15.36 to 4.13 ms.The term I synj is the incoming mixed, excitatory/inhibitory synaptic current originatingfrom the i -th cell onto the j -th, post-synaptic cell: I synj = g ex ( E ex − V j ) s exi + g in ( E in − V j ) s ini , (2)where E ex/in are the reversal potentials for excitatory/inhibitory synapses and 0 ≤ s ex/ini ≤ V j < E ex ) orinhibitory ( V j > E in ). For the slow synapses with PSP summation we employ a first-orderdynamic synapse [8, 21, 35]. The dynamic evolution of its activation rate is governed by thefollowing equation ds i dt = α (1 − s i ) f ∞ ( V i ) − βs i , f ∞ = 11 + e − ν ( V i − θ ) , (3)where θ is the synaptic threshold, whereas α and β are coefficients weighting the raise anddecay terms, respectively. To model the static synapses without PSP summation we employ cell cell FIG. 1. (color online). gHCO neural circuit with inhibitory (denoted with (cid:73) (cid:74) ) and excitatory( • ) synapses reciprocally coupling two oscillatory cells. ≤ s i = f ∞ ( V i ) ≤
1, with θ being below the spike-level.To illustrate the contrasting properties of these synapse models, we refer to Fig. 2, showingthe bursting voltage traces V (red) and V (blue) and the synaptic activation dynamics, fast s in ( t ) (gray) and slow s ex ( t ) (black) at the edge of the I c bursting interval. Observe that theneurotransmitter release rate s in ( t ) of the fast FTM synapse (1) is maximized as soon asthe voltage V ( t ) in the pre-synaptic cell overcomes the synaptic threshold θ in (indicated bythe grey lines in panels a,b), (2) remains constant regardless of the spike frequency, and (3)vanishes with the burst termination. In contrast, the low spike frequency (panels a,c) barelyactivates the slow synapse (see s ex ( t )) that at high spike frequency (panels b,d) exhibits theprofound PSP build up; the ascending rate is ruled by α >
0, and the exponential decaydue to β > θ . III. PARAMETER CALIBRATION
The neuron and synapse models (1-3) are calibrated to physiologically plausible valuesto meet the above requirements (a)-(e) and to ensure a smooth and reversible transitionfrom anti-phase to in-phase bursting occurring in the gHCO as the spike frequency changesdue to I c -variations. Just to clarify things, let us consider the dynamics of the gHCO withfast FTM inhibitory and slow excitatory synapses. Moreover, the corresponding synapticthresholds are set at θ in = −
30 and θ ex = 25mV, respectively. As such, the inhibitorysynapses without PSP summation (de)-activate quickly and their strength remain constantduring each burst regardless of the spike frequency. In contrast, the slow excitatory synapsesexhibit PSP summation that becomes stronger with an increase of the spike frequency.Figure 2 shows that at the low end I c = − .
43 of the bursting region, near the transition tothe hyperpolarized quiescence, the gHCO neurons oscillate in anti-phase with the smallestnumber of spikes per burst and lowest spike frequency (panels (a, c)), whereas on the oppositeside at I c = 0 .
13 the neurons burst in phase with a larger number of spikes per burstand with much higher spike frequency (panels (b, d)). Changing the value of I c changesthe strength of the excitatory synapses, and hence the proportion between inhibition andexcitation that repel the gHCO neurons or attract them to each other, respectively. The phase-lag ∆ (defined on mod 1) between burst initiations in the neurons [20, 36, 38] allows5uantifying the phase-locked states produced by the gHCO. In case of the synchronous orin-phase bursters, ∆ = 0 (or ∆ = 1). When they burst in alternation, with ∆ = 0 .
5, wesay that they are in anti-phase. The intermediate values of ∆ correspond to “winner-less”patterns transitional between the in- and anti-phase states generated by the gHCO.The bifurcation analysis of the system (1–3) was carried out using the computationaltoolbox CEPAGE [24]. Since we want the gHCO to transition from anti-phase regime toin-phase regime varying I c , we need the proportion between inhibition and excitation to besignificantly different for the two values of I c at the edges of its range. To this end, we seekmaximum difference in the mean values of s exi (over one period) at the two extreme valuesof I c , i.e., -0.43 (anti-phase pattern) and 0.13 (in-phase bursting). We set the numericalvalues of θ ex , α and β according to this principle, running a set of simulations over a
0 0.51 -60-300 30 60
FIG. 2. (color online) Asymptotic anti-phase (a) and synchronous (b) bursting voltage traces V (red) and V (blue) at I c = − .
43 and 0.13, resp., in gHCO (1)-(3), superimposed with exci-tatory/inhibitory thresholds θ (horizontal lines) at 10 and -30 mV. (c, d) Synapse dynamics: fastmodulatory s in ( t ) (gray) vs. slowly summating/decaying s ex ( t ) (black). See the Appendix forparameters. θ ex = { , } , 10 evenly spaced values of α ∈ [0 . ,
1] and10 evenly spaced values of β ∈ [0 . , . θ ex indicate voltagelevels representative of two different conditions: at θ ex = 10 each spike appears broader,i.e. V j stays above θ ex for a longer time window; at θ ex = 25 each spike appears narrower,i.e. V j stays above θ ex for a shorter time period. We choose the parameter setting thatprovides maximum difference in the mean values of s exi for the two extreme values of I c (seeAppendix B). The synaptic conductances g in/ex are set to obtain anti-phase synchronizationfor low spike frequency, condition in which the mean value of s exi is minimum, and in-phase synchronization for high spike frequency, condition in which the mean value of s exi ismaximum.The results are summarized in Fig. 3, and reveal the dependence of the phase-lag ∆on the I c -current, and hence explicitly on the spike frequency within bursts. As expected,at low I c -values between − .
43 and − .
40, the fast reciprocal inhibition within the gHCOdominates and makes its neurons burst in alternation with ∆ = 0 .
5. As the I c -currentis increased, the spike frequency raises, which in turn makes the slow excitatory synapsessum up faster and stronger on average. With larger I c values, the reciprocal excitationgradually prevails over the reciprocal inhibition, which gives rise to the smooth onset of the -0.4 -0.3 -0.2 -0.1 0 0.1 0 0.51 FIG. 3. Bifurcation diagram showing how the phase-lag ∆ between the gHCO neurons is affectedby the current I c ; here, 30 initial ∆-values were sampled evenly between 0.05 and 0.95 for each ofthe 50 I c -values. Parameters listed in the Appendix. I c -current. Wenote also that this diagram has been obtained by making a multi-shooting for each parametervalue. This is a direct indication that there is no hysteresis and therefore the absence ofmulti-stability or the coexistence of anti- and in-phase bursting for same parameter values,and that the transition between activity rhythms is continuous and reversible. We wouldlike to re-emphasize that the maximal synaptic conductances g in/ex in Eq. (2) once set arenot changed, and the transition is solely determined by the gradual increase/decrease of themean s exi -value caused by the spike frequency variations in the gHCO neurons. -60-40-20020-60-40-20020 FIG. 4. (color online). (a) Asymptotic bursting voltage trace with undershoot produced by thePlant neuron model [1, 31]. (b) Voltage traces produced by the gHCO with two coupled Plantneurons. See the Appendix for parameters. V. COUNTEREXAMPLES
The proposed gHCO concept can fall apart whenever one or more of the conditions onthe neuron and synapse models are not fulfilled. If the bursting condition (a) is broken,the approach is no longer applicable. Two neurons, spiking in isolation, can burst in al-ternation due to reciprocal inhibition, but not through reciprocal excitation, which makesboth even more synchronously depolarized with a higher frequency. If the neurons under-shoot (condition (b)), which is typical for elliptical bursters [1] (see Fig. 4(a)), the choiceof the inhibitory threshold θ in to warrant evenly constant activation s ini requires additionalconsiderations. Indeed, this choice can result in less robust dynamics of the gHCO, due toinhibition-excitation competition (see Fig. 4(b)). Condition (c), outlining the importanceof being able to control spike-frequency and not only burst duration of the pre-synapticcell, is quite crucial for stable gHCO functions. To point out its significance, we employthe exponential integrate-and-fire (eIF) neuron model [5], where an external current I ext primarily controls the burst duration with insignificant spike-frequency variations, as shownin Fig. 5(a). In this scenario, the activation of both inhibitory and excitatory synapses ismainly determined by the burst duration in the eIF-neurons, and thus I ext -variations canonly cause proportional changes in the average excitatory s exi - and inhibitory s ini -values. Asa result, neither inhibition nor excitation can solely dominate and produce the expectedsolo stable anti-phase or in-phase bursting patterns within the given I ext -range, as shown inFig.5(b). Conversely, changing the parameter g e of the eIF neuron model significantly mod-ifies the spike frequency, and the corresponding bifurcation diagram has the characteristicpitchfork shape, as expected. However, the parameter g e is a conductance, and thus is nota realistic control parameter, according to our guidelines. Condition (d) follows (c), as thesynaptic threshold θ , for the slow synapses, has to be within the spike voltage range of thepre-synaptic neuron and the dynamics is to be slow enough to allow s i ( t ) to grow and thesynapse to exhibit PSP summation. Condition (e) guarantees that the activation of the fastsynapse does not exhibit PSP summation and hence does not change due to spike frequencyvariations in the pre-synaptic neuron. 9
00 750 800 850 900 950 1000 1050 110050607080 3.54.55.5700 750 800 850 900 950 1000 1050 110000.51
FIG. 5. (color online). (a) Mean values (over 5 s) of the IBI (green line) and the burst dura-tion (black line) plotted against I ext for the exponential IF-model [5]. Corresponding bifurcationdiagram for the phase lag ∆ between the cells in the gHCO, in which each cell is an exponentialIF-model (b). V. TOWARDS A LOCOMOTION CPG
As the gHCO often happens to be a CPG building block, we discuss some solutionsensuring that both the phase lags and the burst frequency are consistent for the modeledgaits. For instance, in left-right alternation of the mouse locomotion, a phase lag ∆ = 0 . I c -values and slow bursting at greater I c -values. Therefore, for the gHCO built withsuch models to produce in-phase/anti-phase synchronization at high/low burst frequencies10 FIG. 6. (a) Mean values (over 5 s) of the IBI plotted against g e for the exponential IF-model [5].Corresponding bifurcation diagram for the phase lag ∆ between the cells in the gHCO, in whicheach cell is an exponential IF-model (b). for the desired gaits, the time-scale of the synapses in its circuitry should be swapped:slow inhibitory synapses with PSP summation and FTM-fast excitatory ones without PSPsummation, see the Appendix for details. Moreover, we use a modified version of the first-order synapse to model slow inhibitory synapses. The dynamics of its activation is governedby the following equation ds i dt = α s i (1 − s i ) f ∞ ( V i ) − βs i (4)where the new multiplicative term delays and hence slows down the synaptic activationfor low spike-frequency in the pre-synaptic neuron; the synapse remains inactive near I c = − .
43. The synapse given by Eq. 4 maintains a greater contrast in the mean s i -valuescorresponding to the low and high ends of the bursting I c -range for the given neuron model.11he results are summarized in Fig. 7, representing the bifurcation diagram for this gHCO. Itdemonstrates that the gHCO bursters oscillate robustly in-phase (∆ = { , } ) for negative I c -values and rapidly transition to the stable anti-phase (∆ = 0 . I c is step-wise increased from -0.43 to 0.13. VI. CONCLUDING REMARKS
We developed a generalized HCO-model with a short-term plasticity mechanism, whichaccounts for short timescale gait transitions induced by sensory input or cognitive tasks.The proposed concept is based on simple constraints (i) subjecting models for cells andsynapses and (ii) optimizing the trading-off between physiological plausibility and modelfunctionality. The generality of our approach suggests that it will be applicable for otherbiologically plausible and phenomenological models of endogenous (square-wave) bursters,and for other dynamic synapse models. -0.4 -0.3 -0.2 -0.1 0 0.1 0 0.51
FIG. 7. Bifurcation diagram showing the flat-even phase-lags, ∆ = { , } (in-phase) and ∆ = 0 . I c for the gHCO with slow inhibitoryand fast excitatory synapses; here, 30 initial ∆-values were sampled evenly between 0.05 and 0.95for every I c value out of 50. Parameters listed in the Appendix. CKNOWLEDGMENTS
We would like to acknowledge J. Scully’s contribution to the concept and developmentof the synapse model Eq. (4). A.S.’s research was partially funded by the NSF grant IOS-1455527. M.S. and A.S. conceptualized the work; V.B. and M.L. conducted the experiments. -5005000.5 -0.13-0.1
FIG. 8. (color online) (a) Time evolution of the phase lag ∆ between the gHCO cells (black line)in response to step-wise changes of I c (green dashed lines); I c increased over 25 steps from − .
43 to0 .
13, only the time window in which ∆ transition occurs is shown. (b) Voltage traces progressingfrom in-phase to anti-phase bursting within the time window bounded by the grey vertical lines in(a). ppendix A: Neuron Models1. Thalamic Reticular Neuron Model The thalamic reticular neuron model [13, 28] is defined by the following state equations: dVdt = − I T − I L − I Na − I K − I c + I syn CdCadt = − kI T F d − K T CaCa + K d dydt = y ∞ − yτ y y = { h, m, n, m T , h T } (A1)where V is the membrane potential of the neuron; the ion currents I T (calcium), I Na (sodium), I K (potassium), and I L (leakage) evolve according to the following equations I T = g Ca m T h T ( V − E Ca ) , I L = g L ( V − E L ) ,I Na = g Na m h ( V − E Na ) , I K = g k n ( V − E k ) , which depend on V , on the intracellular calcium concentration Ca and on a set of furtherstate variables (called gating variables ) h , m , n , m T , h T . The differential equations governingthese gating variables have the common structure written above (for the generic gatingvariable y ), where: y ∞ = a y / ( a y + b y ) , τ y = 1 / ( a y + b y ) ( y = { h, m, n } ) a h = 0 . e − V , b h = 4 e − . V − + 1 ,a m = 0 . − V ) e . − V ) − , b m = 0 . V − e . V − − a n = 0 . − V ) e . − V ) − , b n = 0 . e − V m ∞ T = 11 + e − V +527 . , τ mT = 0 .
44 + 0 . e V +2710 + e − V +10215 ,h ∞ T = 11 + e V +805 , τ hT = 62 . . e V +484 + e − V +40750 . In the above equations, h and m are the inactivation and activation variables of the N a + current; n is the activation variable of the K + current; m T and h T are the activationand inactivation variables of the low-threshold Ca current; the leakage current I L hasconductance g L = 0 .
05 [ mScm ] and reversal potential E L = −
78 [ mV ]; I Na and I K are the14ast N a + and K + currents responsible for the generation of action potentials, with con-ductances g Na = 100 [ mScm ] and g k = 10 [ mScm ] and reversal potentials E Na = 50 [ mV ] and E k = −
95 [ mV ]; I T is the low-threshold Ca current that mediates the rebound burst re-sponse, with conductance g Ca = 1 .
75 [ mScm ] and reversal potential E Ca = k RT F log( Ca Ca ); I syn is the synaptic current (Eq. (2) in the paper).When the control current I c is in the range [ − . , .
13] [ µAcm ] the neuron exhibits burstingbehavior. The other parameters are set as follows: C = 1 [ µFcm ] , Ca = 2 [ mM ] , d =1 [ µm ], K T = 0 . mM · ms ] , K d = 0 . mM ]. F = 96 .
489 [
Cmol ] is the Faradayconstant, R = 8 . Jmol · K ] is the universal gas constant and the temperature T is set at309 .
15 [ K ].
2. Exponential Integrate and Fire Neuron Model
The exponential integrate and fire (eIF) neuron model [5] is defined by the following stateequations: dVdt = − g L ( V − E L ) + g e e V − VT ∆ T − u + I ext + I syn Cdudt = a ( V − E L ) − uτ w (A2)where V is the membrane potential of the neuron; u is the adaptation variable; g L = 30 [ nS ]is the leakage conductance and E L = − . mV ] is the leakage reversal potential; I syn isthe synaptic current (Eq. (2) in the paper).When the conductance g e is set at 110 [ nS ], the external current I ext is varied in the range[690 , pA ] (Fig. 5). When the external current I ext is set at 800 [ pA ], the conductance g e is varied in the range [20 , nS ] (Fig. 6). For this range of parameter values, the neuronexhibits bursting behavior. The other parameters are set as follows: C = 2007 . pF ] , V T = − . mV ] , ∆ T = 2 [ mV ] , τ w = 285 . ms ] , a = 4 [ nS ].
3. Plant Neuron Model
The Plant neuron model [1, 31] is defined by the following state equations:15 dVdt = − I T − I L − I Na − I K − I KCa + I ext + I syn CdCadt = ρ ( K c x ( V Ca − V ) − Ca ) dydt = y ∞ − yτ y y = { h, n, x } (A3)where I T = g T x ( V − E I ) , I L = g L ( V − E L ) ,I Na = g I m ∞ h ( V − E I ) , I K = g K n ( V − E K ) ,I KCa = g KCa
CaCa + 0 . V − E K ) ,m ∞ = . − V s ) e − Vs . − V s ) e − Vs + 4 e − Vs ,h ∞ = 0 . e − Vs . e − Vs + e − Vs ,τ h = 12 . . e − Vs + e − Vs ,n ∞ = . − V s ) e − Vs − . − V s ) e − Vs − . e − Vs ,τ n = 12 . . − V s ) e − Vs − . e − Vs ,x ∞ = 1 e . − V − + 1 ,V s = 127 V
105 + 8265105where V is the membrane potential of the neuron; Ca is the intracellular calcium concen-tration; x is the activation variable of the slow inward Ca current; h is the inactivationvariable of the N a + current; n is the activation variable of the K + current; I L is the leak-age current, with conductance g L = 0 .
003 [ nS ] and reversal potential E L = −
40 [ mV ]; I Na and I K are the fast inward N a + and outward K + currents, respectively, with conductances g I = 8 [ nS ] and g K = 1 . nS ] (these values ensure undershoot, see paper) and reversalpotentials E I = 30 mV and E K = −
75 [ mV ]; I T is the slow inward tetrodotoxin-resistant16 a current, with conductance g T = 0 .
01 [ nS ] and reversal potential E T = 30 [ mV ]; I KCa is the outward Ca sensitive K + current, with conductance g KCa = 0 .
03 [ nS ] and reversalpotential E K ; I syn is the synaptic current (Eq. (2) in the paper).The external current I ext is set to 0 .
028 [ µA ]. The other parameters are set as follows: C =1 [ µFcm ] , ρ = 0 . mV − ] , K c = 0 . mV − ] , V Ca = 140 [ mV ] , τ x = 235 [ ms ]. Appendix B: Synapse Parameter Values
In Table I, column A lists the parameter values used for the gHCO with the thala-mic reticular neuron model, first-order dynamic excitatory synapses and static inhibitorysynapses (Figs. 2 and 3). Column B lists the parameter values used when simulating thegHCO with the thalamic reticular neuron model, modified first-order dynamic inhibitorysynapses (Eq. 4) and static excitatory synapses (Figs. 7 and 8 in the paper). Column Clists the parameter values used for the gHCO with the eIF neuron model when varying I ext ,first-order dynamic excitatory synapses and static inhibitory synapses (Fig. 5). Column Elists the parameter values used for the gHCO with the eIF neuron model when varying g e ,first-order dynamic excitatory synapses and static inhibitory synapses (Fig. 6). ColumnE lists the parameter values used for the gHCO with the Plant neuron model, first-orderdynamic excitatory synapses and static inhibitory synapses (Fig. 4). [1] Deniz Ala¸cam and Andrey Shilnikov. Making a swim central pattern generator out of latentparabolic bursters. Int. J. Bifurcat. Chaos , 25(07):1540003, 2015.[2] Jessica Ausborn, Abigail C Snyder, Natalia A Shevtsova, Ilya A Rybak, and Jonathan ERubin. State-dependent rhythmogenesis and frequency control in a half-center locomotor cpg.
J. Neurophysiol. , 119(1):96–117, 2018.[3] Carmelo Bellardita and Ole Kiehn. Phenotypic characterization of speed-associated gaitchanges in mice reveals modular organization of locomotor networks.
Curr. Biol. , 25(11):1426–1436, 2015.[4] Tiaza Bem and John Rinzel. Short duty cycle destabilizes a half-center oscillator, but gapjunctions can restabilize the anti-phase pattern.
J. Neurophysiol. , 91(2):693–703, 2004. ABLE I. Parameter values.A B C D E α ex β ex θ ex
25 -30 -40 -40 -42 g ex E ex
60 60 20 20 50 α in - 0.5 - - - β in - 0.02 - - - θ in -30 25 -48.5 -48.5 -53 g in E in -80 -80 -110 -110 -80 ν
10 10 10 10 10[5] Romain Brette and Wulfram Gerstner. Adaptive exponential integrate-and-fire model as aneffective description of neuronal activity.
J. Neurophysiol. , 94(5):3637–3642, 2005.[6] T Graham Brown. On the nature of the fundamental activity of the nervous centres; togetherwith an analysis of the conditioning of rhythmic activity in progression, and a theory of theevolution of function in the nervous system.
J. Physiol. , 48(1):18–46, 1914.[7] Pietro-Luciano Buono and Martin Golubitsky. Models of central pattern generators forquadruped locomotion i. primary gaits.
J. Math. Biol. , 42(4):291–326, 2001.[8] Dean V Buonomano. Decoding temporal information: a model based on short-term synapticplasticity.
J. Neurosc. , 20(3):1129–1141, 2000.[9] Vittorio Caggiano, Roberto Leiras, Haizea Go˜ni-Erro, Debora Masini, Carmelo Bellardita,Julien Bouvier, V Caldeira, Gilberto Fisone, and Ole Kiehn. Midbrain circuits that set loco-motor speed and gait selection.
Nature , 553(7689):455–460, 2018.[10] Ronald L Calabrese. Half-center oscillators underlying rhythmic movements.
Nature , 261:146–148, 1995.[11] N Dale and A Roberts. Dual-component amino-acid-mediated synaptic potentials: excitatorydrive for swimming in xenopus embryos.
J. Physiol. , 363(1):35–59, 1985.
12] Simon M Danner, Natalia A Shevtsova, Alain Frigon, and Ilya A Rybak. Computational mod-eling of spinal circuits controlling limb coordination and gaits in quadrupeds.
Elife , 6:e31050,2017.[13] Alain Destexhe, Diego Contreras, Terrence J Sejnowski, and Mircea Steriade. A model ofspindle rhythmicity in the isolated thalamic reticular nucleus.
J. Neurophysiol. , 72(2):803–818, 1994.[14] Anca Doloc-Mihu and Ronald L Calabrese. A database of computational models of a half-center oscillator for analyzing how neuronal parameters influence network activity.
J. Biol.Phys. , 37(3):263–283, 2011.[15] Sten Grillner. Biological pattern generation: the cellular and computational logic of networksin motion.
Neuron , 52(5):751–766, 2006.[16] Charuni Gunaratne, Akira Sakurai, and Paul S. Katz. Variations on a theme: species differ-ences in synaptic connectivity do not predict central pattern generator activity.
J. Neurophys-iol. , 118:1123–1132, 2017.[17] Ronald M. Harris-Warrick and Jan-Marino Ramirez. Neural networks for the generation ofrhythmic motor behaviors. In
Neurobiology of Motor Control , chapter 8, pages 225–262. 2017.[18] Auke Jan Ijspeert. Central pattern generators for locomotion control in animals and robots:a review.
Neur. Netw. , 21(4):642–653, 2008.[19] Eugene M Izhikevich. Neural excitability, spiking and bursting.
Int. J. Bifurcat. Chaos ,10(06):1171–1266, 2000.[20] Sajiya Jalil, Dane Allen, Joseph Youker, and Andrey Shilnikov. Toward robust phase-lockingin melibe swim central pattern generator models.
Chaos , 23(4):046105, 2013.[21] Sajiya Jalil, Igor Belykh, and Andrey Shilnikov. Spikes matter for phase-locked bursting ininhibitory neurons.
Phys. Rev. E , 85(3):036214, 2012.[22] Ole Kiehn and Kimberly Dougherty. Locomotion: circuits and physiology. In
Neurosciencein the 21st Century , pages 1337–1365. Springer, 2016.[23] Maxime Lemieux, Nicolas Josset, Marie Roussel, S´ebastien Couraud, and Fr´ed´eric Bretzner.Speed-dependent modulation of the locomotor behavior in adult mice reveals attractor andtransitional gaits.
Front. Neurosci. , 10:42, 2016.[24] Matteo Lodi, Andrey Shilnikov, and Marco Storace. CEPAGEs: a toolbox for central patterngenerator analysis. In
Proc. IEEE ISCAS , pages 1–4, 2017.
25] Matteo Lodi, Andrey Shilnikov, and Marco Storace. Design of synthetic central patterngenerators producing desired quadruped gaits.
IEEE Trans. Circ. Syst. I , 65(3):1028–1039,2017.[26] Matteo Lodi, Andrey L Shilnikov, and Marco Storace. Design principles for central patterngenerators with preset rhythms.
IEEE Trans. Neural Netw. Learn. Syst. , 2020.[27] Yaroslav I Molkov, Bartholomew J Bacak, Adolfo E Talpalar, and Ilya A Rybak. Mechanismsof left-right coordination in mammalian locomotor pattern generation circuits: a mathematicalmodeling view.
PLoS Comp. Biol. , 11(5), 2015.[28] Roman Nagornov, Grigory Osipov, Maxim Komarov, Arkady Pikovsky, and Andrey Shilnikov.Mixed-mode synchronization between two inhibitory neurons with post-inhibitory rebound.
Commun. Nonlinear Sci. , 36:175–191, 2016.[29] M Pinco and A Lev-Tov. Synaptic transmission between ventrolateral funiculus axons and lum-bar motoneurons in the isolated spinal cord of the neonatal rat.
J. Neurophysiol. , 72(5):2406–2419, 1994.[30] Carla MA Pinto and Martin Golubitsky. Central pattern generators for bipedal locomotion.
J. Math. Biol. , 53(3):474–489, 2006.[31] Richard E Plant. Bifurcation and resonance in a model for bursting nerve cells.
J. Math.Biol. , 11(1):15–32, 1981.[32] Akira Sakurai and Paul S. Katz. Command or obey? homologous neurons differ in hierarchicalposition for the generation of homologous behaviors.
J. Neurosc. , 39(33):6460–6471, 2019.[33] David Somers and Nancy Kopell. Rapid synchronization through fast threshold modulation.
Biol. Cybern. , 68(5):393–407, 1993.[34] Kaoru Takakusaki. Neurophysiology of gait: from the spinal cord to the frontal lobe.
Mov.Disord. , 28(11):1483–1491, 2013.[35] X-J Wang. Fast burst firing and short-term synaptic plasticity: a model of neocortical chat-tering neurons.
Neurosc. , 89(2):347–362, 1999.[36] Jeremy Wojcik, Justus Schwabedal, Robert Clewley, and Andrey L Shilnikov. Key bifurcationsof bursting polyrhythms in 3-cell central pattern generators.
PloS One , 9(4), 2014.[37] Junzhi Yu, Min Tan, Jian Chen, and Jianwei Zhang. A survey on cpg-inspired control modelsand system implementation.
IEEE Trans. Neural Netw. Learn. Syst. , 25(3):441–456, 2013.
38] Le Zhao and Alain Nogaret. Experimental observation of multistability and dynamic attractorsin silicon central pattern generators.
Phys. Rev. E , 92(5):052910, 2015., 92(5):052910, 2015.