Scaling law for the transient behavior of type-II neuron models

0.43 MB·English

by  Miguel Angel Durán Roa

Scaling law for the transient behavior of type-II neuron models
∗
M. A. D. Roa and M. Copelli
Laborat´orio de F´ısica Te´orica e Computacional, Departamento de F´ısica,
Universidade Federal de Pernambuco, 50670-901 Recife, PE, Brazil
O. Kinouchi
Faculdade de Filosofia, Ciˆencias e Letras de Ribeir˜ao Preto, Universidade de S˜ao Paulo,
Avenida dos Bandeirantes 3900, 14040-901, Ribeir˜ao Preto, SP, Brazil
N. Caticha
7
Instituto de F´ısica, Universidade de S˜ao Paulo 05508-090, S˜ao Paulo, SP, Brazil
0
0
We study the transient regime of type-II biophysical neuron models and determine the scaling
2
behavior of relaxation times τ near but below the repetitive firing critical current, τ ≃ C(Ic −
n I)−∆. For both the Hodgkin-Huxley and Morris-Lecar models we find that the critical exponent
a is independent of the numerical integration time step and that both systems belong to the same
J
universality class, with ∆ = 1/2. For appropriately chosen parameters, the FitzHugh-Nagumo

4 model presents the same generic transient behavior, but the critical region is significantly smaller.
2 Wepropose an experiment that may reveal nontrivial critical exponentsin thesquid axon.

] PACSnumbers: 05.45.-a,89.75.Fb,05.90.+m
C
N
. I. INTRODUCTION croscopic details involved. Full characterization of these
o
bifurcation routes is important for a deeper understand-
i
b ing of how they affect the firing behavior and possible
- Thesearchforbiophysicalmodelsforinformationpro- implementation of neural codes. For example, it is now
q
[ cessing systems has led to a variety of model neurons acknowledgedthatbistablebehaviorandthesmallrange
which describe the dynamics of the membrane poten- of firing frequencies in neurons that undergo subcritical

2 tial. Collective properties arise from their interaction Hopf bifurcations prevent a robust use of a pure fre-
v through several possible architectures and types of cou- quency code.
9
plings. These can be chemical, voltage-gated synapses,
1 The transient behavior of neuron models has not re-
simpler proteic electrical connections, or even just elec-
0 ceived much attention in either experimental or theoret-
trical ephaptic interactions arising between neighboring
6 ical studies. In this paper we study the divergence of
0 nerve fibers. Maybe the most striking feature of a neu- transient times in a class of conductance-based models
6 ron is the threshold of the stimulus that separates spik-
andshowthatitfollowsauniversalcriticalbehavior. We
0 ing from nonspiking regimes. Spiking neurons generate
proposeanexperimentthatmaytestourtheoreticalpre-
/ taletale signatures which have served as the basis for
o dictionsanddiscusshowneuronscouldemploytransients
i frequency-dependent neural codes, an idea that can be for computational purposes.
b
traced back to the work of Adrian in the 1920s [1]. Al-
-
q though of paramount importance to neural dynamics,
: spike frequencies do not tell the complete story. Sub-
v II. TRANSIENTS IN TYPE-II MODELS
tle computations may arise from subthreshold dynamics
i
X suchasforexamplesintheearlystagesofthemammalian
r visual system, olfactory bulb, and cortex. In many cases The Hodgkin-Huxley (HH) model is a biophysically
a
thekeytotheinformationdynamicsliesinthetransients, motivated system of four coupled nonlinear differential
eitherbelowthe currentthresholdto generateactionpo- equations that describe the dynamics of the membrane
tentials orthe thresholdto generateinfinite sequencesof potential V of the squid giant axon (e.g., [2, 3]):
spikes. Inthispaperweinvestigatetransientspiketrains
of single model neurons since this might have a bearing
on the collective behavior, i.e., computational capabili- dV
C = G m3h(E V)+G n4(E V)
ties, of subthreshold assemblies of neurons. Na Na K K
dt − −
A dynamical system approach reveals that universal +GL(EL V)+I(t),
−
bifurcation scenarios for the firing behavior appear irre- dxi
= α (V)(1 x ) β (V)x . (1)
spective of the specific membrane ion channels and mi- dt xi − i − xi i
x stand for the three gate variables x = m,h, and n
i i
describing the activation of ionic channels and α , β
xi xi
∗Electronicaddress: [email protected] are voltage-dependent transition rates [3]:
2
Hodgkin-Huxley
2.5 0.1V
α (V) = − ,
m e(2.5−0.1V) 1
β (V) = 4e−V/18, − I = 6.1 µA/cm2
m
αh(V) = 0.07e−V/20, V) 50 mV
m
βh(V) = e(3−0.11V)+1, V ( I = 6.25 µA/cm2
0.1 0.01V
α (V) = − ,
n e(1−V) 1
−
βn(V) = 0.125e−V/20, (2) I = 6.262 µA/cm2
0 50 100 150 200 250 300 350 400 450 500
where voltages are expressed in mV, rates in ms−1,
time (ms)
C = 1 µF/cm2 is the specific membrane capacitance,
GNa = 120 mS/cm2, GK = 36 mS/cm2, and GL = FIG. 1: Examples of transient behavior near Ic for the HH
0.3 mS/cm2 are ionic conductances per unit area, and model. The constant step current is applied at t = 10 ms.
E =115 mV, E = 12 mV, and E =10.6 mV are TheestimatedcriticalcurrentisIc =6.26422125685 µA/cm2
Na K L
reversalpotentials. − for an integration timestep dt=0.01 ms.
The HH system plays a fundamental role in the field
ofneurophysiologyandcomputationalneurosciencesince
it defines the class of conductance-based models. As a
function of a constant applied current I, the HH model
undergoesasubcriticalHopfbifurcationatI =I ,above
H
dV V +1
which the fixed point solution (membrane potential at C = 0.5G 1+tanh (E V)
Ca Ca
rest) is no longer stable and trajectoriesare attracted to dt (cid:20) (cid:18) 15 (cid:19)(cid:21) −
a stable limit cycle, leading to repetitive firing (infinite +G w(E V)+G (E V)+I(t),
K K L L
− −
train of action potentials). dw
= 0.1cosh(V/60)[1+tanh(V/30) 2w] , (3)
The sudden jump to a periodic behavior with nonzero dt −
frequency f is analogousto afirst-orderphase transition
andis referredto as type-II behaviorinthe neuroscience when, for example, the values of the parameters are
literature [4]. As in first-order phase transitions, coex- chosen as G = 1.1 mS/cm2, G = 2.0 mS/cm2,
Ca K
istence also appears in type-II behavior. Just below the G = 0.5 mS/cm2, E = 100 mV, E = 70 mV,
L Ca K
Hopf bifurcation, the stable fixed point coexists with a and E = 50 mV. Then, large transient time−s are also
L
stable limit cycle, their basins of attraction being sepa- observed (F−ig. 2).
ratedbyanunstablelimitcycle[4]. Bothlimitcyclesare
The FitzHugh-Nagumo (FHN) system
created in a saddle-node (or “fold”) bifurcation of cycles
at I = I < I . In our analogy with equilibrium phase
c H
transitions, I would correspond to a spinodal point. If
c
dV
the fixed point at I = 0 is perturbed by the application = V(V a)(1 V) w+I,
of a constant current near but below I , several spikes dt − − −
c
mayappearbeforethe systemreturnsto the newresting dw
= ǫ(V γw), (4)
state (Fig. 1). dt −
In the original version of the Morris-Lecar (ML)
model [5], a system with two coupled nonlinear ordinary has been proposed as a low-dimensional toy model that
differential equaions (ODEs) used to describe action po- represents the type-II behavior of the HH and other ex-
tentials in a barnacle motor fiber, the relevant bifurca- citable systems. We verified that the transient behavior
tionatthe onsetofrepetitive firingisa saddle-nodeone. here reported is not seen with the usual parameters [6].
That means that the spiking frequency varies continu- However, the FHN model can reproduce the HH tran-
ously from I as f (I I )β, with β = 1/2, which sient behavior if one chooses parameters near a = 0.5,
c c
is similar to a mean∝field−second-order phase transition γ =4.2, and ǫ=0.01 (Fig. 3).
behavior if we think of f as the order parameter. This In this paper we show that the long relaxation times
transition is called type-I behavior in the neuroscience in type-II models are a consequence of the changes in
literature [4] and does not present a slow transient phe- phase space which occur near the creation of the limit
nomenon similar to Fig. 1. However, the Morris-Lecar cycles. Moreover, this leads to a scaling relation whose
systemhasalsobeenusedtodescribecellswhichpresent exponent can be predicted, in agreementwith numerical
a type-II behavior [3, 4], which occurs for the equations simulations.
3
5
10
Morris-Lecar
∆ = 0.47
I = 24.8411 µA/cm2
4
10
50 mV
mV) ms) HH
V ( I = 24.8413 µA/cm2 τ (
3
10
5
T =10 ms
max 4
T =10 ms
I = 24.841337 µA/cm2 2 max
10
-9 -8 -7 -6 -5 -4 -3 -2 -1
10 10 10 10 10 10 10 10 10
0 50 100 150 200 250 300 350 400 450 500
time (ms) Ic - I (µA/cm2)
FIG. 2: Examples of transient behavior near Ic for the FIG. 4: Relaxation times for the HH model as a function of
ML system. The constant step current is applied at thedistancetocriticalcurrentfordifferentintegrationtimes:
t = 10 ms. The estimated critical current is Ic = Tmax = 104 ms (open circles) and Tmax = 105 ms (filled
24.84134676279 µA/cm2 for an integration time step dt = circles).
0.01 ms.
5
10
dt=0.01
dt=0.1
FitzHugh-Nagumo
∆ = 0.49
4
10
I = 0.10254
V (arb. units) 0.5 (arb. units) I = 0.10254471 τ (ms) 103 ML
2
10
I = 0.102544718 -9 -8 -7 -6 -5 -4 -3 -2 -1
10 10 10 10 10 10 10 10 10
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Ic - I (µA/cm2)
time (arb. units)
FIG. 5: Relaxation times for the Morris-Lecar model as a
FIG. 3: Examples of transient behavior near Ic for the FHN function of the distance to critical current for different inte-
system. The constant step current is applied at t=10. The gration times: dt = 0.01 ms (filled circles) and dt = 0.1 ms
estimatedcriticalcurrentisIc =0.1025447183127 foraninte- (open circles).
gration timestepdt=0.01(allquantitiesinarbitraryunits).
nents for these three type-II biophysical neuron models,
III. SCALING LAW
finding very good agreement between them.
We integrated the equations using a standard fourth-
Therelaxationtimeτ maybedefinedasthetimeuntil orderRunge-Kuttaalgorithmanddetermined τ by mea-
the last spike, or the time until the membrane voltage suringthetimeintervalfromtheonsetofthecurrentstep
stays within a small distance from the resting potential to a near stop of the flow [x˙ <10−5, where x˙ is the ve-
(these twotimes areverysimilarnearI ). Whenweplot locity vector in phase spac|e:|x˙ = (w˙,v˙) for the ML and
c
τ asafunctionofI I wefindapowerlawdivergenceof FHN systems, and x˙ = (w˙,m˙,h˙,n˙) for the HH model].
c
the relaxationtime,−τ C(I I)−∆,where ∆issimilar As opposed to the Hopf bifurcation, the fold bifurcation
c
≃ −
to a dynamic criticalexponent. The ∆ exponent charac- cannot be obtained analytically, so I was estimated nu-
c
terizes the “critical slowing down” behavior near the bi- merically after integration of the ODEs up to a (long)
furcation. We expect that ∆ is a universalexponent but maximum time T . The determination of the criti-
max
we are not aware that this exponent has been measured cal current is sensitive to T , but in practice this only
max
for neuron models. Here we report the measured expo- limits the range of validity of the power law (see Fig. 4).
4
5
10 0.6
a) I < I < I
c H
∆ = 0.48 0.5
s) 0.4 UPO
nit
b. u 104 FHN w 0.3 SPO
r
τ (a 0.2 SFP
0.1
103 0
10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 -60 -50 -40 -30 -20 -10 0 10 20 30
I - I (arb. units)
c V
FIG. 6: Relaxation times for the FitzHugh-Nagumo model 0.6
as a function of the distance to critical current. Note that b) I < I
c
the power law becomes visible only very close to the fold 0.5
bifurcation (Ic−I ∼10−9).
0.4
w 0.3
We have employed T = 105 ms and dt = 0.01 ms,
max
unless otherwise stated. The estimated critical currents 0.2
I (T ,dt) quoted in the figure captions are very pre-
c max
cisely determined given Tmax and the integration step 0.1 SFP (I=0)
size dt.
0
The criticalexponentis determined fromthe plotτ vs
-60 -50 -40 -30 -20 -10 0 10 20 30
I I. Wefound∆=0.47fortheHHsystem(Fig.4)and
c
−
∆ = 0.49 for the ML model (Fig. 5), irrespective of the V
size of the integration step dt. This suggests a universal
exponent ∆ = 1/2. Obtaining long transients in the FIG.7: PhaseportraitsoftheML modelfor I slightly above
FHN model has proved numerically more difficult, since (a) and below (b) Ic. The long transient is dominated by
the phenomenon occurs only very close to I (Fig. 6). the time it takes to pass through the region where the limit
c
Nonetheless, we have obtained the exponent ∆=0.48. cyclesareabouttoemerge. SPO(UPO)=Stable(Unstable)
Figure 7(a) shows that for I . I < I an unstable Periodic Orbit, SFP = Stable Fixed Point. Horizontal axes
c H
(V) in mV and vertical axes (W) are dimensionless.
limit cycle and a stable one coexist and surround a sta-
ble fixed point. The fixed point for I = 0 lies outside
both limit cycles [Fig. 7(b)]. Therefore, when the cur-
rent is abruptly changed to I . I , the fixed point is
c
displaced to a region within a ghost limit cycle, and the beenobservedbecausethemodelhasbeenstudiedinthe
transient is completely dominated by the time it takes absence of an external current.
for the system to overcome it. The ghost is a natural It should be clear that the current step is just a sim-
consequence of the system being immediately below a ple way of putting the system outside the ghost limit
saddle-node bifurcation of cycles, and can be character- cycle, but the scaling law for the transient is not re-
ized by the vanishingly small flow component normal to strictedtothissomewhatartificialprotocol(eventhough
the half-stablelimitcycle thatiscreatedatI =Ic. Itef- it is very common in both theory and experiments).
fectivelyworksasaone-dimensionalextendedbottleneck Close to the fold transition, any short-lived perturba-
through which the system must pass before reaching the tion that is strong enough to make the system cross
fixed point (see Fig. 8 for a caricature). If one considers the ghost limit cycle will give rise to long transients
an analogous system in polar coordinates [7] θ˙=f(r,θ), back to the fixed point. We exemplify this with a bi-
r˙ = µr+r3 r5, where f(r,θ) > 0 , r > 0, it is clear ologically plausible example in the HH model. Suppose
− ∀
that for µ.µ = 1/4 the time for the system to over- the system is somehow maintained close to criticality at
c
come the ghost at−r = 1/2 scales as τ (µ µ)−1/2 I . I (this could be achieved by several different pos-
c c
∼ −
(the Hopf bifurcation occurring only at µ = 0). Solu- sible mechanisms, so we just fix I). In addition, assume
H
tions for the transient behavior have been obtained by the neuron undergoes a fast excitatory postsynaptic po-
Tonnelier[8]forMcKean’spiecewiselinearversionofthe tential(EPSP)simulatedbyaninjectedsynapticcurrent
FHN model [9]. However, the scaling behavior has not I (t)=g (t)(E V):
syn syn Na
−
5
θ(t′)(g t′/τ2)exp( t′/τ ),whereθistheHeavisidefunc-
V SPO tion, gm = 1s mS/c−m2, τs = 0.5 ms, and t′ = t t ,
m s EPSP
−
where t =50 ms is the time the EPSP is initiated.
EPSP
UPO Results are shown in Fig. 9, whose similarity with Fig. 1
UFP
attests to the robustness of the effect (notice however
that, differently from Fig. 1, in Fig. 9 the resting mem-
brane potential is at the I = 0 fixed point before the
6
perturbation). This opens interesting possibilities from
UPO
the point ofview of neuronalcomputation: the lengthof
τ
SFP thetransientresponseofaneuron“probed”byanEPSP
SPO could code for its internal state of excitability, that is,
for how close it is to I . Naturally, this transient cod-
c
ing would work only if the system is close to the fold
bifurcation. We thereforecouldhaveanotherexample in
I neuroscience of optimal information processing at criti-
cality [11, 12, 13, 14, 15, 16, 17].
I I Itis interestingtopointoutthatthis bottleneck effect
c H is analogousto whatoccursfor type-I neuronsabovethe
saddle-nodebifurcationthatleadsto repetitivefiring. In
FIG. 8: Schematic bifurcation diagram for type II neuron thatcase,however,theghostresultsfromtheanihilation
models. The fold bifurcation occurs at Ic, while the Hopf of fixed points (not limit cycles) and the period T of the
bifurcation occurs at IH. Owing to the onset of the current limitcyclesdivergesas(I Ic)−1/2. Thisprovidesacom-
step (dot-dashed arrow), the fixed point for I = 0 becomes plementary scenario conn−ecting both classes of neurons:
the initial condition in a new phase portrait. The transient thetransientoftype-IImodelsbelowI divergeswiththe
c
τ (solid arrow) to reach the new fixed point is governed by
same exponent as the period of type-I models above I ,
theghost limit cyclewheretheUPOand theSPOannihilate c
that is, ∆=β.
each other (see Fig. 7).
IV. CONCLUDING REMARKS
gsyn 1 mS/cm2 beWhaeviwoerrfeorutnraabnlseietnotfisnindnaeduersocnrimptoidoenlsoofrththisesacsasloincgialtaewd
40 45 50 55 60 dynamic critical exponent in the literature. Some kind
of critical slowing down for subthreshold oscillations has
I = 6.1 µA/cm2
V) been reportedexperimentallyin the squid axon[18], but
m 50 mV
V ( these authorsexaminedthe vicinity ofa parametricsub-
criticalHopfbifurcation,notasaddle-nodebifurcationof
I = 6.25 µA/cm2 cycles induced by external currents. We propose that a
similar experiment with high-precision injected currents
near I could be used to check the power law found in
c
I = 6.262 µA/cm2 the computational model. Since the area of the giant
squid axon is of order 1 cm2, to examine the critical
0 100 200 300 400 500
∼
time (ms) regime requires that current fluctuations should be less
than 103pA(seeFig.4). Evenifthefullcriticalregime
∼
FIG. 9: Long transients appear if the system initially at rest seems to be hard to achieve,the initial divergence in the
with I . Ic is perturbed by an additional EPSP (see text transient lifetime may provide an experimental check of
for details). Solid lines are membrane potentials, while the our predictions.
dashed line is thesynaptic conductance.
We emphasize that both in standard experiments and
inoursingle-compartmentmodelspaceclampingisused.
It might happen though, as in spin systems, that for
the extendedrealsystemwithoutvoltageclamp,orfor a
compartmental model with a large number of compart-
dV ments, the exponents may differ from the values here
C = G m3h(E V)+G n4(E V)
dt Na Na− K K − reported, changing the universality class to one not de-
+G (E V)+I +I (t). (5) scribed by a “mean field” approach. Interestingly, this
L L syn
−
would mean that collective properties within a nontriv-
The fast change in the synaptic conductance is ial universality class could be observed at the level of a
′
given by Rall’s alpha function [10]: g(t) = singleaxon. Sowhetherthisresultcanbe verifiedexper-
6
imentally hinges on the effects that spatially extended jection” in type-II neuron models. For instance, in the
neurons may have on the robustness of this picture. phenomenonofcoherenceresonance[23]noiseitselfplays
Furthermore, noise could always play a role. In stud- this reinjecting role. However,the excitable systems em-
ies of type-I intermittency in simple maps, the length of ployedareusuallynotcloseenoughtothe foldtransition
the “laminar phase” l in a chaotic regime is analogous to exhibit long transients, so it would be interesting to
to the transient in thhisiwork and diverges as l ǫ−1/2 investigate the effects of the scaling laws we report here
h i∼
because of a zero-dimensionalbottleneck as the distance in the resonance curves. These theoretical issues should
ǫ to a tangent bifurcation tends to zero [19, 20]. In the be dealtwith before engagingin anexperimentalsearch.
presenceofadditivenoisewithamplitudeg [21],thescal- Acknowledgments: The research was supported
ingchangesto√ǫ l f(g2/ǫ3/2),wheref isauniversal by FACEPE, CNPq and special program PRONEX.
h i∼
function (see also [22] for recent extensions). It is con- M.A.D.R. thanks the organizers of LASCON, the Latin
ceivable that similar scaling relations could be obtained American School on Computational Neuroscience, for fi-
for τ provided that the chaotic phase of intermittency nancial support. M.C. thanks J. R. Rios Leite and J. R.
theory could be replaced by some mechanism of “rein- L. de Almeida for inspiring discussions.
[1] E. D.Adrian, J. Physiol. (London) 61, 49 (1926). Natl. Acad.Sci. USA 97, 3183 (2000).
[2] D. Johnston and S. M.-S. Wu, Foundations of Cellular [12] J.M.BeggsandD.Plenz,J.Neurosci.23,11167(2003).
Neurophysiology (MIT Press, Cambridge, 1995). [13] D. R.Chialvo, Physica A 340, 756 (2004).
[3] C.Koch,Biophysics of Computation (OxfordUniversity [14] V. M. Egu´ıluz, D. R. Chialvo, G. A. Cecchi, M. Ba-
Press, NewYork, 1999). liki,andA.VaniaApkarian,Phys.Rev.Lett.94,018102
[4] J. Rinzel and B. Ermentrout, in Methods in Neuronal (2005).
Modeling: From Ions to Networks, edited by C. Koch [15] L. S. Furtado and M. Copelli, Phys. Rev. E 73, 011907
and I. Segev (MIT Press, Cambridge, MA, 1998), pp. (2006).
251–292. [16] O. Kinouchiand M. Copelli, Nat. Phys. 2, 348 (2006).
[5] C. Morris and H. Lecar, Biophys. J. 35, 193 (1981). [17] D. R.Chialvo, Nat. Phys. 2, 301 (2006).
[6] J. Rinzel, in Research Notes in Mathematics Nonlinear [18] G. Matsumoto and T. Kunisawa, J. Phys. Soc. Jpn. 44,
Diffusion,editedbyW.E.Fitzgibbon andH.R.Walker 1047 (1978).
(Pittman Press, London, 1977), pp.186–212. [19] P.MannevilleandY.Pomeau,Phys.Lett.75A,1(1979).
[7] S.H.Strogatz,NonlinearDynamicsandChaos: withAp- [20] Y. Pomeau and P. Manneville, Commun. Math. Phys.
plications to Physics, Biology, Chemistry and Engineer- 74, 189 (1980).
ing (Addison-Wesley,Reading, MA,1997). [21] J.E.Hirsch,B.A.Huberman,andD.J.Scalapino,Phys.
[8] A.Tonnelier, SIAMJ. Appl.Math. 63, 459 (2002). Rev. A 25, 519 (1982).
[9] M. P. McKean, Adv.Math. 4, 209 (1970). [22] H. L. D. de S. Cavalcante and J. R. Rios Leite, Phys.
[10] H.R.Wilson,Spikes, Decisions and Actions: Dynamical Rev. Lett.92, 254102 (2004).
Foundations of Neuroscience (Oxford University Press, [23] A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775
Oxford,1999). (1997).
[11] S. Camalet, T. Duke, F. Ju¨licher, and J. Prost, Proc.