Quantitative modeling and analysis of bifurcation-induced bursting

1.1 MB·

by  J. E. Rubin

Quantitative modeling and analysis of bifurcation-induced bursting
J.E. Rubin
Department of Mathematics, University of Pittsburgh,
301 Thackeray Hall, Pittsburgh, PA 15260, USA
B. Krauskopf and H.M. Osinga
Department of Mathematics, University of Auckland,
Private Bag 92019 Auckland 1142, New Zealand
(Dated: January 23, 2017)
Modelingandparameterestimationforneuronaldynamicsareoftenchallengingbecausemanypa-
7 rameterscanrangeoverordersofmagnitudeandaredifficulttomeasureexperimentally. Moreover,
1 selecting a suitable model complexity requires a sufficient understanding of the model’s potential
0 use, such as highlighting essential mechanisms underlying qualitative behavior or precisely quan-
2 tifying realistic dynamics. We present a novel approach that can guide model development and
n tuning to achieve desired qualitative and quantitative solution properties. Our approach relies on
a thepresenceofdisparatetimescalesandemploystechniquesofseparatingthedynamicsoffastand
J slow variables, which are well known in the analysis of qualitative solution features. We build on
0 these methods to show how it is also possible to obtain quantitative solution features by impos-
2 ing designed dynamics for the slow variables in the form of specified two-dimensional paths in a
bifurcation-parameter landscape.

]
S PACSnumbers: 02.30.Oz,05.10.-a,87.10.Ed,87.19.ll
D Keywords: multipletimescales,fast-slowdecomposition,neuronal bursting,bifurcation
.
h
t I. INTRODUCTION fastsubsystemthatarevisitedoncesuchbifurcationsare
a
m crossed[4, 10,12]. Moreprecisely,bursting patterns will
be qualitatively the same if the associated paths traced
[ Pulsinginlasers[24],mixed-modeoscillationsinachemi-
by the slow variables cross the same bifurcations of the
calreaction[15,22],andburstinginneurons[8]arejusta
1 fastsubsysteminthesameorder. However,thisapproach
v fewexamplesamongtheubiquitouspatternsthatexhibit does not reveal quantitative features of the oscillation,
5 features evolving on disparate time scales. One classical whichmaybequiteimportantforthe particularapplica-
8 approach to understanding such patterns is the method
tion associated with the pattern.
6 of fast-slow decomposition of the underlying mathemat-
In this paper, we propose a novel extension to the ap-
5 ical model, where each variable in the model is classi-
0 plication of fast-slow decomposition, with two aims in
fied as fast or slow. All slow variables are then fixed, so
. mind. First, we extend the method to the quantitative
1 that the fast variables define the so-called fast subsys-
realm, so that it can be used to guide the quantitative
0 tem for which the slow variables are treated as parame-
fitting of data arising from experiments. Second, we il-
7
ters [14, 17–19]. The slow dynamics impose a particular
1 lustrate that fast-slow decomposition can be used as a
paththat the slowvariablessweepout inthe bifurcation
: modeling tool for systems where not enough details are
v landscape of the fast subsystem. The full system dy-
known to model the dynamics of the slow variables ac-
i namics can now be approximated by assuming that the
X curately. The underlying idea of our approach is to pre-
fast variables drift along a sequence of attractors, punc-
r scribe a simple form of slow dynamics, leading to a par-
a tuated by occasional rapid transitions between different
ticular class of imposed paths for which the associated
attractors at certain bifurcation events, while the slow
parameters can be modulated. Importantly, we consider
variables trace the path imposed by the slow dynamics.
two slow variables, which provides access to a param-
Inparticularwhenasingleslowvariableisused,theclas-
eterized continuum of paths that can be used to tune
sicalmodelingapproachconsidersthelocusofzerospeed,
quantitative features of anoscillationpattern. The com-
or nullcline of the slow variable with respect to the fast
bination of imposed paths and easily obtained quantita-
variables, and includes an additional ordinary differen-
tive information about the fast subsystem can then be
tial equation accordingly [4, 6, 12]; another approach is
used to steer the fast-variable outputs to a quantitative
toviewtheslowdynamicsasanon-automomousexternal
agreementwithdesiredtargets. Theresultingsystemcan
force [10, 11, 16].
serveasamodelofthephenomenonofinterestinitsown
The method of fast-slow decomposition offers a sim- right, but also as a guide for subsequent parameter esti-
ple and effective way of classifying different oscillatory mationinthefullmodel,ifoneisavailable,andasatool
patterns for two-time-scale systems. Each type of pat- for the development of a model for the slow components
tern is uniquely described by the sequence of bifurca- of the system, if one is not available.
tions that are encountered as the slow variables trace Inwhatfollows,weconsiderthespecific fieldofneuro-
a selected path, and by the particular attractors of the physiology and focus on patterns generated by neurons.
2
Neuronsengageinactivitypatternsknownasburstingin The quantitative features of DB bursting may be im-
a variety of settings, including sleep, novelty detection, portantforthebiologicalfunctionofaneuron;forexam-
generationofrepetitive movements,releaseofhormones, ple, the spiking phase within a burst may be associated
and certain pathological conditions [5, 8]. In very gen- with release of substances such as hormones or with ac-
eralterms,burstingreferstoanytimecourseofthemem- tivation of a particular muscle group. Moreover, it may
branepotentialthatfeaturesactivephasesofconsecutive not be obvious that two bursting solutions, which are
high-frequency oscillations alternating with intervals in qualitatively similar from the point of view of fast-slow
which oscillations are much smaller, much more infre- decomposition, really share what biologists would con-
quent, or absent altogether. This simple characteriza- sidertobethesamefeatures. Forexample,ifoscillations
tion, however, encompasses a striking diversity of burst- become sufficiently small, they will be undetectable or
ingpatternsthatariseacrossdifferentneuronsinvarious swampedbynoise,whileonlycertainfrequenciesofoscil-
contexts [4, 12, 16]. lations may suffice to achieve a biological purpose; thus,
Classic analysis of bursting dynamics, as pioneered by twosolutionswithoscillationsofdifferentamplitudesand
Rinzel [18, 19], utilizes the fast-slow nature of such sys- frequenciesmaymeritadistinctionthatisnotpresentin
tems and is based on singularity theory [4, 10, 16]. The qualitative analysis of fast-slow decomposition.
keyideais toinvestigateburstingpatternsandto design
minimalbursting models via the analysisofthe underly- What determines the quantitative features of a burst-
ing fast subsystem. In particular, systems with a single ing pattern in a fast-slow system? We present a case
slow variable havebeen exploredextensively in this way. studyillustratingthatanaturalextensionofthemethod
Forsuchsystems,animposedpathassociatedwithape- offast-slowdecompositioncanprovideusefulinformation
riodic bursting pattern is necessarily a line segment and for quantitatively fitting data on bursting arising from
the slow variable oscillates back and forth over a given biological experiments. Here, we go beyond the consid-
rangeof values. The bursting patterns obtainedfor such eration of which bifurcation curves are crossed and take
systemsareentirelycharacterizedbythenumberandor- into account details of the path in the plane traced by
derofbifurcationsofthefastsubsystemthatareencoun- the twoslowvariablesinbetweenthese crossings,aswell
tered along the imposed path; for example, a complete as corresponding details about the attractors of the fast
classification of bursting patterns was attempted in [12]. subsystem encountered along this path.
There are cases of bursting that involve two slow vari-
This paper is organized as follows. In the next sec-
ables,however. Whatisknownasparabolicburstingisa
tion, we presentthe relevantanalysisofthe DB bursting
notable case, for which the qualitative pattern of burst-
modelfrom[21]andcompareitsdynamicswiththatofa
ing has been explained in terms of slow motion across
four-dimensional reduced model. The full details of the
bifurcationcurvesofthefastsubsystem[20];moreover,a
seven-dimensional model can be found in [21], but the
canonicalunderlying model associated with this bifurca-
system of equations is also given in the Appendix. Sec-
tion structure has been derived via nonlinear coordinate
tion III illustrates how imposed paths for the two slow
changes [7, 23]. Our addition to this set of ideas is that
variables in the four-dimensional model can be used to
quantitative features can be included in the designstage
obtain quantitatively similar bursting patterns as exhib-
if two (or more) slow variables are used — by consider-
ited by the full seven-dimensional model. This idea is
ing anentire family of paths thatencounter the required
then utilized in Sec. IV to explore how other aspects of
sequence of bifurcations.
the fast-subsystem dynamics, which are related to but
As the leading example for this paper, we consider a
go beyond the analysis of its bifurcation diagram, deter-
seven-dimensional model of bursting in neurons of the
mine quantitative features of the bursting patterns. We
respiratorybrainstem,system[21]withfivefastandtwo
end with a discussion in Sec. V.
slow variables. This model exhibits a particular form of
neuronal bursting, which we refer to as depolarization
block or DB bursting; an example of the corresponding
voltage time course is shown in Fig. 1(b). This type
of bursting pattern has also been observed in other
neuronal models [3, 9]. Eachcycle within a DB bursting
pattern consists of:
II. THE DB BURSTING MODEL AND ITS
(i) a silent phase of quiescence; DYNAMICS
(ii) the emergence of sustained voltage spiking with
gradually increasing spike frequency;
(iii) a fairly abrupt significant increase in spike fre- Rubin et al. [21] presented a seven-dimensional model
quency and attenuation of spike amplitude; of bursting in neurons of the respiratory brain stem.
(iv) an approach towards a depolarization block state The model is based on Hodgkin–Huxley formalism and
of steady, elevated voltage; is given by the following system of ordinary differential
(v) a re-emergence of spiking; and equations involving five fast and two slow variables; see
(vi) a return to quiescence. the Appendix for a complete description of the various
3
Rubinet al. [21]usedfast-slowdecompositionto high-
(a)
light the underlying bifurcation set of (1) that dictates
5.6
the qualitative form of the DB bursting pattern. They
Na found that there are two types of bifurcations, namely,
a saddle-node bifurcation on an invariant cycle, de-
5.4 SNIC AH notedSNIC,andanAndronov–Hopfbifurcation,denoted
AH; Figure 1(a) shows the curves SNIC and AH in the
(Ca,Na)-plane. For values (Ca,Na) in between the
5.2 curves SNIC and AH, an attracting periodic orbit of (1)
exists togtherwith a singleequilibrium fromwhich itbi-
furcatesatAH;the equilibriumalsoexiststo the rightof
AH, where it is stable. For values (Ca,Na) to the left
5
0 0.5 1 Ca 1.5 of SNIC, three equilibria exist; one of these is stable and
corresponds to hyperpolarized voltages near the resting
(b)
potential. During one period of the DB bursting pat-
20
tern in the full system (1)&(2), the path taken by the
V] slow variables Ca and Na is such that bifurcations are
m 0 crossed four times; this is illustrated in Fig. 1(a), where
[
v the DB bursting pattern is shown in projection onto the
-20 (Ca,Na)-plane. Our hypothesis is that the quantitative
features of this particular path taken by the two slow
variables Ca and Na control the specific pattern of the
-40
time course of v shown in panel (b).
The stages of the DB bursting pattern listed in Sec. I
-60 can be explained in terms of the trajectory path in the
0 1000 2000 3000 4000 5000 6000 (Ca,Na)-plane. Startingfromapointontheperiodicor-
t [ms] bit locatedtop-left in the (Ca,Na)-plane, the trajectory
lies near the family of hyperpolarizedrest states and the
FIG.1. DBburstingintheseven-dimensionalsystem(1)&(2). DBburstingpatternisinthequiescentorsilentphase(i);
Panel(a)showstheburstingoscillationinprojectionontothe upon crossing the curve SNIC, sustained voltage spiking
(Ca,Na)-planeoverlaidonthebifurcationsetofthefastsub- emerges (ii); as soon as the curve AH is crossed, the fre-
systeminthisregion,givenbythelociofsaddle-nodebifurca- quency of the spiking increases significantly, while the
tions on an invariant cycle and Andronov–Hopf bifurcations,
amplitude decreases during stages (iii) and (iv); spiking
labelled SNIC and AH, respectively. Panel (b) shows thecor- re-emergesasthepathcrossesAHagain(v),whichisfol-
responding time course of voltage v. lowed by a return to silence after crossing SNIC for the
second time (vi). Hence, the active phase of each burst,
namely, the epochs of voltage spiking and elevated volt-
functions and parameter values used:
age in stages (ii)–(v), encompasses the segment on the
v′ = −1c [IL(v)+IK(v,n) side ofthe curveSNICthat containsAH;the attenuation
 +INa(v,m,h)+Isyn(v,s) of spike amplitude and approach towards depolarization
mn′′ == ττmn11((vv)) [[nm+∞∞I((CvvA))N−−(vnm,],C],a)+Ipump(Na)], (1) biavnallsogrosciakaonsfotsfAocscyHticaa;ltgeceeoadsmnc(dpiuiairA)rvenaesdnwrdoiotf(nhsioav[vd2)–d1oH]lc.eoc-BnuporafdrbbreieeftubtworifceuaaerntncdiaottnChieorwnetiswtoshomnbcaaurnnoriss[n3ts---]
 hs′′ == τh1(τv1s) [[(h1∞−(vs))−s∞h(],v)−ks], idanigmrFeoeddrnyusnctiheaodminsaimlcpsaom.pdoeedrl.,elwIanesatonhuirnackdtauotafaltsahepetp,oloiucrtapttauiortng,oetfitothuwetopusueltdvefboner-
′
Ca = ε [k s−k (Ca−Ca )],
IP3 Ca b (2) replaced by experimentally recorded data.
′
(cid:26)Na = α [−I (v,Ca)−I (Na)].
CAN pump
Here, system (1) comprises equations for the five fast
variables v, n, m, h, and s, and system (2) comprises III. IMPOSED PATHS FOR A
FOUR-DIMENSIONAL REDUCED MODEL
equations for the two slow variables Ca and Na. Hence,
the five-dimensional fast subsystem depends on two pa-
rameters that yield curves in the (Ca,Na)-parameter OursimulationsshowthatafairlysimilarDB-likeburst-
plane along which bifurcations occur. These curves can ing pattern can be obtained from a reduction of the
readilybecomputedwithstandardsoftwarepackages,for seven-dimensional model of Rubin et al. [21] to a four-
example, with XPPAUT [6]. dimensional version. This reduction is achieved by a
4
of traversals of bifurcation curves. Despite the overall
(a)
similarity between the systems, comparison with Fig. 1
Na
reveals some differences with the DB bursts from the
6.2 seven-dimensional model, which we found to persist de-
spiteextensiveexplorationofparameterspace: intheDB
SNIC bursts of the reduced model, stage (ii) of sustained volt-
5.8
age spikes is brief, the voltage spikes are less attenuated
intheapproachtowarddepolarizationblockinstage(iv),
5.4 andthere-emergenceofspikingisprolongedinstage(v).
These effects are linked in Fig. 2(a) with the behavior of
AH the slow variable Ca; the projection of the bursting so-
5
lutionofthefour-dimensionalsystemontothe(Ca,Na)-
0 0.1 0.2 Ca 0.3 planeshowsthatduringthespikingphase,thetrajectory
dancesbackandforthacrossthecurveAH,failingtocross
(b)
over it completely.
20 The comparison between Figs. 1 and 2 indicates that
V] the model reduction from seven to four dimensions in-
m 0
[ duces an unwanted effect on the behaviour of the slow
v
-20 variablesCaandNa. Ratherthantryingtofittheslow-
variable dynamics, however, our approach is to ignore
-40
their biological relevance completely. That is, we will
demonstrate that, given a particular bifurcation struc-
-60
ture of the fast subsystem, we can generate a desired
-80 quantiative pattern by imposing a path on the associ-
ated slow variables. For this purpose, unwanted effects
-100
0 500 1000 1500 2000 2500 3000 inaparticularparametertuningarenotrelevant,aslong
t [ms] as the bifurcation structure is present. To this end, we
replace the biological (Ca,Na)-dynamics with imposed
FIG. 2. DB bursting in the reduced four-dimensional sys- paths in the (Ca,Na)-plane, so that the combined sys-
tem (3)&(2). Panel (a) shows theburstingoscillation in pro- temhastherelevantqualitativeandquantitativefeatures
jectionontothe(Ca,Na)-planeoverlaidontheSNICandAH ofDBburstingasgivenbythe seven-dimensionalmodel;
curves. Panel (b) shows the corresponding time course of see also [10, 16] for similar ideas with one-dimensional
voltage v;compare with Fig. 1.
imposed paths.
More specifically, we consider a parameterized path P
inthe(Ca,Na)-planeintheformofanellipsewithprin-
quasi-steady-state approximation, made by setting m =
cipal axes along the Ca- and Na-axes. Such an ellipse
m∞(v) and s=s∞(v)/(s∞(v)+k), and by the classical
reasonablyresemblesthepathshowninFig.1(a),yetitis
stepofreplacinghwith 1−1.08n[13,18]. The fastsub-
definedby only five parametersthatstill provideenough
system is then reduced to the following two-dimensional
freedom to adjust quantitative features of the dynamics.
model:
A sixth parameter is the speed ε with which the ellipse
v′ = −1 [I (v)+I (v,n) is traced. Hence, P is defined as
c L K
 +INa(v,m∞(v),1−1.08n)
P = P(Ca ,Na ,d,Ca ,Na ,ε)
 ++IIsCyAnN(v(,vs,∞C(av))+/(sIp∞u(mvp)(+Nak))]), (3) Cac(t)=cCac +0 (C0a0−Cac) cos(εt)
 n′ = τn1(v) [n∞(v)−n]. :=Na(t)=Nac +− d((NNaa00−−NNaacc))csoins((εεtt)),
Figure 2 shows DB bursting for the reduced sys- + d1(Ca0−Cac) sin(εt).
tem (3)&(2). The DB bursting pattern of (3)&(2) is 
again projected onto the (Ca,Na)-plane in panel (a), Here,Cac andNac definethecenteroftheellipse,disits
along with the bifurcation curves SNIC and AH for the aspect ratio, and (Ca0,Na0) is a chosen initial point at
two-dimensionalfast subsystem (3), and the time course time t = 0. The ranges for Ca and Na are the intervals
for v is shown in panel (b). Due to the approximations [Cac−δ, Cac+δ]and[Nac−d1δ, Nac+1dδ],respectively,
used,anadjustmentofsomeoftheparametervalueswas where δ = (Ca −Ca )2+d2(Na −Na )2. For the
0 c 0 c
required for achieving the characteristics of DB burst- paths usedpin this paper, we fix Na0 = Nac and use
ing of the target system (1)&(2); see the Appendix for Ca to tune the path width and d to control its aspect
0
details. One key ingredient is that we have obtained a ratio. Weassume0<ε≪1,suchthatP istracedinthe
similarbifurcationstructureofthefastsubsystemfor(3) counter-clockwisedirectionandCaandNaevolveslowly
asfor(1). Anotherisachievingtheappropriatesequence in time.
5
It is convenient for our purposes that the evolution
(a)
along an ellipse P can also be considered as the solution
Na
of the system of ordinary differential equations SNIC
6
′
Ca = −ε d(Na−Na ),
c
(4)
(cid:26)Na′ = ε 1(Ca−Ca ),
d c
5
for the initial conditions Ca and Na . By combining
0 0
system (4) with the fast subsystem (3) of the reduced AH
four-dimensional model, we obtain the four-dimensional
driven model (3)&(4) that is the main subject of this 4
paper. Notethattheformulation(4)sufficesforourpur-
poseshere,eventhoughitsperiodicorbits—the ellipses 0 0.4 0.8 Ca 1.2
— are not isolated. In case continuation of the overall
periodicorbitofthedrivensystem (3)&(4)isrequired,a (b)
20
family ofisolatedattracting ellipticalperiodic orbits can
be obtained in similar fashion from a scaled version of V]
the Hopf normal form [6]. m 0
[
v
-20
IV. RESULTS
-40
TheapproachofimposingapathP forthetwoslowvari-
ables Ca and Na allows us to control and explore their
qualitative and quantitative effects systematically. Our -60
main result is that we are able to quantify different fea-
0 200 400 600 800 1000 1200 1400
turesoftheDBburstingpatternsintermsoftheunderly- t [ms]
ingfastsubsystem;thesefindingsgobeyondthestandard
approach of determining the bifurcation diagram of the
FIG. 3. DB bursting patterns for the seven-dimensional
fast subsystem. driven system (1)&(4) with Cac =0.7, Nac =5.35, Ca0 =0
andε=0.009, andthreedifferentvaluesford,namely,d= 1
2
(blue), d=2 (purple) and d=20 (red). Panel (a) shows the
A. Bursting with an imposed (Ca,Na)-ellipse three paths P in the (Ca,Na)-plane overlaid on the loci for
SNIC and AH, and panel (b) shows the three corresponding
time courses for v.
We expect that a driven system, where the two slow
variables Ca and Na evolve as determined by an im-
posed elliptic path, can exhibit DB bursting if the fast
subsystem gives rise to a region in the (Ca,Na)-plane andAH inthespecifiedsequenceneededforDB bursting
bounded by bifurcation curves SNIC and AH. We now and,asillustratedinpanel(b),allyieldDBburstingpat-
confirm that claim by showing that the qualitative fea- terns that exhibit the qualitative features described for
tures of DB bursting are obtained by choosing a path P stages(i)–(vi);notethatthepathwithd= 12 yieldssome
that crosses these curves in the specific sequence SNIC, tonicspikingattheonsetofeachactivephasethatisalso
AH, AH, and SNIC. seen in the DB bursting of system (1)&(2), and which is
We first explore this for a driven version of the seven- a quantitative feature not seen for the other paths.
dimensional model, that is, by combining the fast sub- Qualitatively similar DB bursting can be obtained for
system (1) with the driven system (4). Figure 3 shows the four-dimensional driven system (3)&(4). Since the
three different DB bursting patterns for system (1)&(4), loci of SNIC and AH for the fast subsystem (3) have
usingthreedifferentimposedpathsfortheslowvariables shifted slightly in the (Ca,Na)-plane, slightly different
Ca and Na. The imposed paths in the (Ca,Na)-plane pathsmustbeimposed. Figure4showstheeffectofthree
are shown in panel (a) overlaid on the loci for SNIC and differentellipsescenteredatCa =0.15andNa =5.85,
c c
AH of the fast subsystem (1); the corresponding three with Ca = 0 and ε = 0.004. Again, only the aspect
0
time coursesforv areshowninpanel(b). Foreachpath, ratio d was varied, where we chose d = 1 (blue), d = 1
5
we used Ca = 0.7, Na = 5.35, with initial condition (purple) and d = 50 (red). As for Fig. 3, the associated
c c
Ca = 0 (recall that we always set Na = Na ) and paths are overlaid in panel (a) on the loci for SNIC and
0 0 c
constant speed ε = 0.009. Only the aspect ratio d was AH in the (Ca,Na)-plane, and the three corresponding
varied,namely, d= 1 (blue) for the fattest ellipse shown time courses for v are shown in panel (b). Note that
2
in panel (a), which was increased to d = 2 (purple) and only a short segment of each of the three imposed ellip-
then d=20 (red). All three paths cross the curves SNIC tic paths lies in the region to the right of AH, which
6
ing the aspect ratio d of the elliptic path, as well as the
(a)
speed ε at which this path is traced. For example, com-
Na SNIC parison between Figs. 3 and 4 suggests that the speed
with which Ca and Na traverse the same path leads to
6 quantitative differences in the bursting. We next pursue
the idea oflinkingmodel parametersassociatedwith the
imposed path in (Ca,Na) space to quantitative features
of the DB pattern.
5 AH
B. Quantitative features of bursting beyond the
basic bifurcation structure
0 0.1 0.2 Ca 0.3
40 We now explore how the dynamics of the fast subsys-
(b)
tem can be analysed to determine quantitative features
20
of bursting solutions for the full system. Here, we fo-
V]
m 0 cus exclusively on the reduced four-dimensional driven
[ system (3)&(4).
v
-20 Wefirstreturntotheconsiderationofthetraversalrate
of the imposed path. Figure 5 shows time courses of v
-40
forthreequantitativelydifferentburstingpatternsforthe
-60 four-dimensionaldrivensystem(3)&(4)withε=0.002in
panel(a),ε=0.006inpanel(b),andε=0.01inpanel(c)
-80
for an imposed elliptic path with the parameters Ca =
c
-100 0.15, Nac =5.85, Ca0 =0 and d=0.1. Consistent with
0 500 1000 1500 2000 2500 3000 the previous subsection, when ε is small, each phase of
t [ms] the burst is elongated due to the smaller traversal rate;
note the different t-scales in the three panels. Moreover,
FIG.4. DBburstingpatternsforthefour-dimensionaldriven as ε increases, the faster traversal of the same path in
system (3)&(4) with Cac = 0.15, Nac = 5.85, Ca0 = 0 and the (Ca,Na)-plane yields less attenuation of the voltage
ε = 0.004, and three different values for d, namely, d = 1 oscillationsduringthedepolarizationblockphaseofeach
5
(blue), d=1 (purple) and d=50 (red). Panel (a) shows the burstcycle. As is particularlyclearinFig.5(a),a slower
three paths P in the (Ca,Na)-plane overlaid on the loci for traversalintensifies this attenuation.
SNIC and AH, and panel (b) shows the three corresponding
We dig deeper into the source of quantitative differ-
time courses for v.
ences in patterns by next considering the time intervals
between spikes in stage (ii) of a DB bursting pattern,
during which repetitive, high-amplitude voltage spikes
regulates stages (iv) and (v) of DB bursting. Conse- occur. These spikes arise from a family Γ of stable
0
quently, the approach towards and away from the depo- periodic orbits of the fast subsystem that exist in the
larization block state of steady, elevated voltage is rela- (Ca,Na)-region in between the bifurcation curves SNIC
tively brief; compare Figs. 3(b) and 4(b). Also, the path and AH. The time intervals between spikes can vary sig-
with d = 1 (blue) yields particularly deep hyperpolar- nificantly across successive spike pairs within a bursting
5
izations between bursts, and maximal spike attenuation solution, and the time until emergence of the first spike
during bursts, which is also the case for the path with after th esecond crossing of SNIC can differ appreciably
d = 1 (blue) in Fig. 3. The difference in the range for across parameter values and paths. Our hypothesis is
2
the t-axes in Figs.3(b) and 4(b) is entirely controlledby that these quantitative features observed in stage (ii)
the different choices for ε. of DB bursting are controlled by characteristics of the
Thus, as expected, the qualitative features of the (Ca,Na)-dependent family Γ .
0
DB bursting pattern for both the full seven-dimensional Specifically,wenotethateachperiodicorbitinΓ hasa
0
driven system (1)&(4) and the reduced four-dimensional well-definedperiodthatdependsonthechoiceforCaand
driven system (3)&(4) are entirely determined by the Na. The periods encountered along a path P are deter-
type of bifurcations exhibited by their respective fast mined by the locationof this imposedpath in the region
subsystems and the specific order in which these bifu- bounded by SNIC and AH. We find that the progression
ractionsareencountered. As illustratedinFigs.3 and4, of periods encountered alongP strongly shapes the time
we can easily control such properties by choosing Na intervals between spikes in stage (ii) of the resulting DB
c
largeenoughand deciding onappropriatevalues for Ca burstingpattern. Figure6(a)showsthecurvesSNICand
c
and Ca . We also observe, however, that quantitative AHinthe(Ca,Na)-planetogetherwithcontoursofequal
0
differences in the bursting patterns are obtained by tun- period for the periodic orbits of the family Γ . The pe-
0
7
(a) (a) 10.1
7
20
Na
V] AH
m 0 6 SNIC
[
v
-20
5
-40 15.0
12.0
4
-60
0 500 1000 1500 t [ms] 0 0.2 Ca 0.4
(b) 40
(b)
V] 0
20 m
V] v[-40
m 0
[
v 0 100 200 300 t [ms]
-20 (c) 40
-40 V] 0
m
[
v -40
-60
0 100 300 500 t [ms] 0 100 200 300 t [ms]
(c)
FIG. 6. Interspike intervals of the DB bursting pattern are
significantly affected by the different periods of the periodic
20
orbits in the family Γ0 of the fast subsystem that are en-
V] countered along the imposed path. Panel (a) shows SNIC
m 0 and AH in the (Ca,Na)-plane along with contours (grey) of
[
v equal period for periodic orbits in Γ0. Overlaid are a verti-
-20 cally wide (blue) path with Cac =0.15, Nac =5.2, d =0.1,
and Ca0 = 0 and a horizontally elongated (red) path with
-40 Cac = 0.1, Nac = 5.1, d = 1, and Ca0 = −0.1. The cor-
responding time courses of v over the first half period are
shown in panels (b) and (c), respectively; their crossings of
-60 the curves SNIC and AH are marked with green and black
0 100 200 300 t [ms] circles, respectively.
FIG. 5. The same path in the (Ca,Na)-plane can yield
quantitative differences in bursting features, depending on
Na = 5.1, d = 1, and Ca = −0.1. The corresponding
the speed at which it is traversed. Shown are the time c 0
time traces of v over the first half period are shown in
courses of v for the four-dimensional driven system (3)&(4)
panels(b)and(c)ofFig.6,respectively. Thefilledblack
resulting from the same imposed elliptic path with speeds
ε=0.002inpanel(a),ε=0.006inpanel(b),andε=0.01in andgreencirclesindicatethemomentwhenthesolutions
panel(c). TheotherparametersforthepathareCac =0.15, cross SNIC and AH, respectively.
Nac =5.85, Ca0 =0, and d=0.1. Thefat(blue)pathcrossesSNICroughlyorthogonalto
the period contours and thus progresses quickly to a re-
gionwith relativelylow periods. The interspikeintervals
riodsassociatedwiththe contoursprogressivelydecrease in the correspondingtime trace of v, shown in Fig. 6(b),
as we move farther from SNIC, starting from 40 for the are relatively short, with little change over the period of
curve closest to the SNIC, although close to AH some the burst. In contrast, the thin (red) path crosses SNIC
non-monotonicity sets in (e.g, note the turning point in a direction that is aligned with the high-period con-
for the labeled contour with period 10.1). This panel tours close to the curve SNIC. Correspondingly, there is
also includes two elliptic paths with ε = 0.009: the fat alongerdelayfromthe crossing(blackcircle)tothefirst
(blue) path is given by Ca = 0.15, Na = 5.2, d = 0.1, spike in the time trace of v, shown in Fig. 6(c), and the
c c
and Ca = 0, and the thin (red) path by Ca = 0.1, first few interspike intervals are longer than any of the
0 c
8
about the equilibria associated with this elevated state.
(a) 0.03
0.01 NearthecurveAH ofAndronov–Hopfbifurcations,these
7 SNIC AH equilibria all have complex-conjugate pairs of eigenval-
-0.01
Na ues. Their realparts are negative to the rightof AH and
-0.03
positivetotheleftofit. Figure7(a)shows15contoursin
-0.05
6 the (Ca,Na)-plane of these real parts Re(λ), uniformly
distributedfrom−0.05to0.09;thecurveAHcorresponds
tothecountourRe(λ)=0. Weselectthreeellipticpaths
5 that are all centered at Cac =0.19 and Nac =5.75 and
use Ca = 0.04 and ε = 0.004. The only difference in
0
0.09 thesecontoursisintheiraspectratios,whichared=0.1
0.07 (blue),d=0.2(purple)andd=0.4(red). Theirintersec-
0 0.2 Ca 0.4 tionswithAH aremarkedwithsquaresandtriangles,in-
dicating the first andsecond crossings,respectively. The
(b)
corresponding time courses for v are shown in Fig. 7(b),
0 where time is shifted so that the first crossing occurs at
V] t=0 (green square). Here, the burst with d=0.2 (pur-
m
[ ple, top) is plotted in actual coordinates, while the ones
v with d = 0.1 (blue, bottom) and d = 0.4 (red, middle)
-100 havebeenshifteddownby200and100mV,respectively.
As shown in Fig. 7(b), stage (iv) of the DB burst-
ing patterns is rather similar acrossthese parameterval-
ues, even though the times until the second crossing of
-200
AH (triangles) are different. Observe that the times be-
tween AH crossings for the bursts with d = 0.2 (purple,
-400 -200 0 200 400 600 800 1000 top) and d = 0.4 (red, middle) are similar, while the
t [ms]
burstwithd=0.1(blue,bottom)exhibitsasignificantly
longerstage(iv). Fig.7(a)illustratesthattheburstwith
FIG. 7. Contraction to and expansion away from the de- d=0.4 (red) does not cross the contour Re(λ)=−0.05,
polarization block state depends on the real part Re(λ) of which indicates that the net contractiontowards the de-
theeigenvaluesoftheassociated equilibriainthefastsubsys-
polarization block state is not as strong as for the other
tem that are encountered along the imposed path. Panel (a)
two bursts. Consequently, even though the bursts with
shows SNIC and AH in the (Ca,Na)-plane along with con-
d = 0.2 (purple) and d = 0.4 (red) spend about the
tours (grey) of Re(λ). Overlayed are paths with Cac =0.19, same time to the right of AH, the time needed to ex-
Nac = 5.75, Ca0 = 0.04 and ε = 0.004, and varying d =0.1
(blue), d = 0.2 (purple) and d = 0.4 (red). The correspond- pandawayfrom the depolarizationblock state is shorter
ingtimecoursesofvareshowninpanel(b)withtimeshifted for the burst with d = 0.4 (red), resulting in a shorter
so that the first crossing occurs at t=0 (green square), and time until the onsetof stage (v); see Fig. 7(b). Note fur-
voltageshifted downby200and100mV for thetimecourses therthatthe longertime spenttothe rightofAHforthe
withd=0.1(blue)andd=0.4(red),respectively;thesecond burstwithd=0.1(blue)resultsinastrongercontraction
crossing of AH is marked with triangles. towards the depolarization block state and subsequent
slowerexpansionawayfrom it. Hence, this case includes
virtually no additional spikes before the initiation of the
silent phase at the end of stage (v).
ones shown in Fig. 6(b), with a gradual compression of
In fact, the times spent between crossings of bifurca-
spiketimesoverthecourseoftheburst. Notethatwese-
tion curves and the rates of contraction and expansion
lected parameter values such that the overalltimes from
encounteredduring these times, as indicated by the con-
the start (black circles) to the end (green circles) of the
tours in Fig. 7(a), combine in such a way that the to-
active phase, in between SNIC and AH, are quite similar
tal number of additional spikes exhibited in stage (v),
for both DB bursts, so the longer interspike intervals in
before the return to the silent phase, depends non-
Fig. 6(c) are not simply due to a slower passagethrough
monotonicallyonthe aspect ratioof the imposed elliptic
stage (ii) of the DB bursting pattern.
path. ThetimecoursesinFig.7(b)areorderedtoempha-
Next,weturntostages(iv)and(v)oftheDBbursting size this non-monotonicity; they reveal that the largest
pattern and quantify the approaches towards and away number of additional spikes occur for the (purple) ellip-
from the depolarization block state of steady, elevated ticpathwithd=0.2. Althoughthesespikesstartearlier
voltage, which occur after the first and second crossing within the burst, for larger d, such as d = 0.4 (red), the
of AH, respectively. The observed rate of contraction relatively short time from the second AH crossing to the
(expansion) is related to the stable (unstable) eigenval- secondSNICcrossinglimitsthenumberofspikesthatcan
uesobtainedfromthe linearizationofthefastsubsystem occur within stage (v). For smaller d, such as d = 0.1
9
(blue),theneedforadditionalexpansiontoovercomethe quencywithinthebursts;thisrelationcanbeanalyzedby
strong contraction from stage (iv) is so strong that not computing contours of equal period in the region where
evenafullspikecanbefiredwithinstage(v)beforeSNIC (attracting) periodic orbits of the fast subsystem exist.
is crossed and the active phase ends. Moreover, we found that abrupt transitions from slow
spikes (at the onset of spiking) to faster spikes can be
predictedthisway. Similarly,theeigenvaluesoftheequi-
V. DISCUSSION libria of the fast subsystem that are encountered along
thepathareresponsiblefortheobservedcontractionand
expansionrates associatedwith the depolarizationblock
We introducedanddemonstratedaconceptualapproach
stateofDBbursting;detailscanbeanalyzedbycomput-
to modeling and analysis of fast-slow dynamics that al-
ingcontoursofequalrealpartsofsucheigenvaluesinthe
lows one to determine and extract quantitative informa-
relevant regions of the (Ca,Na)-plane. These contours
tion of relevance to the appplication at hand. Our ap-
explain an apparent and initially counter-intuitive non-
proachappliestofast-slowsystemswithatleasttwoslow
monotonicityinthenumberofspikesthatariseattheend
variables. The qualitativefeaturesofthe systemdynam-
oftheDBphase,justbeforetheonsetofthesilentphase,
ics arethen givenby the bifurcationstructure ofthe fast
with respect to the aspect ratio d of the imposed ellipti-
subsystem and how it is encountered by a closed path
cal path in the (Ca,Na)-plane. We find similar results
in the space of slow variables. Obtaining information on
(not reported here) with the models proposed in [3] and
thefastsubsystemintheregionsthataretraversedalong
in [9]. Given their connection to the well-established re-
this path, suchas periods ofperiodic orbits andcontrac-
lationshipbetweenbifurcationstructuresofthe fastsub-
tion rates associated with invariant objects, allows one
system and qualitative features of the bursting pattern,
to characterize how specific quantitative features of the
webelievethatourideaswillextendnaturallytomodels
observed solution arise. Moreover, we specify a family
with different underlying bifurcation diagrams of their
of paths in the state space of slow variables suitable for
associated fast subsystems; for example, our approach
the construction of an associated driven system, which
might be of use for modeling and analysis of the sound
we use to investigate systematically these quantitative
patterns produced by songbirds [1, 2].
solution characteristics.
As our specific test-case example we considered depo-
larization block or DB bursting as displayed by a seven- The linkage that we have established between quan-
dimensional model from [21] for neurons of the respira- titative bursting features and path of traversal through
tory brain stem. Here, the concentrations Ca of calcium the (Ca,Na)-plane opens the door to designing an im-
andNaofsodiumintheneuronevolveonamuchslower posedpathfortheslowvariablesCaandNatoachievea
time scale than the voltage potential v across the cell particular quantitative outcome. This ability to control
membrane, which is the observable in experiments. We solutionfeaturesisveryusefulfromamodeldevelopment
derived a four-dimensional reduction of that model with perspective, especially in situations where it is advanta-
onlytwofastvariablesandalsoconsidereditsdrivenver- geous to work with a lower-dimensionalmodel reduction
sion, where the family of paths was chosen to be ellipses that must be tuned to capture complicated dynamics.
in the plane of Ca and Na. DB bursting is organized The design of paths can guide model development and,
qualitatively by a curve of saddle-node bifurcations on moreover, it can be used for parameter estimation and
aninvariantcycle (SNIC) anda curveofAndronov–Hopf fitting for the slow dynamics, which may be particularly
bifurcations(AH)ofthefastsubsystem,andellipticpaths useful in light of difficulties in experimentally measuring
naturally capture the crossing sequence arising in the quantities associated with slow variables in neural mod-
biological model. The bursting pattern traces a family els. Once a family of paths has been chosen, the path or
of attracting periodic orbits of the fast subsystem, and pathswiththemostsuitablequantitativefeaturescanbe
also exhibits contractiontowardsand subsequentexpan- selected. We used ellipses here, but other, more compli-
sionawayfromtheequibriathatundergotheAndronov– cated paths may be considered, for example, to achieve
Hopf bifurcation. Although there are some quantita- a closer fit between the voltage dynamics of the reduced
tive discrepancies between the active phase dynamics in four-dimensionalmodelandthetargetseven-dimensional
the seven- and four-dimensionalbiological models, a key model. An important point is that the family of paths
point is that by tuning imposed paths, we can achieve a is specified by a resonably low number of parameters.
more complete quantitative agreement. In future work, we plan to consider the utility of our
More specifically, we found that the amount of time proposedapproachfor parameter estimationof fast-slow
an imposed path spends in different regions of the bursting models by coupling it with optimization tech-
(Ca,Na)-plane,correspondingtopassagebetweenbifur- niques. One challengeistodefine suitablefunctions that
cation curve crossings, is only one factor in determin- encapsulate, in terms of relevant quantitative features,
ing timing-related quantitative features of the resulting the distance of a given (periodic) output from the dy-
bursting pattern. In fact, the periods of the periodic or- namics under consideration; in applications, the latter
bits of the fast subsystem that are encountered along may be generated by an actual experiment, rather than
the path play a crucial role in shaping the spike fre- a higher-dimensional mathematical model.
10
ACKNOWLEDGMENTS where X ∈ {n,m,h}; the time scale τs(v) = τs for s
is taken constant. There are many parameters in this
model, and their values as used in this paper are listed
This work was partially supported by NSF awards
in Table I.
DMS 1312508 and DMS 1612913. JR thanks the Uni-
Thereducedtwo-dimensionalfastsubsystem(3)isob-
versity of Auckland for hosting an extended visit while
tained from the five-dimensional fast subsystem (1) by
part of this work was performed.
taking quasi-steady-state assumptions for m and s and
setting h = 1−1.08n. We use the same parameter val-
uesasinTableI,exceptthosespecifiedinTableII,which
APPENDIX were adapted to recover the bifurcation structure of the
fast subsystem and general slow subsystem behavior as
in system (1)&(2).
For completeness,we provideherethe precisedefinitions
of the functions used in system (1)&(2); the parameters
TABLE I. Parameter values used for the functions in the
are given in Table I. Recall that system (1)&(2) is the seven-dimensional model (1)&(2).
samemodelofDB bursting inneuronsofthe respiratory
conductances reversal half
brain stem as described in [21] and full details can also
(nS) potentials (mV) activations
be found there.
gL = 3.0 EL = −60.0 θh = −30.0mV
The equation for voltage v in system (1) includes the gNa = 150.0 ENa = 85.0 θm = −36.0mV
following currents: gK = 30.0 EK = −75.0 θn = −30.0mV
gsyn = 2.5 Esyn = 0.0 θs = 15.0mV
IL(v) = gL(v−EL), gCAN = 4.0 ECAN = 0.0 kCAN = 0.9µM
I (v,n) = g n4(v−E ), slopes time
K K K
(mV or µM) constants (ms) scaling constants
I (v,m,h) = g m3h(v−E ),
Na Na Na σh = 5.0 th = 15.0 kNa = 10.0 mM
Isyn(v,s) = gsyns(v−Esyn), σm = −8.5 tm = 1.0 Nab = 5.0 mM
g (v−E ) σn = −5.0 tn = 30.0 kCa = 22.5 ms−1
ICAN(v,ca) = CAN CAN , σs = −3.0 τs = 15.0 Cab = 0.05µM
1+exp((Ca−kCAN)/σCAN) σCAN = −0.05 kIP3 = 1200 µMms−1
I (Na) = r (φ(Na)−φ(Na )),
pump pump b other
C = 45.0pF k = 1.0 rpump = 200.0pA
where ε = 7×10−4 α = 6.6×10−5mMpA−1ms−1
Na3
φ(Na)= .
Na3+k3 TABLEII.Parametervaluesusedforthefunctionsinthetwo-
Na dimensional fast subsystem (3) that are different from those
used in thefive-dimensional fast subsystem (1).
The equations for n, m, h, ands all have a similar form,
with gK = 15.0nS θs = 10.0mV
gCAN = 10.0nS kCAN = 0.25µM
1.0 σs = −8.0mV kCa = 60.0ms−1
X∞(v)= 1.0+exp((v−θX)/σX), kIP3 = 1700 µMms−1
k = 10.0 rpump = 1500.0pA
where X ∈{n,m,h,s}, and ε = 0.005
t
X
τ (v)= ,
X
cosh v−θX
(cid:16) 2σX (cid:17)
[1] A. Amador and G. B. Mindlin, Eur. Phys. J. B 87, 300 Singapore, 2005).
(2014). [6] B. Ermentrout, Simulating, Analyzing, and Animating
[2] R.G.Alonso,M.A.Trevisan,A.Amador,F.Gollerand Dynamical Systems (SIAM,Philadelphia, PA,2002).
G. B. Mindlin, Front.Comput. Neurosci. 9, 41 (2015). [7] G. B. Ermentrout and N. Kopell, SIAM J. Appl. Math.
[3] E. Barreto and J. Cressman, J. Biol. Phys. 37, 361 46, 233 (1986).
(2011). [8] B.ErmentroutandD.H.Terman,MathematicalFounda-
[4] R. Bertram, M. J. Butte, T. Kiemel, and A. Sherman, tions of Neuroscience (Springer, NewYork,NY, 2010).
Bull. Math. Biol. 57, 413 (1995). [9] B. Ermentrout, J. Dutta-Moscato, and D. Pinto, in
[5] S. Coombes and P. C. Bressloff, Bursting: The Gene- Bursting: TheGenesisofRhythmintheNervousSystem,
sis of Rhythm in the Nervous System (World Scientific, edited byS.Coombes and P. Bressloff (World Scientific,