Mixed mode synchronization and network bursting of neurons with post-inhibitory rebound
Roman Nagornov, Grigory Osipov, Maxim Komarov, Arkady Pikovsky, Andrey Shilnikov
MMixed mode synchronization and network bursting of neurons withpost-inhibitory rebound
Roman Nagornov, a) Grigory Osipov, Maxim Komarov,
1, 2
Arkady Pikovsky,
1, 2 and Andrey Shilnikov
1, 3 Department of Control Theory, Nizhni Novgorod State University, Gagarin Av. 23, 606950, Nizhni Novgorod,Russia Department of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str 24/25, Potsdam,Germany Neuroscience Institute and Department of Mathematics and Statistics, Georgia State University, 100 Piedmont Str.,Atlanta, GA 30303, USA (Dated: 22 September 2018)
We study activity rhythms produced by a half-center oscillator – a pair of reciprocally coupled neurons thatare capable of producing post-inhibitory rebounds. Network (coupling-induced) burstings possessing twotime scales, one of fast spiking and another of slow bursting and quiescent periods, are shown to exhibitnontrivial synchronization properties. We consider several network configurations composed of both endoge-nously bursting, tonic spiking and quiescent neurons, as well as various mixed combinations. We show thatsynchronization at low frequency of bursting can be accompanied by complex, sometimes chaotic, patters offast spikes.
This study is focused on the mechanismsof rhythmogenesis and robustness of anti-phasebursting in half-center-oscillators (HCOs) con-sisting of two reciprocally inhibitory coupled neu-rons. There is a growing body of experimen-tal evidence that a HCO is a universal buildingblock for larger neural networks, including centralpattern generators (CPGs) controlling a varietyof locomotion behaviors in spineless animals andmammals. It remains unclear how CPGs achievethe level of robustness and stability observed innature. There has been a vastly growing consen-sus in the community of neurophysiologists andcomputational researchers that some basic struc-tural and functional elements are likely sharedby CPGs of both invertebrate and vertebrate an-imals. In this study we consider several config-urations of HCOs including coupled endogenousbursters, tonic spiking, and quiescent neurons,that become network bursters only when coupledby fast inhibitory synapses through the mecha-nism of post-inhibitory-rebound (PIR). The goalis better understanding the PIR mechanism as akey component for robust anti-phase bursting ingeneric HCOs.
I. INTRODUCTION
Synchronization of coupled oscillators is a fundamentalphenomenon in nonlinear systems that has been observedin a wide range of diverse applications . The mathe-matical concept of synchronization, first introduced anddeveloped for periodic oscillators , has further been gen-eralized for other aperiodic systems, including ones with a) Electronic mail: [email protected] chaotic dynamics. In life sciences, of a keen interest issynchronization or phase locking among oscillators withmultiple time scales. They may include mixed-mode andslow-fast relaxation-type oscillators , whose synergeticinteraction can give rise to the onset of a variety of syn-chronization patterns . In neuroscience, a plethora ofrhythmic motor behaviors with diverse time scales, suchas heartbeat, respiration, chewing, and locomotion onland and in water are produced and governed by neuralnetworks called Central Pattern Generators (CPGs) .A CPG is a microcircuit of neurons that is able to au-tonomously generate an array of polyrhythmic burstingpatterns, underlying various motor behaviors.Endogenous (self-sustained) bursting and network(coupling-induced) bursting are composite oscillatory be-haviors, featuring active phases during which a neuronor a group of neurons generates trains of fast action po-tentials, which are alternated with long interburst inter-vals during which it remains inactive or quiescent, untila new cycle of bursting occurs. In this paper we examinesynchronization of bursting patterns emerging throughinteractions of two interneurons coupled reciprocally byfast inhibitory synapses. This study has been driven bytwo major motivations: first, a general one concerningquestions on synchronization of mixed-mode oscillators.The second is a neuroscience related one, aimed at un-derstanding of intrinsic mechanisms of rhythmogenesisin CPGs composed, often symmetrically, of such smallnetworks of interneurons, as outlined below. It is stillunclear how CPGs can achieve the level of synergy androbustness to produce a plethora of rhythmic patternsobserved in nature.Recent experimental and theoretical studies have dis-closed a distinct role of CPGs in generation of adap-tive and coordinated motor activity of animals . Animportant feature of CPGs is their ability to producevarious types of rhythmic bursting activity, what causesflexible and adaptive locomotion of an organism. To ro- a r X i v : . [ n li n . C D ] D ec bustly govern motor patterns, CPGs are able to flexiblyadjust their oscillatory properties (such as bursts dura-tion, frequency of spiking, phase relations of bursts) dueto feedback from sensory inputs, for example, in responseto changes of environment . Up to a certain extent,the flexibility of CPG behaviors may be attributed to itsmultistability (of several coexistent attractors represent-ing different bursting rhythms in a phase space of thedynamical system) allowing for fast switching betweenoperating modes .From the theoretical point of view, a CPG is modeledas a small network of coupled oscillatory, or quiescent,interneurons, each described by a system of nonlinearordinary differential or difference equations (dynamicalsystem) . The study of CPGs allows one to developa general understanding of synchronization patterns inmixed-mode oscillators, applicable to systems of differentphysical and biological origin.There is a growing body of experimental evidence thata universal building block of most identified CPGs isa half-center oscillator (HCO) . A HCO is a pair oftwo reciprocally inhibitory interneurons bursting in al-ternation. Such a pair can be composed of endogenouslybursting interneurons, as well as of endogenously tonicspiking or quiescent ones that start bursting only whenthey are coupled in a network. Theoretical studies have indicated that formation of an anti-phase burst-ing rhythm is always based on some slow subsystem dy-namics. The slow subsystem is strongly associated withthe slow membrane current, such as persistent sodiumcurrent or slow calcium-dependent current; following we term currents associated with slow-varying concen-trations and gating variables as slow ones. There arethree basic mechanisms to generate alternating burstingin the HCO: release, escape, and post-inhibitory rebound(PIR). The first mechanism is typical for endogenouslybursting neurons . The other two mechanisms underlienetwork bursting in HCOs comprised of neurons whichare depolarized or hyperpolarized quiescent interneurons,respectively, in isolation The PIR mechanism uses reciprocal inhibition to keepnetwork bursting. It triggers an onset of a single or aseries of action potentials in the post-synaptic neuron af-ter it has been prolongedly hyperpolarized and abruptlyreleased from inhibition generated by the pre-synapticneuron during an active, tonic spiking phase of burst-ing. After that, the neurons of the HCO swap theiropposite roles to repeat the PIR mechanism. The PIRpromotes the action potential generation after a periodof sufficiently strong hyper-polarizing (inhibiting) input,as illustrated in Fig. 1. PIR is often caused by a low-threshold activated calcium current in neurons and theirbiophysically plausible models.This paper examine the PIR contribution to formation,synchronization and robustness of bursting rhythms ininhibitory coupled neurons. Here, we employ a modifica-tion of the Hodgkin-Huxley-type model introduced in ,to plausibly describe the PIR mechanism. Depending on its parameters, the model is known to produce an arrayof typical neuronal activities such as excitable dynamics,periodic spiking and bursting.Below we first examine the conditions that confirm andreproduce the PIR mechanism in the neurons. We willargue that the PIR is a pivotal component that promotesan alternating bursting rhythm in the HCO made of in-trinsically spiking and excitable neurons. Finally, we willshow that the PIR mechanism makes anti-phase networkbursting even more robust and easier to initiate in thecase of endogenously bursting neurons.The paper is organized as follows: in Section II we in-troduce the neuronal model, and discuss its dynamicalproperties in Section III. Section IV is focused on syn-chronization properties of a pair of endogenous bursters,while Section IV examines anti-phase bursting onset incoupled tonic-spiking neurons. Sections VI and VII dis-cusses the PIR-mechanisms of the HCO, and synergeticdynamics of mixed neurons: individually tonic spikingand hyperpolarized quiescent. V ( m V ) g Ca =1 . mS/cm g Ca =0 . mS/cm
900 1000 1200 2000 t ( ms ) I e x t ( µ A / c m ) FIG. 1. Abrupt release of a hyperpolarized pulse of the exter-nal current I ext (bottom panel) triggering a post-inhibitoryrebound of bursting in the quiescent neuron (top panel). Theparameter g Ca stated in the legend denotes maximal conduc-tance of the current I T . When the parameter g Ca is relativelysmall, the effect of the I T current is negligible what causesabsence of PIR in the dynamics of the neuron (dashed greencurve in the top panel). II. BASIC MODEL
In this study we use a neuronal model of the Hodgkin-Huxley type that has been introduced in . The basicequations describing the dynamics of the membrane po-tentials of the neurons read as follows (prime denotes thetime derivative, i = 1 , C m V (cid:48) i = I ( i ) ext − f i ( V i ) − N (cid:88) j =1 ,j (cid:54) = i I syn ( V i , V j ) ,f i ( V i ) = I ( i ) leak + I ( i ) Na + I ( i ) K + I ( i ) K [ Ca ] + I ( i ) T . (1)The variable V i ( t ) describes evolution of the membranepotential of the i -th neuron. The first two terms on theright-hand side of Eq. (1) govern intrinsic dynamics ofthe neuron: I ( i ) ext stands for a constant external currentapplied to the neuron, while term f i ( V i ) represents thesum of intrinsic ionic currents. The sum of synaptic cur-rents N (cid:80) j =1 ,j (cid:54) = i I syn ( V i , V j ) describes the coupling interac-tions between the neurons. The full details of currentrepresentation are given in the Appendix below.There are two principal control parameters in this sys-tem: the external current, I ( i ) ext , governs the activity typein an isolated, uncoupled neuron. The second parameter,the maximal conductance G of the synaptic current givenby Eq. (13), controls the coupling strength, and thus thenetwork dynamics of the coupled neurons.The low-threshold activated calcium current I T , mod-eled by Eq. (7) and regulated by the maximal conduc-tance g Ca , causes the PIR mechanism in the neurons (1),this effect is illustrated in Fig. 1. Rebound bursts canbe triggered by injection of a hyper-polarizing pulse ofexternal current I ( i ) ext (bottom panel in Fig. 1). If the im-pact of I T current is strong enough, the neuron generatesbursts of spikes, after recovering from the hyperpolarizedquiescent state (top panel in Fig. 1). III. DYNAMICS OF AN ISOLATED NEURON
Let us first consider the dynamics of an isolated neu-ron at I syn = 0 in Eq. (1). We treat the external cur-rent I ext as the primary bifurcation parameter in bothcomputational and biological experiments, which allowsone to examine and plausibly model various neuron ac-tivities and transitions between them. For the same pur-pose, we employ interspike intervals (ISIs), which are dis-tances between adjacent spikes generated by the neuron,as measurements characterizing the neuron’s dynamicalregimes. Figure 2 serves as an illustration of how the ISIschanges with variations of I ext .First we describe the neuron dynamics at large I ext values. For I ext > µA/cm the neuron produces tonicspiking activity (corresponding to a stable periodic orbitin the phase space of the model (1)). Here, the valueof ISI (which is the period of the stable orbit) is in-versely proportional to the I ext -value. The period de-creases with increase of I ext , and vice-versa. As one cansee from Fig. 2, near a critical value I ext ≈ µA/cm theISI shows an “unbounded” growth, which is an indica-tion of a bifurcation of the stable periodic orbit with an I ext ( µA/cm ) I S I ( m s ) FIG. 2. Bifurcation diagram for an isolated neuron: plot-ting inter-spike intervals (ISI) against the external current I ext (this horizontal axis is broken in the middle to high-light the ranges of nontrivial behaviors) reveals windowsof bursting ( − . µA/cm < ∼ I ext < ∼ . µA/cm ), quies-cence (0 . µA/cm < ∼ I ext < ∼ µA/cm ), and tonic spiking( I ext > ∼ µA/cm ) activity in the neuron. The low ISI branchcorresponds to short time intervals between fast spikes, andthe top branch represent long interburst intervals betweenconsecutive spike trains. arbitrarily long period. Detailed examinations indicatethe occurrence of a homoclinic saddle-node bifurcation underlying the transition from tonic spiking activity tohyperpolarized quiescence, represented by a stable equi-librium state at low values of the membrane potential V .This stable equilibrium state (a node with real and nega-tive characteristic exponents) persists within the param-eter window 0 . µA/cm < ∼ I ext < ∼ µA/cm , where theneuron remains ready for PIRs. Small perturbations ofthe quiescent state of the neuron have no pronounced ef-fects. Relatively strong perturbations can trigger a spike,after which the neuron comes back to the over-dampedquiescence state. With decreases of I ext from the thresh-old I ext ≈ µA/cm , the steady state becomes a focus.The stable focus loses the stability through a supercrit-ical Andronov-Hopf bifurcation at I ext ≈ . µA/cm .These oscillations are not seen in Fig. 2 because theyare under the threshold of spike generation. The ampli-tude of the oscillations increases rapidly as I ext is furtherdecreased. At I ext ≈ . µA/cm the shape of the pe-riodic orbit in the phase space changes via the additionof “turns, ” the so-called spike-adding mechanism .Such a periodic orbit is associated with bursting activ-ity, with turns corresponding to spikes within a burst.Bursting that includes two time scales can be recognizedin Fig. 2 through two characteristic branches: the bot-tom branch corresponds to small ISI values due to fastspiking, while the top branch is due to long interburstintervals. More details on geometry of bursting can befound in . Decreasing the external current I ext makesthe bursting neuron more depolarized and increases thenumber of spikes per burst. Having reached some maxi-mal value of spikes per burst at I ext ≈ − . µA/cm ,decreasing I ext gives rise to bursts with fewer spikes.The value I ext ≈ − . µA/cm corresponds to the oc-currence of another saddle-node bifurcation of equilibriasince the time interval between two consecutive burstsbecomes arbitrarily large. For I ext < − . µA/cm theneuron remains quiescent and excitable.The examination of bifurcations have been repeated forseveral values of the maximal conductance, g Ca , regulat-ing the low-threshold Ca current responsible for PIRin the neuron. We have found that this current bears aninsignificant affect on the intrinsic dynamics of the indi-vidual neurons, mainly because the membrane potentialdoes not decrease below the threshold value to activatethe current. As we will see below, this is not the casefor coupled neurons, where the parameter g Ca has a pro-nounced effect on the collective dynamics. IV. SYNCHRONIZATION OF TWO BURSTINGNEURONS
In this section and sections V-VII below, we explore arepertoire of rhythmic bursting outcomes generated by aHCO of two neurons, coupled reciprocally by fast, non-delayed inhibitory synapses. It is known that such cou-pled, endogenous bursters can produce a range of syn-chronous rhythmic outcomes with various fixed phase-lags due to overlapping spike interactions, see andreferences therein.We will first examine the cooperative dynamicsin the network of two coupled neurons ( N = 2in Eqs. (1)) with different I (1 , ext within the interval[ − . µA/cm , . µA/cm ]. This range of I ext corre-sponds to endogenous bursting in both neurons. Dueto the difference, ∆ I ext = I (2) ext − I (1) ext , the temporal andquantitative characteristics of endogenous bursters suchas their period, duty cycles, the spike numbers per burst,are different. The strength of coupling is quantified bythe maximal conductance, G , of the inhibitory synapticcurrent I syn . In what follows, we show that while ∆ I ext remains relatively small, increasing the coupling strengthshall give rise to an onset of synchrony between the neu-rons with the same bursting period. However, this willno longer be true for larger ∆ I ext values. A. Nearly identical bursters
While | ∆ I ext | < ∼ . µA/cm , synchronous burstingcan already occur at a relatively weak coupling. Here,the burst periods of both neurons become equal, how-ever the spike number per burst in the neurons may notbe the same. This is shown in Fig. 3, which depicts thedependence of the mean phase-lag between oscillationsof the low-threshold Ca -currents in the neurons, on G ( mS/cm ) ∆ ϕ g Ca =1 . mS/cm g Ca =1 . mS/cm g Ca =2 . mS/cm FIG. 3. Mean phase lag, ∆ ϕ plotted against the couplingstrength G for different values of the maximal conductance g Ca of the low-threshold Ca -current (larger g Ca values pro-mote stronger PIRs in the neurons); here I (1) ext = 0 . µA/cm and I (2) ext = 0 . µA/cm . the coupling strength G for several values of the con-ductance g Ca . Note that the post-inhibitory reboundsin the neurons under consideration are due to the slow-est, low-threshold Ca -current I T , whose magnitude iscontrolled by the maximal conductance g Ca . It is worthnoticing that PIRs occur more reliably with increasing g Ca . G ( mS/cm ) F r e q u e n c y ( k H z ) ∆ I =0 . µA/cm ∆ I =0 . µA/cm ∆ I =0 . µA/cm ∆ I =0 . µA/cm ∆ I =0 . µA/cm FIG. 4. Bifurcation diagram for antiphase synchronizationregimes: averaged frequency of bursting is plotted against thecoupling strength for several increasing ∆ I ext values. Doubleoverlapping branches are the indication of bursting dichotomywith two slow different frequencies in the two coupled neurons.Note a large plateau of the pronounced 4:3 frequency lockingat ∆ I ext = 0 . µA/cm , collapsing into anti-phase synchronylocked at an 1:1 ratio at higher G values. Nevertheless, for a broad range of initial condition, theneurons will predominantly burst in anti- phase ratherthan in-phase. This implies that the HCO can be bi-stable. This observation may shed light on properties ofmultifunctional CPGs composed of coupled neurons in-cluding coupled HCOs, which can switch between multi-ple functional states associated with different locomotiontypes. While in-phase bursting appears to be atypical forneurons with fast inhibitory coupling like in our case, thecoexistence of anti-phase and in-phase regimes have beenreported elsewhere . However, in-phase synchrony canprevail whenever both HCO neurons are driven exter-nally by another inhibiting burster or HCO .Next, we discuss the effects that the PIR activity ofthe neurons may have on antiphase bursting. Figure 3discloses how the coupling threshold for the onset of syn-chronous bursting depends on mismatches between theneurons of the HCO, as well as on the phase-locked fre-quency of slow I T -oscillations. The increasing phases,defined on modulo 2 π , of the neurons are reset at thevery beginning of each burst. Figure 3 shows the evolu-tion of the phase lag, ∆ ϕ , between the oscillations by theslow gating variables, h (1 , T , for inactivation of the low-threshold Ca -current, which is is plotted against thecoupling constant G for various g Ca -values. We foundin most cases that the phase lag tends eventually to π ,i.e. the anti-phase bursting occurs in the HCO, with thecoupling strength above a critical value. Depending on I ext , the actual ∆ ϕ -value may deviate from π . One cansee from Fig. 3 that increasing g Ca (promoting strongerpost-inhibitory rebound activity) leads to quite symmet-ric antiphase bursting even at small values of the couplingstrength. This is an explicit manifestation of the contri-bution of the slow low-threshold Ca -current in foster-ing the antiphase synchronization between the burstingneurons in the HCO. B. Bursters far from identical
Increasing ∆ I ext between coupled oscillators leads toan array of pronounced synchronization effects. Figure 4represents the slow (bursting) frequencies plotted againstcoupling strength, G , at several ∆ I ext values. One ob-serves the presence of overlapping branches at low G -values. This indicates that the coupled neurons generatebursting activities at different frequencies, until the cou-pling is increased over a threshold value at which bothbranches merge. This threshold value becomes higherwith increases of ∆ I ext , as the individual neurons be-come more and more different. We note from this figurethat for the largest case ∆ I ext = 0 . µA/cm , the HCOneurons become locked at the 4:3 frequency ratio, priorto the occurrence of the ultimate 1:1 frequency lockingat higher coupling values.Figure 5 illustrates transitions to synchrony for thelargest case ∆ I ext = 0 . µA/cm . The left panels demon-strate the established antiphase bursting in voltage tracesproduced by the HCO at various values of the coupling t ( ms ) V ( m V ) stneuron ndneuron h T h T stneuronspikes ndneuronspikes t ( ms ) V ( m V ) stneuron ndneuron h T h T stneuronspikes ndneuronspikes t ( ms ) V ( m V ) stneuron ndneuron h T h T stneuronspikes ndneuronspikes t ( ms ) V ( m V ) stneuron ndneuron h T h T stneuronspikes ndneuronspikes FIG. 5. Left column shows progressions of phaselocked voltage traces of two coupled neurons at∆ I ext = 0 . µA/cm , I (1) ext = 0 . µA/cm for G =0 . mS/cm ; 0 . mS/cm ; 0 . mS/cm ; 0 . mS/cm .Right panels demonstrate the Lissajous curves drawn byslow variables h (1) T and h (2) T , which correspond to (top)quasi-periodic dynamics at G = 0 .
02; a 4:3-frequency lockingregime at G = 0 . mS/cm ; an 1:1 chaotic locking at G = 0 . mS/cm ; (bottom) an 1:1 periodic locking at G = 0 . mS/cm . Circles and squares mark spike events inbursting. strength, G , below and above the threshold. The rightpanels represent the so-called Lissajous curves, which areparametrically traced down by the slow variables h (1) T and h (2) T of both neurons. These curves help to interpret thecorresponding types of frequency locking including quasi-periodic dynamics. On these curves we mark timing ofindividual spikes.It is worth noticing that such a large ∆ I ext makesthe phases of slow oscillations shift closer to π/
2, ratherthan to π . With a π/ G ≈ . mS/cm . More-over, our simulation indicated that there is a small win-dow, ∆ G ≈ . × − mS/cm , of hysteresis occur-ring at the transition. Figure 6 reveals a pecular fea-ture of the synchronous state, existing at strong cou-pling G ∈ [0 . mS/cm , . mS/cm ]: while thespike number in bursts generated by neuron 2 remainsnearly constant, the number of spikes in bursts by neu-ron 2 deviates quite strongly before it becomes fixed withstronger coupling.From the two bottom insets in Fig. 5 one can observethat depending on coupling strengths the 1:1 locking fea-tures chaotic and periodic modulations via the slow gat-ing variables of the neurons. In the chaotic case, thenumber of spikes per burst in either neuron or bothneurons will vary while the alternating bursting in theHCO persist locked at an 1:1-ratio within a windows[0 . mS/cm , . mS/cm ]. This is an explicit mani-festation of the phenomena of chaotic phase synchroniza-tion . G ( mS/cm ) s p i k e s / b u r s t st neuron nd neuron FIG. 6. Bifurcation diagram showing spike number per burstgenerated by the neurons plotted against the coupling con-stant G ; other parameters are same as in Fig. 5. Value G ≈ . mS/cm is a threshold towards the 1:1-frequencylocking; G ≈ . mS/cm corresponds to phase slipping inthe chaotic 1:1 locking state of alternating bursting. V. ANTIPHASE BURSTING OUT OF SPIKINGNEURONS
In this section we study emergence of alternating burst-ing in a HCO made of two coupled neurons that aretonically spiking in isolation at I (1) ext = 5 µA/cm withrelatively small ∆ I ext = 0 . µA/cm . Our goal is toreveal how increasing the inhibitory coupling strengthtransforms such tonic spikers into network bursters, asillustrated in Fig. 7.Weak coupling gives rise to periodic antiphase tonicspiking in the network at G = 0 . mS/cm (Fig. 7(a)).With further increases of the coupling strength, the HCO t ( ms ) V ( m V ) stneuron ndneuron (a) t ( ms ) V ( m V ) stneuron ndneuron (b) t ( ms ) V ( m V ) stneuron ndneuron (c) t ( ms ) V ( m V ) stneuron ndneuron (d) FIG. 7. Voltage oscillatory activity generated by the HCO at g Ca = 1 mS/cm and I ext = 5, and ∆ I ext = 0 . µA/cm :(a) antiphase spiking at G = 0 . mS/cm ; (b) chaotic spik-ing activity at G = 0 . mS/cm ; (c) forced sub-thresholdoscillations in neuron 2 due to fast spiking in neuron 1 at G = 2 mS/cm , (d) network bursting at G = 4 mS/cm . neurons evolve through a window of irregular, unreliablespiking at G = 0 . mS/cm , towards an asymmetric net-work dynamics. In this robust regime, either neurongenerates tonic spiking activity, while the other one isforced to produce sub-threshold oscillations of quite alarge amplitude, such as ones shown in Figs. 7(b,c) at at G = 2 mS/cm . A dramatic increase in the inhibitorycoupling strength beyond G ≈ mS/cm finally forcesthe neurons to begin bursting in alternation.Such bursting activity is often referred to as a networkbursting, which is the result of strongly reciprocal inter-actions of two individually tonic spiking neurons. Fig-ure 2 provides an explanation why inhibitory couplingmust be so strong to induce network bursting in the givenmodel. There is a wide gap between the I ext parametervalues corresponding to bursting (at the low end) and totonic spiking (on the high end). Within this gap, theneurons remain hyperpolarized quiescent, and hence thenetwork bursting can only occur through PIR mecha-nisms. A strong flux of inhibitory current is requiredto originate from the presynaptic, tonic spiking neuronin order to make the postsynaptic neuron drive over thegap to become a network burster. Then, the neuronsswap the roles of the driving and driven neurons and soforth. This assertion is further supported by Fig. 8 plot-ting inter-spike intervals (ISIs) in the driven postsynap-tic neurons (neurons 1 and 2 in Insets (a) and (b), resp.)against the coupling strength for three values of g Ca . Wecan see the corresponding gap in the right panel, withinwhich neuron 2, through strongly inhibitory force, gen-erates large sub-threshold oscillations, partially becauseof the network asymmetry due to ∆ I ext = 0 . µA/cm .Figure 8 also discloses the qualitative role of the slowlow-threshold Ca -current I T ; namely, increasing thecorresponding rebound parameter g Ca lowers the thresh-old of coupling-induced bursts, and shrinks the parame-ter interval of hyperpolarized quiescence in the neurons. G ( mS/cm ) I S I ( m s ) g Ca =1 mS/cm g Ca =1 . mS/cm g Ca =2 . mS/cm (a) G ( mS/cm ) I S I ( m s ) g Ca =1 mS/cm g Ca =1 . mS/cm g Ca =2 . mS/cm (b) FIG. 8. Bifurcation diagram of inter-spike intervals for cou-pled neuron 1 (a) and neuron 2 (b) plotted against the cou-pling constant G for three increasing values of g Ca ; here I ext = 5 µA/cm and ∆ I ext = 0 . µA/cm . VI. PIR MECHANISM FOR ANTIPHASE BURSTING G ( mS/cm ) I S I ( m s ) g Ca =1 . mS/cm g Ca =1 . mS/cm g Ca =2 . mS/cm (a) G ( mS/cm ) I S I ( m s ) g Ca =1 . mS/cm g Ca =1 . mS/cm g Ca =2 . mS/cm (b) FIG. 9. Bifurcation diagram of ISIs evolutions for coupledneuron 1 (a) and neuron 2 (b) plotted against the couplingconstant G ; here I (1) ext = 2 µA/cm and ∆ I ext = 0 . µA/cm correspond to hyperpolarized quiescent neurons in isolation. In this section we examine the emergence of networkantiphase bursting through the PIR mechanism. Here, itis imperative that both neurons remain hyperpolarizedwithin a parameter window quiescence 0 . µA/cm < ∼ I ext < ∼ µA/cm , as seen from Fig. 2. It is also nec-essary for PIR network bursting that the initial statesof the neurons must be different: one tonic spiking andone hyperpolarized quiescent. An alternative is an ap-plication of negative pulse of current to trigger a PIR ina targeted neuron. In addition to the above constrains,the coupling strength must exceed a certain threshold,see Fig. 9.It is seen from this plot that relatively weak couplingcannot initiate PIR network bursting. Rather some sub-threshold oscillations are generated in the post-synapticneuron, which cease as soon as the pre-synaptic neu-ron ends its active spiking phase and becomes hyper-polarized quiescent as well. For each set of the parame-ters there is a threshold for the coupling strength beyond which the HCO robustly produces anti-phase bursting.Figure 9 also demonstrates that the PIR mechanism ofnetwork bursting becomes more reliable with increasing g Ca , which lowers the threshold value of inhibitory cou-pling.Regularization of PIR network bursting is quite sensi-tive to variations of the coupling strength. By inspectingthe top ISI-branch in Fig. 9 one can identify the stabil-ity windows of G -values corresponding to plateau-like ISIintervals alternating with the windows of fluctuations ofISI values. Fluctuations are not marginal – about 10%off the mean value. The PIR network bursting is illus-trated in Fig. 10(a). Fig. 10(b) reveals a quasi-periodicmodulation of the slow gating variables; here the periodicLissajous curve has a shape of figure-eight. t ( ms ) V ( m V ) stneuron ndneuron (a) h T h T stneuronspikes ndneuronspikes (b) FIG. 10. (a) Antiphase network bursting through the PIRmechanism is filled out densely by alternating voltage traceswithout quiescent gaps unlike the case of coupled endoge-nous bursters through the release mechanism. (b) Lissajouscurve of a figure-eight shape traced down by the slow gat-ing variables of the models at G = 1 . mS/cm and g Ca =1 . mS/cm ; other parameters the same as in Fig. 9). Circlesand squares in panel (b) mark occurrences of spikes. VII. COUPLING OF NEURONS IN DIFFERENT MODES
In this section we consider the network dynamics ofcoupled neurons operating in different modes. Specifi-cally, the neurons are set close to “quiescence-spiking”and “quiescence-bursting” transitions, respectively. AsFig. 2 suggests, we choose the following values of the ex-ternal drive: I (1) ext = 4 . µA/cm and I (2) ext = 3 . µA/cm .This means that in isolation neuron 1 produces tonicspiking activity, while neuron 2 is hyperpolarized andready for the PIR mechanism. Our network simula-tions of two such neurons are summarized in Figs. 11and 12. We can conclude that while coupling remainsweak, neuron 1 tonically spikes (Fig. 12(a)). Abovea threshold, G ≈ . mS/cm , neuron 2 begins spik-ing as well. Fig. 12(b)) depicts alternating spiking at G = 1 . mS/cm with the locking ratio 2:1, which meanstwo spikes in neuron 1 are followed by a single spike inneuron 2. The network begins to burst in antiphase,first with irregularly varying spike numbers, as G exceeds2 mS/cm , as illustrated in Fig. 12(c). Further increasein the coupling strength above G = 2 . mS/cm , regu-larizes network bursting that becomes stable and similarto the pattern shown in Fig. 10. G ( mS/cm ) I S I ( m s ) st neuron nd neuron FIG. 11. Bifurcation diagram of ISI plotted against the cou-pling parameter G indicates the threshold beyond which thenetwork made of the tonically spiking neuron 1 and the hy-perpolarized quiescent neuron 2 begins generating anti-phasebursting; g Ca = 1 . mS/cm . t ( ms ) V ( m V ) stneuron ndneuron (a) t ( ms ) V ( m V ) stneuron ndneuron (b) t ( ms ) V ( m V ) stneuron ndneuron (c) FIG. 12. Established voltages traces generated by the HCOneurons: (a) tonic spiking in neuron 1 and inhibition inducedsubthreshold oscillations in neuron 2 at G = 1 mS/cm ; (b)periodic antiphase spiking with a 2:1 locked ratio at G =1 . mS/cm ; (c) chaotic antiphase bursting near the thresholdat G = 2 mS/cm . Figure 13 represents the bifurcation diagram for thenetwork of an endogenous burster – neuron 2 at I (2) ext =0 . µA/cm , coupled with neuron 1 – hyperpolarized at I (1) ext = 1 µA/cm . Because the quiescent neuron 1 ini-tially remains below the synaptic threshold, anti-phasebursting in the network may only start when the PIRmechanism is induced by the endogenous burster. Thisexplains a qualitative resemblance of the diagram inFig. 13 with that in Fig. 11 for the HCOP made of tonic-spiking and quiescent neurons. G ( mS/cm ) I S I ( m s ) st neuron nd neuron FIG. 13. Bifurcation diagram of ISI plotted against the cou-pling parameter G indicates the threshold beyond which thenetwork made of the bursting neuron 2 and the hyperpolar-ized quiescent neuron 1 begins generating anti-phase bursting; g Ca = 1 . mS/cm . VIII. CONCLUSIONS
This study focused on the mechanisms of rhythmogen-esis of anti-phase bursting in HCO networks consisting oftwo reciprocally inhibitory coupled neurons. Such HCOsare primary building blocks for larger neural networks,including CPGs controlling a plethora of locomotion be-haviors in spineless animals and mammals. There is agrowing consensus in the neuroscience community thatCPGs share some universal principles of their function-ing.In this study, for CPG neurons we used the biophys-ically plausible Hodgkin-Huxley-type model. Its mainfeature is post-inhibitory rebound dynamics: a quiescentneuron is able to produce a single or a series of spikesafter it has been quickly released from inhibition by an-other pre-synaptic neuron, or by a hyperpolarized pulseof external current. The PIR mechanism allows a pairof such neurons to stably generate anti-phase burstingin a network initiated externally by inhibitory perturba-tion(s).We considered several configurations of HCOs includ-ing coupled endogenous bursters. We also discussedHCOs involving tonic spiking and quiescent neurons thatbecome network bursters only when coupled by fast in-hibitory synapses. In our examination of synchronizationproperties of bursting we found that, in all consideredcases, the network can reliable achieve synchrony in anti-phase bursting. We also described some partial configu-rations leading to incomplete synchronization where theneurons become partially synchronized, through slow-varying currents. Meanwhile, cross-correlations in theirfast voltage dynamics are not always obvious. These fastvoltage dynamics may give rise to the emergence of syn-ergistically complex states, including chaos in the neuralnetwork.We found that while enhancing the PIR mechanismdoes not lead to drastic changes in the dynamics of theindividual neurons, it causes significant modulation ofnetwork dynamics. Specifically, we found that the win-dows of anti-phase bursting rhythms can be extendedin the parameter space of the network, when increasingthe PIR mechanisms in individual neurons become moreprominent. This suggests that the PIR is a key compo-nent for robust and stable anti-phase bursting in HCOs.In the future, we plan to examine specific networks con-sisting of several HCOs, which have been identified inswim CPGs of several sea mollusks.
ACKNOWLEDGMENTS
The research is supported by the grant (the agreementof August 27, 2013 N 02.B.49.21.0003 between The Min-istry of education and science of the Russian Federationand Lobachevsky State University of Nizhni Novgorod,Sections 2-4) and by the Russian Science Foundation(Project No. 14-12-00811, Sections 5-7). A.P. thanksA. Politi for interesting discussions and acknowledges theG.Galilei Institute for Theoretical Physics (Italy) for thehospitality and the INFN for partial support during thecompletion of this work. M.K. thanks Alexander vonHumboldt foundation for support. A.S. also acknowl-edges the support from NSF grants DMS-1009591, RFFI11-01-00001, and RSF grant 14-41-00044. We thank A.Kelley for helpful suggestions.
IX. APPENDIX: CONDUCTANCE BASED MODEL
The model in this study is adopted from Ref. . Thedynamics of the membrane potential, V is governed bythe following equation: C m V (cid:48) = I ext − I T − I leak − I Na − I K − I K [ Ca ] − I syn ; (2)here, C m is the specific membrane capacity, I ext is the ex-ternal current in nA , I T is the slow low-threshold Ca -current, I leak is the leakage current, I Na is the N a + -current, I K is the K + -current, I K [ Ca ] is the slow Ca -activated K + -current which moderates bursting activityin the model, and I syn is the synaptic current from otherneurons. Leak current I leak is given by I leak = g Na ( V − E leak ) (3)with E leak = − mV being reversal potential for leakcurrent, and maximal conductance g L = 0 . mS/cm .Dynamics of fast N a + -current I Na = g Na m h ( V − E Na )is described by the following equations: m (cid:48) = 0 . − V ) e . − V ) − − m ) − . V − e . V − − mh (cid:48) = 0 . · e − V (1 − h ) − e − . V − + 1 h, (4)where g Na = 100 mS/cm is the maximal conductance of N a + -current, E Na = 50 mV is the Nernst equation po-tential for N a + -current, m and h are the gating variablesdescribing activation and inactivation of the current. Dy-namics of fast K + -current I K is described by I K = g K n ( V − E K ) ,n (cid:48) = 0 . − V ) e . − V ) − − n ) − . e (10 − V )40 n (5)with n being the gating activation variable; here E K = − mV and g K = 10 mS/cm . Dynamics of intraneu-ronular concentration of calcium ions [ Ca ] is described[ Ca ] (cid:48) = − kI T F d − K T [ Ca ][ Ca ] + K d , (6)where the first term is a inflow through thin membranedue low-threshold Ca -current, and the second is a con-tribution of Ca ion-pump, F = 96 , C/mol is theFaraday constant, d = 1 µm is the membrane thick-ness, k = 0 . I T ] = µA/cm and [ Ca ], that has dimension of millimole , K T = 10 − mM/ms − and K d = 10 − mM are the speedconstants from the Michaelis–Menten approximation ofion-pump kinetics, they are taken the way that one canobserve fast output of calcium from a neuron in millisec-ond time scale.Dynamics of slow low-threshold Ca -current I T is mod-eled as follows: I T = g Ca m T h T i ( V − E Ca ) , (7)with g Ca = 1 . mS/cm , and the gating activation m and inactivation h are given by m (cid:48) = − m − m ∞ ( V ) τ m ( V ) , h (cid:48) = − h − h ∞ ( V ) τ h ( V ) , (8)with m ∞ ( V ) = 11 + e − V +527 . ,τ m ( V ) = 0 .
44 + 0 . e V +2710 + e − V +10215 ) ,h ∞ ( V ) = 11 + e V +805 ,τ m ( V ) = 22 . . e V +484 + e − V +40750 ) . (9)These constants were taken at temperature 36 o C andextraneuronular concentration calcium concentration[ Ca ] = 2 mM . Equilibrium potential E Ca dependingon intraneuronular concentration of Ca are found fromthe Nernst equation: E Ca = ¯ k R · T F ln (cid:18) [ Ca ] [ Ca ] (cid:19) ; (10)here R = 8 . J K/mol o , T = 309 . o K , dimension-less constant ¯ k = 1000 for E Ca is measured in milli-volts, extraneuronular concentration of calcium ions is[ Ca ] = 2 mM .Dynamics of Ca -activated K + -current is determinedthrough the following equations: I K [ Ca ] = g K [ Ca ] m K [ Ca ] ( V i − E K ) , where m (cid:48) = − m − m ∞ ([ Ca ]) τ m ([ Ca ]) , (11)and τ m ([ Ca ]) = 148[ Ca ] + 0 . ,m ∞ ([ Ca ]) = 48[ Ca ] Ca ] + 0 . , (12)with g K [ Ca ] = 10 mS/cm , E K = − mV m being thegating activation variable.Synaptic currents are models using the fast-thresholdmodulation paradigm : I syn ( V i , V j ) = GS ( V j − θ syn )( V i − E syn ) , (13)0where G is the maximal conductance of synaptic cur-rent flowing from pre-synaptic j-th neuron into the post-synaptic i-th neuron. For inhibitory coupling we set E syn = − mV ; the synaptic activity function S ( V ( j ))is given by S ( V ) = 11 + e − V − θ ) , (14)with the synaptic threshold θ = 20 mV set in a middleof fast spikes. A. Pikovsky, M. Rosenblum, and J. Kurths,
Synchronization: AUniversal Concept in Nonlinear Sciences , Cambridge NonlinearScience Series (Cambridge University Press, 2003). A. Shilnikov, L. Shilnikov, and D. Turaev, “On some mathemat-ical topics in classical synchronization: a tutorial,” J Bifurcationsand Chaos , 2143–2160 (2004). M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. M.Osinga, and M. Wechselberger, “Mixed-mode oscillations withmultiple time scales,” SIAM Rev. , 211–288 (2012). S. Jalil, I. Belykh, and A. Shilnikov, “Fast reciprocal inhibitioncan synchronize bursting neurons,” Physical Review E (2010),10.1103/PhysRevE.81.045201. I. Omelchenko, M. Rosenblum, and A. Pikovsky, “Synchroniza-tion of slow-fast systems,” European Physical Journal-SpecialTopics , 3–14 (2010). S. Jalil, I. Belykh, and A. Shilnikov, “Spikes matter for phase-locked bursting in inhibitory neurons,” Physical Review E (2012), 10.1103/PhysRevE.85.036214. J. Wojcik, J. Schwabedal, R. Clewley, and A. Shilnikov, “Keybifurcation of bursting polyrhythms in 3-cell central pattern gen-erators,” PLoS ONE , e92918 (2014). A. Selverston, ed.,
Model Neural Networks and Behavior (Springer, Berlin, 1985). E. Marder and R. Calabrese, “Principles of rhythmic motor pat-tern generation,” Physiological Reviews , 687–717 (1996). S. Gillner and P. Wallen, “Central pattern generators for locomo-tion, with special references to vertebrates,” Ann. Rev. Neurosci. , 233–261 (1985). M. I. Rabinovich, P. Varona, A. I. Selverston, and H. D. I.Abarbanel, “Dynamical principles in neuroscience,” Rev. Mod.Phys. , 1213–1265 (2006). K. M. Shaw, D. N. Lyttle, J. P. Gill, M. J. Cullins, J. M. Mc-Manus, H. Lu, P. J. Thomas, and H. J. Chiel, “The significanceof dynamical architecture for adaptive responses to mechanicalloads during rhythmic behavior,” Journal of Computational Neu-roscience , 1–27 (2014). J. T. C. Schwabedal, A. B. Neiman, and A. L. Shilnikov, “Robustdesign of polyrhythmic neural circuits,” Phys. Rev. E , 022715(2014). E. M. Izhikevich,
Dynamical Systems in Neuroscience (MITPress, Cambridge, mass., 2007). A. A. V. Hill, S. Van Hooser, and R. L. Calabrese, “Half-centeroscillators underlying rhythmic movements,” in
The Handbook ofBrain Theory and Neural Networks , edited by M. A. Arbib (TheMIT Press, 2003). X.-J. Wang and J. Rinzel, “Alternating and synchronous rhythmsin reciprocally inhibitory model neurons,” Ann. Rev. Neurosci. , 233–261 (1985). S. Daun, J. E. Rubin, and I. A. Rybak, “Control of oscillationperiods and phase durations in half-center central pattern gener-ators: a comparative mechanistic analysis,” Journal of Compu-tational Neuroscience , 3–36 (2009). A. Destexhe, D. Contreras, T. Sejnowski, and M. Steriade, “AModel of Spindle Rhythmicity in the Isolated Thalamic ReticularNucleus,” Journal of Neurophysiology , 803–818 (1994). D. Perkel and B. Mulloney, “Motor Pattern Production in Recip-rocally Inhibitory Neurons Exhibiting Postinhibitory Rebound,”Science , 181–183 (1974). F. Skinner, N. Kopell, and E. Marder, “Mechanisms for oscilla-tion and frequency control in reciprocally inhibitory model neu-ral networks,” Journal of Computational Neuroscience , 69–87(1994). J. Angstadt, J. Grassmann, K. Theriault, and S. Levasseur,“Mechanisms of postinhibitory rebound and its modulation byserotonin in excitatory swim motor neurons of the medicinalleech,” Journal of Comparative Physiology A-Neuroethology,Sensory, Neural and Behavioral Physiology , 715–732 (2005). N. Kopell and G. Ermentrout, “Mechanisms of phase-locking andfrequency control in pairs of coupled neural oscillators,” Hand-book on Dynamical Systems , 3–54 (2002). L. Shilnikov, A. Shilnikov, D. Turaev, and L. Chua,
Methodsof Qualitative Theory in Nonlinear Dynamics, Parts I and II (World Scientific Publ., 1998,2001). P. Channell, G. Cymbalyuk, and A. Shilnikov, “Origin ofbursting through homoclinic spike adding in a neuron model,”Phys.Rev.Lett. , 134101 (2007). P. Channell, I. Fuwape, A. Neiman, and A. Shilnikov, “Vari-ability of bursting patterns in a neuron model in the presence ofnoise.” Journal Computational Neuroscience , 527–542 (2009). A. Shilnikov, “Complete dynamical analysis of a neuron model,”Nonlinear Dynamics , 305–328 (2012). S. Jalil, D. Allen, J. Youker, and A. Shilnikov, “Toward robustphase-locking in melibe swim central pattern generator models,”Chaos , 046105 (2013). I. Belykh and A. Shilnikov, “When weak inhibition synchronizesstrongly desynchronizing networks of bursting neurons.” PhysRev Lett , 078102 (2008). A. Shilnikov, R. Gordon, and I. Belykh, “Polyrhythmic syn-chronization in bursting networking motifs,” Chaos (2008),10.1063/1.2959850. N. Kopell and D. Somers, “Rapid synchronization through fastthreshold modulation,” Biol. Cybern.68