Table Of ContentModeling of Spiking-Bursting Neural Behavior Using Two-Dimensional Map
Nikolai F. Rulkov
Institute for Nonlinear Science, University of California, San Diego, La Jolla, CA 92093-0402
(February 8, 2008)
A simple model that replicates the dynamics of spiking and spiking-bursting activity of real bio-
logical neurons is proposed. The model is a two-dimensional map which contains one fast and one
slow variable. The mechanisms behind generation of spikes, bursts of spikes, and restructuring of
the map behavior are explained using phase portrait analysis. The dynamics of two coupled maps
which modelthebehavioroftwoelectrically coupled neuronsisdiscussed. Synchronizationregimes
for spiking and bursting activity of these maps are studied as a function of coupling strength. It is
demonstrated that the results of this model are in agreement with the synchronization of chaotic
2 spiking-bursting behaviorexperimentally found in real biological neurons.
0
0 PACS number(s): 05.45.+b, 87.22.-q
2
n
a I. INTRODUCTION chaotic spiking-bursting oscillations the system of dif-
J ferential equations describing each neuron should be at
5 Understanding dynamical principles and mechanisms least a three-dimensional system (see for example [7]).
behindthecontrolofactivity,signalandinformationpro- Further simplification of phenomenological models for
]
D cessing that occur in neurobiological networks is hardly complex dynamics of the neuron can be obtain using
C possiblewithoutnumericalstudiesofcollectivedynamics dynamical systems in the form of a map. Despite the
of large networks of neurons. These simulations need to low-dimensional phase space of this nonlinear map, it is
.
n take into account the complex architecture of couplings abletodemonstratealargevarietyofcomplexdynamical
i
l amongindividualneuronsthatissuggestedbydatafrom regimes.
n
biological experiments. One of the complicating factors The use of low-dimensional model maps can be use-
[
in understanding the simulation results is the complex- ful for understanding the dynamical mechanisms if they
1 ityoftemporalbehaviorofindividualbiologicalneurons. mimic the dynamics of oscillations observed in real neu-
v
This complexity is due to the large number of ionic cur- rons, show correct restructuring of collective behavior,
6
0 rents involved in the nonlinear dynamics of neurons. As and are simple enough to study the reasons behind such
0 a result, realistic channel-based models proposed for a restructuring. Thedevelopmentofamodelmapwhichis
1 single neuron is usually a system of many nonlinear dif- capable of describing various types of neuralactivity, in-
0 ferential equations (see, for example [1–5] and the re- cluding generation of tonic spikes, irregular spiking, and
2 view of the models in [6]). The strong nonlinearity and both regularand irregularbursts of spikes, is the goalof
0
/ high dimensionality of the phase space is a significant this paper.
n obstacle inunderstanding the collective behaviorof such Itisknownthatconstructingalow-dimensionalsystem
i
l dynamical systems. The dynamical mechanisms respon- ofdifferentialequationwhichiscapableofgeneratingfast
n
sible for restructuring collective behavior of a network spikes bursts excited on top of the slow oscillations, one
:
v of channel-based neuron models are difficult or impossi- needs to consider a system which has both slow and fast
Xi ble to analyze because these mechanisms are well hid- dynamics (see for example [7–11]). Using the same ap-
den behind the complexity of the equations. If there is proach one can construct a two-dimensional map, which
r
a a chance to identify possible dynamical mechanisms be- can be written in the form
hindaparticularbehaviorofthenetworkwithouttheuse
x =f(x ,y +β ) , (1a)
of complex models, then this knowledge can be used to n+1 n n n
guide through numerical study of this behavior with the yn+1 =yn µ(xn+1)+µσn , (1b)
−
high-dimensional channel-based models. Based on these
where x is the fast and y is the slow dynamical vari-
results one can specify the actual dynamical mechanism n n
able. Slow time evolution of y is due to small values of
occurred in the network and understand the biological n
the parameter µ=0.001. Terms β and σ describe ex-
processes which contribute to it. n n
ternalinfluences appliedto the map. Theseterms model
One way to reveal the dynamical mechanisms is to
thedynamicsoftheneuronundertheactionoftheexter-
use a simplified or phenomenological model of the neu-
nal DC bias current I and synaptic inputs Isyn. The
ron. However, many experiments indicate the restruc- DC n
term σ can also be used as the control parameter to
turing of collective behavior utilizes a variety of dy- n
select the regime of individual behavior.
namical regimes generated by individual dynamics of
A model in the form of two-dimensional map, whose
a neuron. When neural dynamics include the regime
form is similar to (1) was used in the study of dynami-
1
calmechanismsbehindtheemergenceandregularization this section is to illustrate dynamical mechanisms that
of chaotic bursts in a group of synchronously bursting control the response, and to show how the relation be-
cells coupled through a mean field [12,13]. The effect of tweenthe two inputs ofthe mapinfluence the properties
anti-phase regularization was also modeled before with of the response. Section IV presents the results ofa syn-
one-dimensionalmaps [14]. In both these models the os- chronization study in two coupled maps.
cillations during the burst were described by a chaotic
trajectory. Map (1) improves these models by adding a
II. INDIVIDUAL DYNAMICS OF THE MAP
feature that enables one to mimic the dynamics of in-
dividual spikes within the burst. This is achieved us-
ing a modification of the shape of the nonlinear function Consider the regimes of oscillations produced by the
f(x,y)whichisnowadiscontinuousfunctionoftheform individual dynamics of the map. In this case the inputs
β = β and σ = σ are constants. Note that if β is a
n n
α/(1 x)+y, x 0
constant, then it can be omitted from the equations us-
− ≤
f(x,y)=α+y, 0<x<α+y (2) ingthe changeofvariabley +β ynew. Therefore,the
1, x α+y n → n
individualdynamicsofthemapdependsonlyonthecon-
− ≥
trol parameters α and σ, and the map can be rewritten
where α is a control parameter of the map. The depen-
in the form
dence of f(x,y) on x computed for a fixed value of y is
shown in Fig.1. In this plot the values of α and y are
x =f(x ,y ), (3a)
set to illustrate the possibility of coexistence of limit cy- n+1 n n
y =y µ(x +1)+µσ . (3b)
cle, Pk, corresponding to spiking oscillation in (1a), and n+1 n− n
fixed points x and x . Note that when y increases or
ps pu Typical regimes of temporal behavior of the map are
decreasesthe graphoff(x,y)movesupordown,respec-
shown in Figs 2 and 3. When the value of α is less then
tively, except for the third interval x α+y, where the
≥ 4.0then,dependingonthevalueofparameterσ,themap
values of f(x,y) always remain equal to -1.
generatesspikesorstaysinasteadystate(seeFig2). The
frequencyofthespikesincreasesasthevalueofparameter
σ is increased (see Fig 2).
2
(a)
1
1
n
x 0
1
+ 0
n P −1
x
k
−1 (b)
1
x
pu n
−2 x x 0
ps
−1
−2 −1 0 1 2
(c)
x
n 1
n
x 0
FIG.1. The shape of the nonlinear function f(x,y), plot-
ted for α = 6.0 and y = −3.93, is shown with a solid line.
−1
Thedashedlineillustrates asuper-stablecycle Pk ofthefast
map (3a), where the value of y is fixed. The stable and un- 0 500 1000 1500 2000
n
stable fixed points of the map are indicated by xps and xpu,
respectively.
FIG.2. Waveforms of spiking behavior generated by map
The paper is organizedasfollows: SectionII considers (3) with α=4. Panel (a) shows the transition to the regime
the features ofthe fastand slowdynamics which explain of silence for σ = −0.01. The regimes of continuous tonic
the formationofvarioustypesofbehaviorintheisolated spiking are computed for σ=0.01 (b) and σ=0.1 (c).
map. Italsodiscussesthebifurcationsresponsibleforthe
qualitative change of activity. The dynamics of response For α>4 the map dynamics are capable of producing
generated by the map after it is excited or inhibited by bursts of spikes. The spiking-bursting regimes are found
an external pulse is discussed in Section III. The goalof in the intermediate region of parameter σ between the
2
regimes of continuous tonic spiking and steady state (si- Considering the fast and slow dynamics together, one
lence). The spiking-bursting regimes include both peri- can see that, if x is in the stable branch S (y), then
s ps
odicandchaoticbursting. Afewtypicalbursting-spiking the map (3) has a stable fixed point. The stable fixed
regimes computed for different values of σ are presented point corresponds to the regime of silence in the neu-
in Fig 3. ral dynamics. The oscillations in the map dynamics will
Due to the two different time scales involved in the appear when x >1 √α. This is the threshold of exci-
s
−
dynamics of the model, the mechanisms behind the re- tation, which corresponds to the bifurcation values of σ
structuringofthedynamicalbehaviorcanbeunderstood given by
throughthe analysisof the fast andslowdynamics sepa-
σ =2 √α. (6)
rately. Inthiscase,timeevolutionofthefastvariable,x, th −
is studied with the one-dimensional map (3a) where the Tounderstandthedynamicsofexcitedneuronswithin
slow variable, y, is treated as a controlparameter whose this model one needs to consider branches that corre-
value drifts slowly in accordance with equation (3b). It spond to the spiking regime. To evaluate the location
follows from (3b) that the value of y remains unchanged of the spiking branches, consider the mean value of x
n
only if x=xs given by computedfortheperiodictrajectoryofthefastmap(3a)
as a function ofy. Here y is treatedas a parameter. Use
x = 1+σ. (4)
s − of such an approximationis quite typical for an analysis
involvingfastandslowdynamics. Itworkswellforsmall
If x<x , thenthe value of y slowlyincreases. If x>x ,
s s values of parameter µ.
then y decreases.
It follows from the shape of f(x,y) that for any value
of y, map (3a) generates no more than one periodic tra-
2 (a) jectory P , where k is the period (i.e. x = x ). The
k n n+k
1 cycles P always contain the point x = 1 and the
k n
xn 0 index k increases stepwise (k k +1) as −y decreases.
−1 Since the trajectoryP always→has a point in the flat in-
k
−2 terval of f(x,y), all these cycles are super-stable except
3 (b) for bifurcation values of y for which the trajectory con-
2
tains the point x = 0 (see Fig 1). The location of the
n
1
n spiking branch, S , of ”slow” motion in the phase
x 0 spikes
plane (x ,y ) can be estimated as the mean value of x
−1 n n
−2 computed for the period of cycle Pk,
3 (c)
k
2 1
x = f(m)( 1,y), (7)
n 1 mean k X −
x 0 m=1
−1
where k is the period of P , and f(m)(x,y) is the m-th
−2 k
iterateof(3a),startedatpointxandcomputedforfixed
0 500 1000 1500 2000
n y. Thespikingbranchof”slow”dynamicsevaluatedwith
(7) is shown in Fig4. One can see that this branch has
many discontinuities caused by the bifurcations of the
FIG.3. Typicalwaveformsofthespiking-burstingbehavior
super-stable cycles P .
k
generated bythemap (3) forthefollowing parameter values:
To complete the picture of fast and slow dynamics of
α = 4.5, σ = 0.14 (a); α = 6.0, σ = −0.1 (b); and α = 6.0,
the model for α>4, one needs to consider the fast map
σ=0.386 (c).
bifurcation associated with the formation of homoclinic
orbit h originating from the unstable fixed point, x .
From(3a)onecanfindtheequationforthecoordinate pu pu
This homoclinic orbitoccurs when the coordinate of x
of fixed points x of the fast map. This equation is of pu
p
become equal to -1. It can be easily shown that such a
form
situation can take place only if α > 4. The homoclinic
α
y =x , (5) orbit forms at the value of y where the unstable branch
p
− 1−xp Spu(y) crosses the line x = −1, see Fig.4b. When the
map is firing spikes and the value of y gets to the bi-
where x 0. Equation (5) defines the branches of slow
p ≤ furcation point, the cycle Pk merges into the homoclinic
motion in the two-dimensional phase space (x ,y ) (see
n n orbit, disappears, and then the trajectory of the map
Fig 4). The stable branch S (y) exists for x <1 √α
and the unstable branch S p(sy) exists withinp1 √−α jumps to the stable fixed point xps.
pu − ≤ Typical phase portraits of the model, obtained un-
x 0
p ≤ der the assumptions made above, are presented in Fig.4.
3
Fig.4a shows the typical behavior for 2 < α 4. Here, getheranddisappear. Beforethis bifurcationthe system
≤
only two regimes are generated. The first regime is the is in x and, therefore, y increases. The termination of
ps
stateofsilence,whentheoperatingpoint(OP),givenby the burst is due to the bifurcation of the fast map as-
the intersection of x = 1+σ and one of the branches, sociated with the formation of the homoclinic orbit h .
s pu
−
isonthebranchS (y). Thesecondregimeistheregime Here, the limit cycle of the spiking mode merges into
ps
of tonic spiking, when the operating point is selected on the homoclinic orbit and disappears. After that the fast
the spiking branch S . subsystem flips to the stable fixed point x . Then the
spikes ps
process repeats (see arrows in Fig 4b).
It is clear from Fig.4a that when operating point is
5 4
(a) 6 set on Sps(y) the model will be in the regime of silence.
0
When the operating point is set on the branch S
spikes
S
spikes the model produces tonic spiking, unless the point is set
Spu(y) close to the formation of homoclinic orbit hpu. One can
see that, at the vicinity of this bifurcation, the branch
x −1
OP Sspikes becomes densely folded. As a result, the behav-
ior of y, which is governed by the mean value of x ,
n
x
s can become extremely sensitive to small perturbations
and even lead to instability caused by high-gain feed-
−2
S (y) back. This is oneofthe reasonsfor the irregular,chaotic
ps
spiking-burstingbehaviorwhich occursinthe mapwhen
−3.5 −3 −2.5 the operating point is set close to the area of the branch
y
S wherethis branchis densely folded. The detailed
spikes
andrigorousanalysisofchaoticdynamicscannotbedone
(b) Sspikes 8 7 6 withintheapproximationsmadeaboveandrequiremore
0 precisecomputationofS whichisbeyondthe scope
spikes
of this paper.
S (y) 8
pu
x
−1 bursts of
OP xs spikes
L
6 σ ts
−2 Sps(y) th
α
−4 −3.95 −3.9 −3.85
y
4 spikes
FIG. 4. Stable (Sps(y), Sspikes) and unstable (Spu(y))
branchesofslowdynamicsof(3)plottedontheplaneofphase silence
variables (y,x). The cases of α=4 and α=6 are presented
in (a) and (b), respectively. Numbers on the branch Sspikes
standforthevalueofkofthecyclesPk. Theoperatingpoints 2−1 −0.5 0 0.5 1
OP are selected to illustrate the regime of silence in (a) and σ
the regime of spiking-bursting oscillations in (b). Arrows in-
dicatethedirectionofslowevolutionalongthebranches,and
switching in thecase of (b).
FIG.5. Bifurcation diagram on theparameter plane (σ,α).
When α > 4 the picture changes qualitatively (see
The results of the analysis presented above are sum-
Fig 4b). Now, stable branches S (y) and S are
ps spikes marized in the sketch of the bifurcation diagram plotted
separated by the unstable branch S (y). If the operat-
pu on the parameter plain (σ,α) (see Fig 5). The bifurca-
ing pointis selectedonS (y), then the phase ofsilence,
pu tioncurveσ correspondstotheexcitationthreshold(6)
th
corresponding to slow motion along S (y), and spiking,
ps where the fixed point of the 2-d map becomes unstable
when system moves along S , alternate forming the
spikes and the map starts generating spikes. Curve L shows
ts
regime of spiking-bursting oscillations. The beginning of
the approximate location of the border between spiking
a burst of spikes corresponds to the bifurcation state of
and spiking-bursting regimes, obtained in the numerical
the fast map where fixed points x and x merge to-
ps pu simulations. Note that separation of spiking and burst-
4
ingregimesisnotalwaysobvious,especiallyintheregime regimeoftonicspikingwithα=5.0,σ =0.33,β =0. To
of chaotic spiking. The regime of spiking-bursting oscil- studytheresponsebehavior,positiveandnegativepulses
lations takes place within the upper triangle formed by of amplitude 0.8 and duration of 100 iterations were ap-
curves σ and L . The regimes of chaotic spiking or plied to the continuously spiking map. In the simula-
th ts
spiking-burstingbehaviorarefoundintherelativelynar- tions presented below the coefficient σe was selected to
row region of the parameters located around L . This be equal to one.
ts
regioncontainscomplexstructureofbifurcations,associ-
atedwithmulti-stableregimes,andisnotshowninFig5. 1
The bifurcation diagram shows the role of control pa-
n 0
rameters in the selection of dynamical features of the I
considered neuron model. Both parameters can be used −1
2
tomimicaparticulartypeofneuralbehavior. Parameter
1
σ can be used to model the external DC current injec- n
x 0
tion that depolarizes or hyperpolarizes the neuron. In
−1
this case σ can be written as
−3.4
σ =σ +I , (8)
u DC n
y
where σ is the parameter that selects the dynamics of −3.5
u
isolated neuron, and IDC is the parameter that models 0 500 1000
n
the DCcurrentinjectedintothe cell. Thechangesinbe-
haviorcausedbyI aresimilarto the actionofthe pa-
DC
rameter I in the well known Hindmarsh-Rose model [7]. FIG.6. Response of the model with σe =1 and βe =0 to
If the individual dynamics of the modelled neuron are apositivepulseofIn. Theparametersofthemapareselected
in the regime of tonic spiking (α=5.0, σ=0.33, β=0).
capable of generating only regimes of silence and tonic
spiking,anddonotsupporttheregimeofburstsofspikes,
then the value of α should be set below 4. In this case,
1
no matter what the level of I injected is, the system
DC
will not show bursts of spikes (see Fig 5). n 0
I
−1
2
III. MODELING OF RESPONSE TO THE
1
INJECTION OF CURRENT xn 0
−1
To study the dynamical regimes of neuronal behavior −3.4
in experiments, biologists change the type of neural ac-
n
tivityusingtheinjectionofelectricalcurrentintothecell y
though an electrode. It was shown above that the injec- −3.5
tion of DC current can be modeled in the map (1) using 0 500 1000
n
parameter σ = σ (see equation (8)). In this case, since
n
the external influence does not vary in time, the role of
parameter β = β is not important, because the behav- FIG.7. Response of the model with σe = 1.0 and βe = 0
n
ior of the map after the transient is independent of the to a negative pulse of In. The parameters of the map are
selected in the regime of tonic spiking (α = 5.0, σ = 0.33,
value of β. However, when the injected current changes
β=0).
in time, it may be useful to consider the dynamics of
parameter β in order to provide more realistic modeled
n
Figures 6 and 7 show the response to the positive and
behavior during the transient. Taking this into account,
negative pulse, respectively, computed with βe = 0. In
the input of the model can be considered in the form:
this case, during the action of a positive pulse, the value
βn =βeIn , σn =σeIn, (9) of yn increases monotonically because of the increased
value of σ (see (1b)). The increase of y pushes the
n n
where I is injected current, and coefficients βe and σe fast map up, see Fig 1. This leads to an increase in
n
are selected to achieve the desired properties ofresponse the frequency of spiking. After the action of the pulse
behavior. ends the value ofy monotonicallydecreasesback to the
n
This Section briefly illustrates how the relation be- original state (see Fig. 6).
tween these coefficients effects the dynamics of response Whenanegativepulseisapplied,itpushestheoperat-
tothepulseofI . Tobespecific,considerthemapinthe ing point down and, if the amplitude of the pulse is suf-
n
5
ficiently large,shutsoffthe regimeofspiking(see Fig.7). thatthe responsedynamicstothe positivepulse changes
Thishappensbecausethetrajectoryofthefastmapgets qualitatively. Indeed, when the pulse current is applied,
to the perturbed operating point which is now on the then, acting through β , it immediately forces the fast
n
stable branch S (y), (see Fig 4b). After the pulse is map to shift up. As a result, the rate of spiking and
ps
overthe spikingdoesnotre-appearimmediately because meanvalueofx increasessharply. Variabley reactsto
n n
the system spends some time drifting along the stable thischangeanddecreasesitsvaluetocompensateforthe
branchofslowmotionsS (y),duringwhichthevariable suddenchangeofthemeanvalueofx . Whentheaction
ps n
y overshootsits originalvalue for the spiking regime. As ofthepulse isover,thevalueofβ returnstoits original
n
a result, after the system switches to the spiking branch value and, due to the updated levels of y , the fast map
n
S ,y monotonicallydriftsdowntotheunperturbed overshoots its original state. As a result, the trajectory
spikes n
operating point. The dynamics of the slow evolutionare of the fast map reaches the stable fixed point. To return
clearly seen in the lower panel of Fig.7. to the originalregime of spiking the system has to go all
the way along the branches of slow motion S (y) and
ps
1 S backtotheoriginaloperatingpoint(seeFig.4b).
spikes
This type of response is observed in real neurobiological
n 0
I
experiments, see for example [15].
−1
Using the same analysis one can understand the new
2
effects in the response to a negative pulse caused by the
1
n action of β . These new effects can be clearly seen com-
x 0 n
−1 paring Fig.9 with Fig. 7.
Thepresentedresultsillustratehow2-dmap(1)along
−3.4
with equation (9) can model a large variety of transient
n
y neural behavior induced by injected current. Due to the
−3.5
simplicity of the model one can clearly see the nonlin-
ear mechanisms behind the response behavior andapply
0 500 1000
n them to select the desired balance between σe and βe to
model a particular type of response.
FIG.8. Responseofthemodelwithσe=1.0andβe =1.0
to a positive pulse of In. The parameters of the map are IV. REGIMES OF SYNCHRONIZATION IN TWO
selected in the regime of tonic spiking (α = 5.0, σ = 0.33,
COUPLED MAPS
β =0).
This section presents the results of studies of syn-
1 chronization regimes in coupled chaotically bursting 2-d
maps. The goal of this study is to reproduce the main
n 0
I regimes of synchronous behavior found in a real neuro-
−1 biological experiment [16]. The experiment was carried
1 outontwo electricallycoupledneurons (the pyloricdila-
n tors, PD) from the pyloric CPG of the lobster stom-
x −1
atogastric ganglion [17]. The regimes found in the ex-
−3 periment were also reproduced in numerical simulations
using Hindmarsh-Rose model [18,19], one-dimensional
−3.4
map model [14], and in experiments with electronic neu-
n
y rons [19].
−3.5
Theequationsusedinthenumericalsimulationsofthe
0 500 1000
n coupled maps are of form
x =f(x ,y +β ),
i,n+1 i,n i,n i,n
FIG.9. Responseofthemodelwithσe=1.0andβe =1.0
y =y µ(x +1)+µσ +µσ , (10)
i,n+1 i,n i,n i i,n
to a negative pulse of In. The parameters of the map are −
selected in the regime of tonic spiking (α = 5.0, σ = 0.33, where index i specifies the cell, and σ is the parameter
i
β =0). thatdefinesthedynamicsoftheuncoupledcell. Thecou-
plingbetweenthecellsisprovidedbythecurrentflowing
Figures 8 and 9 illustrate how response to the pulse from one cell to the other. This coupling is modeled by
of I changes when coefficient β is not equal to zero.
n n
To be brief and specific consider the case of βe = 1.0 βi,n =gjiβe(xj,n xi,n)
−
and σe = 1.0. Comparing Fig.8 with Fig.6 one can see σ =g σe(x x ) (11)
i,n ji j,n i,n
−
6
where i = j, and g is the parameter characterizing the This regime of synchronization of chaotic bursts is also
ji
6
strength of the coupling. The coefficients βe, σe set the observed in the experimental study of coupled PD neu-
balance between the couplings for the fast and slow pro- rons(seeFigure2cin [16]). Itisimportanttoemphasize
cesses in the cells, respectively. In the numerical simu- that, both in this simulation and in the experiment, the
lations the values of the coefficients are set to be equal: regimeofantiphasesynchronizationcharacterizedbythe
βe = 1.0 and σe = 1.0. The other parameters of the onset of regular bursts.
coupled maps (10) that remain unchanged in the simu- Thesimplicityofthismodelenablesonetounderstand
lations have the following values: µ = 0.001, α = 4.9, thepossiblecauseforthe onsetofregularburstinginthe
1
α = 5.0. The coupling between the maps is symmetri- antiphase synchronization. It is shown in Section II that
2
cal, g =g =g. chaotic dynamics of bursts occurs when the operating
ji ij
point of the system appears close to the leftmost area
of the spiking branch S . In this region, S is
2 (a) spikes spikes
denselyfoldedand,ifsystemslowsdowninthisarea,the
1
n timing for the end ofthe burstbecomes very sensitiveto
xi, 0 infinitesimal perturbations.
−1
2 (b) 0
(a)
1
n
xi, 0
−1
−0.5
2 (c) xf
1
n
xi, 0 −1
S (y)
−1 ps
−2
0 500 1000 1500
n −1.5
−3.7 −3.69 −3.68 −3.67 −3.66 −3.65
y
f
FIG. 10. Waveforms generated by the two coupled 2-d
maps (10) with σ1 = 0.240 and σ2 = 0.245. Waveforms of 1.0 (b)
x1,n and x2,n are shown by solid and dashed lines, respec-
tively. Three different regimes of spiking-bursting behavior
are shown. Individual chaotic spiking-bursting oscillations
0.0
of uncoupled maps, g = 0.0, - (a). Regime of synchronized
chaotic bursts, computed with g = 0.043, - (b). Regime of xf
antiphase synchronization, computed with g=−0.029, -(c).
−1.0
First,considerthemainregimesofsynchronizationbe- Sps(y)
tween the maps generating irregular, chaotic bursts of
spikes. To set the individual dynamics of the maps to
−2.0
this regime, the parameters σi, that take into account −3.7 −3.68 −3.66 −3.64
y
theDCbiascurrentinjectedintotheneurons(see[16]for
f
details), are tuned to the following values: σ = 0.240,
1
σ2 = 0.245. When these maps are uncoupled (g = 0) FIG. 11. Attractors computed from the waveforms of the
they produce chaotic spiking-bursting oscillations shown firstcelloperatingintheregimeofuncoupledoscillations-(a)
inFig10a. Whenthe couplingbecomessufficientlylarge and in the regime of antiphase synchronization - (b). Wave-
the slow components of the bursts synchronize while forms for these regimes are shown in Fig. 10a and Fig. 10c,
spikes within the bursts remain asynchronous see the respectively. The points of the attractors are connected by
waveforms in Fig. 10b obtained with g = 0.043. This lines to clearly show theirsequence in time.
regimeofsynchronizationistypicalfornaturallycoupled
PD neurons (see Figure 2a in [16]). This fact is illustrated in Fig. 11a. This figure shows
Introduction of negative coupling, g <0, also leads to theattractorcomputedfromthewaveformsx1,nandy1,n
synchronization of bursts, but in this regime of synchro- ofthe firstcellwhenthe cells areuncoupled(this regime
nization the systems burst in antiphase. Typical wave- is shown in Fig. 10a). To plot the attractor, the wave-
forms produced in this regime are presented in Fig. 10c. forms of x1,n and y1,n were filtered by a fourth order
7
low-pass filter with cutoff frequency 0.6. As a result, throughthe areaofcomplex behaviortowardthe branch
the attractor does not contain sharp spikes and the ap- S . Due to this fast change the duration of bursting
pu
proximate location S , which is close to the center is determined only by the dynamics of the silent phase,
spikes
of filtered spiking oscillations, is easy to see. The level which is regular.
of coordinate x , corresponding to the operating point, It was observed in the experiment that, in the regime
f
is equal to x = 1+σ = 0.76. One can see from oftonicspiking,spikesinPDneuronssynchronizeatlow
s 1
− −
Fig. 11a that when the center of oscillations gets close levels of coupling (i.e. without added artificial coupling,
to that level the trajectory becomes rather complex, as see [16] for details). This property of synchronization in
shown by the dense set of trajectories in the left part of the regimeof tonic spiking is also typicalfor the dynam-
the attractor. ics of the maps considered here (see Fig.12). In these
Thenumericalsimulationsshowthatthiseffectofreg- numerical simulations, the uncoupled maps were tuned
ularizationtakesplaceevenwhenσe =0. Therefore,the to generate spikes with slightly different rates. One can
vertical shift of operating point, x , caused by the slow clearly see the beats between the waveforms plotted in
s
couplingcurrentisnotveryimportantforthiseffect. The Fig.12a. When small coupling, g = 0.008, is introduced
coupling term β directly influences the fast dynamics the beats disappear as soon as the spikes get synchro-
i,n
by shifting the graph of the 1-d map (see Fig. 1) up or nized (see Fig. 12b).
down depending on the sign of β . In the regime of Although the synchronization between continuously
i,n
antiphase synchronization β is positive, when the i-th spiking maps is easily achieved, it is important to un-
i,n
cell is spiking, and negative,when it is silent. Therefore, derstand that the dynamics of such synchronization in
from the viewpoint of fast dynamics, the coupling forces the maps is a much more complex process than it is in
the cells to stay on the current branch of slow motion. the models with continuous time. Due to the discrete
One can understand these dynamics from the graph of time,theperiodicspikescanlockwithdifferentfrequency
f(x,y)and(1a). Thiseffectisalsoseenastheformation ratios. This results in the existence of a complex multi-
of extended shape of attractor plotted for the regime of stable structure of synchronization regimes. This struc-
antiphasesynchronization(seeFig.11b). Notethataxes ture becomes more noticeable in the dynamics of syn-
in Fig. 11a and Fig. 11b have different scales. chronizationas the number of iterations in the period of
spikes decreases.
2 (a)
1 V. DISCUSSION AND OUTLOOK
n
i, 0
x
−1
The simple phenomenological model of complex dy-
namics of spiking-bursting neural activity is proposed.
2 (b)
The model is given by two-dimensional map (1). The
1
n first equation (1a) of the map describes the fast dynam-
xi, 0 ics. Its isolated dynamics is capable of generating sta-
−1 ble limit cycles, which mimics the spiking activity of a
neuron, and a stable fixed point, which corresponds to
0 100 200 300 400 500
n phase of silence. This 1-d map has a region of param-
eters were both these stable regimes coexist. Existence
of such multistable regimes is the reason for generation
FIG. 12. Tonic spiking waveforms generated with of bursts, when the operating point, defined by the dy-
σ1 = 0.653 and σ2 = 0.714. Waveforms of x1,n and x2,n namics of the second equation (1b), is properly set (see
are shown by solid and dashed lines, respectively. Spiking in Fig 4b).
the uncoupled maps, g = 0.0, - (a). Regime of synchronized The shape of the nonlinear function f(x,y) used in
spikes, computed with g=0.008, - (b)
(1a) is selected in the form (2) based on the following
considerations: (i) The fast map should generate limit
Intheregimeofchaoticburststhelevelofx iscloseto
s cycles whose waveforms mimics those of the spikes. (ii)
the stable branch of spiking S , and, therefore, the
spikes Each spike generated by the map always has a single
slow evolution along the branch of silence S is faster
pu iteration residing on the most right interval x α +
thenalongS . Thismeansthatdurationofthephase ≥
spikes y. Therefore, the moment of time corresponding to the
ofsilence inthe cellis shorterthandurationofthe burst
appearanceofatrajectoryonthisintervalcanbeusedto
of spikes. When the cell switches from silence to spiking
define the time of a spike. This feature is important for
it sharply changes the value and the sign of β . This
i,n modeling the dynamics of chemical synapses in a group
changepushesS ofthespikingcelldowntoitsorigi-
spikes of coupled neurons. (iii) The analytic expression of the
nallocation,and,asaresult,quicklydragsthetrajectory
nonlinear function in interval x < 0 should be simple
8
enough to allow rigorousanalysis of bifurcations of fixed VI. ACKNOWLEDGMENT
points of the fast map. (iv) Use of the fixed level, 1,
−
in the rightmost interval of the function simplifies the The author is grateful to M.I. Rabinovich, R. El-
analysisofdynamicsattheendofthebursts. Theendof son, A. Selverston, H.D.I. Abarbanel, P. Abbott, V.S.
a burst is associated with the formation of a homoclinic Afraimovich and A.R. Volkovskii for helpful discussions.
orbit,which correspondsto the case whenthe trajectory This work was supported in part by U.S. Department
fromthisintervalmapsintothe unstablefixedpoint(see of Energy (grant DE-FG03-95ER14516),the U.S. Army
Section II). Research Office (MURI grant DAAG55-98-1-0269), and
Itisclearthattheshapeoffunction(2)canbemodified by a grant from the University of California Institute
to take into account other dynamical features that need for Mexico and the United States (UC MEXUS) and
to be modeled. Due to the low-dimensionality of the the Consejo Nacionalde Ciencia y Tecnologiade M´exico
model, the dynamical mechanisms behind its behavior (CONACYT).
are easy to understand using phase plane analysis. This
allows one to see how introduced modifications influence
the dynamics.
Some modification of the fast map is required to en-
hancethe regionofparameterswherethe modelcanstill
be used to mimic neural dynamics properties. For ex-
ample, from (2) and Fig. 1 one can see that, due to the [1] A.L. Hodgkin and A.F. Huxley, J. Physiol. 117, 500
dynamics of coupling terms, the trajectory of the fast (1952).
[2] T.R. Chay,Physica D 16, 233 (1985).
system can stay in the middle interval of f(x,y) for sev-
[3] T.R. Chay, Biol. Cybernetics 63, 15 (1990); T.R. Chay,
eral iterations. This can happen when the external in-
Y.S.Fan,andY.S.Lee,Int.J. Bifurcation andChaos 5,
fluence monotonically pushes the function up while the
595 (1995)
trajectory is located in the middle interval. In this case
[4] F.Buchholtz,J.Golowash, I.R.Epstain,andE.Marder,
the trajectory will map to the middle interval again and J. Neurophysiol. 67, 332 (1992).
again, increasing the duration of a spike. This artifact [5] D. Golomb, J. Guckenheimer,and S. Gueron, Biological
canbe removedusing one ofthe following modifications. Cybernetics 69, 129 (1993).
Introduction of a sufficient gap between the right end of [6] H.D.I. Abarbanel et al., USP. FIZ. NAUK. 166, 363
this interval and the diagonal (see Fig. 1) will help the (1996).
map to terminate the spike despite the monotonic ele- [7] J.L.Hindmarsh and R.M.Rose, Proc.R.Soc. London B
vation of the function. Alternatively one can introduce 221, 87 (1984); X.-J. Wang,Physica D 62,263 (1993).
[8] J.Rinzel,inOrdinaryandPartialDifferentialEquations,
an additional condition to the fast map (1a) that forces
edited by B.D. Sleeman and R.J. Jarvis, Lecture Notes
the map always to iterate its trajectory from the mid-
in Mathematics Vol. 1151 (Springer, New York, 1985),
dle interval to the rightmost one, despite the dynamics
pp. 304-316.
of y (see (2)). This can be achieved using the function
[9] J. Rinzel, in Mathematical Topics in Population Biol-
f(xn,y) of the following form ogy,Morphogenesis,andNeurosciences,editedbyE.Ter-
amotoandM.Yamaguti,LectureNotesinBiomathemat-
α/(1 x )+y, x 0
n n ics Vol. 71 (Springer, New York,1987), pp.267-281.
− ≤
f(xn,y)=α+y, 0<xn <α+y [10] V.N. Belykh et al., Eur.Phys. J. E 3, 205 (2000).
1, xn α+y or xn−1 >0 [11] E.M. Izhikevich, Int. J. Bifurcation and Chaos 10, 1171
− ≥ (2000)
An important feature of the model discussed here is [12] N.F. Rulkov,Phys, Rev,Lett 86, 183 (2001)
that one can use two inputs, βn and σn, to achieve the [13] G. deVries, Phys. Rev.E 64 051914 (2001).
desired response dynamics. Although these inputs are [14] B.Cazelles,M.Courbage,andM.Rabinovich,Europhys.
notdirectlyrelatedtodynamicsofspecificioniccurrents, Lett. 56 (4), 504 (2001).
they can be used to capture the collective dynamics of [15] D. Kleinfeld, F. Raccuia-Behling,and H.J. Chiel, Bio-
phys.J. 57, 697 (1990)
these currents. Selecting a proper balancebetweenthese
[16] R.C. Elson et al., Phys.Rev.Lett. 81 5692 (1998).
inputs, one can model a large variety of the responses
[17] R.M. Harris-Warrick et al., in Dynamics Biological Net-
that are seen in different neurons. Again, the simplic-
works: The Stomatogastric Nervous System, edited by
ity of the model helps one to understand the dynamical
R.M. Harris-Warrick et al. (MIT Press, Cambridge,
properties of each input and set the proper balance be-
MA,1992).
tween them. [18] H.D.I.Abarbanelet al.,Neural.Comput.81567(1996).
[19] R.D. Pinto et al.,Phys. Rev.E 62 2644 (2000).
9
