A simple model for 1/fα noise ∗ J¨orn Davidsen and Heinz Georg Schuster Institut fu¨r Theoretische Physik und Astrophysik, Christian-Albrechts-Universit¨at, Olshausenstraße 40, 24118 Kiel, Germany 2 (February 1, 2008) 0 0 We present a simple stochastic mechanism which generates pulse trains exhibiting a power law 2 distributionofthepulseintervalsanda1/fαpowerspectrumoverseveraldecadesatlowfrequencies n withαclosetoone. Theessentialingredientofourmodelisafluctuatingthresholdwhichperforms a a Brownian motion. Whenever an increasing potential V(t) hits the threshold, V(t) is reset to the J origin and a pulse is emitted. We show that if V(t) increases linearly in time, the pulse intervals 3 can be approximated by a random walk with multiplicative noise. Our model agrees with recent 1 experimentsinneurobiologyandexplainsthehighinterpulseintervalvariabilityandtheoccurrence of 1/fα noise observed in cortical neurons and earthquakedata. ] h c 5.40.-a, 87.10.+e, 89.75.Da e m Theomnipresenceof1/fα noiseinnatureisoneofthe Such a model is often used to describe single neurons: oldestpuzzlesincontemporaryphysicsstilllackingagen- V(t) is the membrane potential and the emitted pulse - t erally accepted explanation. The phenomenon is charac- is the generated action potential. However, it is usually a t terizedbyacertainbehavioroftherespectivetimesignal: assumed that the threshold is constant in time. Recent s . The powerspectrumS(f)is proportionalto 1/fα atlow investigationhaveshownthatthis isnottrueforcortical t a frequencies f with α 1. Examples include the light of neuronsin vivo [16]andin vitro [17]. Additionally,there m quasars[1],electricalm≈easurements[2],musicandspeech is evidence that the spike trains of auditory neurons [18] - [3],humancognition[4]andcoordination[5],thecurrent and of neurons in the mesencephalic reticular formation d through ion channels [6], network traffic [7], burst errors [19] are not renewal, i.e., successive time intervals be- n in communication systems [8], freeway traffic [9], gran- tweenspikesarecorrelated. Thesefactsareincorporated o c ular flow [10], etc. The time signals of a large number inourmodelinthe simplestpossibleway. Moreover,the [ of these systems [5,7–10] resemble a pulse train consist- effectsofdeadtimeorabsoluterefractoriness,whichlim- ing of individual, largely identical events which occur at its the rate atwhicha neuroncanfire,areautomatically 1 v discrete times. This is especially true for spike trains of included in our model via the lower bound Cl. The up- 8 singlenervecellsforwhich1/fα noisehasbeenobserved perboundC preventsaninfinitetimedifferencebetween u 9 in various brain structures [11–15]. The reported expo- two pulses. 1 nents, which depend both on the presence or absence of 1 a sensory stimulus [11,12] and on the state of the ani- 1.4 0 2 mal (REM sleep vs awake state) [13–15], vary from 0.68 1.2 C(t) 0 to 1.38. The power-law behavior for spike train power t/ spectra lies within the range 0.01 to 10 Hz, extending 1 a typically over 2 decades. In almost all cases the upper m 0.8 limit ofthe observedtime overwhichfractalcorrelations - exist is imposed by the duration of the recording. d 0.6 n InthisLetter,weproposeasimplemechanismforgen- V(t) o erating pulse trains with 1/fα behavior in systems with 0.4 c a threshold-controlled dynamics like, e.g., neurons and : 0.2 v earthquake faults. Our model is based on an integrate- i and-fire (IaF) mechanism and consists of a single unit 0 X characterized by two variables (see Fig. 1): The voltage 0 1 2 3 4 5 ar V(t) and the threshold C(t). Initially, the voltage is be- t low the threshold. Then, the voltageincreasesmonoton- FIG.1. Dynamicsofourmodelforalinearincreasingvolt- ically in time — in the simplest case just linearly — and age V(t) with V0 = 0. The dashed lines represent the lower the threshold evolves according to a Brownian motion and theupperboundary for thefluctuatingthreshold C(t). with diffusion constant D within reflecting boundaries As one can see from Fig. 2, the power spectrum of V < C < C(t) < C . As soon as V(t) has reached the 0 l u the pulse train generated by our model with linear in- threshold, the voltage is reset to V and a pulse of unit 0 creasingvoltageshowsa1/fα decayoverseveraldecades heightis emitted. Inthis way,apulse trainis generated. [20]. The frequency below which white noise behavior Note that the threshold is not reset to its initial value. 1 is observed is determined by Cu Cl and goes to zero 10−2. − for C C . The exponent α is not universal and increause−s ifl →the∞ratio (C V )/√D increases. In the 10−3 (a) l 0 − limit (C V ) √D, we find α = 0.5. This result is l 0 (b) − ≪ −4 explained by the fact that this limit corresponds to the 10 caseinwhichthewaitingtimesbetweenpulsesarepurely ) determined by the diffusive dynamics of the threshold. S(f10−5 (c) Hence, each waiting time is merely the first return time −6 10 ofa Brownianmotionwhich obeys a power-lawdistribu- tion with exponent 1.5. This implies α=0.5 [21]. −7 − 10 Tocompareourresultswithrealneurons,considerthe parameters of curve (b): The dead time of a neuron is −5 −4 −3 −2 −1 0 10 10 10 10 10 10 typically of order of milliseconds and the maximal time f / f difference between two spikes of the order of seconds. Nyq. This is exactly the ratio between Cl and Cu. Moreover, FIG. 2. Power spectrum of the pulse train generated by wecannowidentify oneunitoftime inourmodelwith5 our model. Parameters are (a): V0 = 0, Cl = 0.2, Cu = 40, milliseconds real time. Hence, the 1/fα behavior in (b) √D=0.2. (b) as (a) with Cu =200. (c): V0 =0, Cl =0.02, is found for f < 11 Hz. This and the fact that α 1.02 Cu = 4000, √D = 2.0. This curve is shifted down by 1 reproduces the experimental results very well. ≃ decade. Allcurvesshowaclear1/fα decay. αvariesfrom0.6 Extensive numerical investigations have shown that to 1.1. our model is very stable with respect to variations of the dynamics. Different forms of voltage increase, e.g., The FPDTF ofourmodelwith linear voltageincrease a linear, a squared, or a square-root increase, give sim- can be obtained by mapping the model to an IaF model ilar results. The assumption of a monotonical increase with constant threshold C = (Cu +Cl)/2: The voltage of the voltage can also be dropped. A small amount of V˜(t) is defined as the sum of V(t) and C C(t). This noisecanbe addedto the voltagesignalwithoutaltering means that V˜(t) fluctuates around V(t). It−also implies our findings. Substituting the reflecting boundaries by thatV˜(t)isresettoV0+C C(tk)afterthekththreshold − a confining potential does not change our results, either. crossing. Hence, the correlationsare now encoded in the This points towards a generic behavior. fluctuating reset. In conclusion, V˜(t) behaves almost as To obtain the power spectrum of our model analyti- a Brownian motion with drift to an absorbing barrier cally, consider the signal X(t) = δ(t t ) where t with a reset following each barrier crossing. The only denotesthetimeofoccurrenceofthPekkthp−ulske. Itfollowks difference is that the stochastic process is restricted to as shown in [22] the interval [V(t) + C Cu,V(t) + C Cl] at time t − − which makes an exact mathematical treatment difficult. 2 T Neglecting the restriction for a moment and setting C = 1 S(f)= lim (cid:12) dtX(t)exp−i2πft(cid:12) , (1) 0,weobtainfromthe associatedFokker-Planckequation T→∞2T (cid:12)(cid:12)Z (cid:12)(cid:12) the FPTDF (cid:12)−T (cid:12) (cid:12) (cid:12) = lim 1 (cid:12) expi2πf(tk+q−t(cid:12)k). (2) g (τ V˜(0))= −V˜(0) exp (V˜(0)+τ)2, (5) T→∞2T XX 1 | √2πDτ3/2 − 2Dτ k q With I¯=limT→∞ 21T(kmax−kmin+1), this leads to and g1 ≡ 0 for τ < 0. Here, V˜(0) is the initial dis- tance. Eq. (5)is,ofcourse,justanapproximationtothe S(f) I¯ expi2πf(tk+q−tk) , (3) FPTDFofourmodel. Thiscanalreadybe seenfromthe ≈ Xq D E fact that the FPTDF of our model is identically zero for τ < (C V ) and τ > (C V ). However, as will be l 0 u 0 where denotes the average over the ensemble and − − h···i shownbelowthis approximationprovestobe veryuseful over k. Hence, we need to know the probability distri- and Eq. (5) can even be simplified to a Gaussian with butions Ψ (τ)dτ of the time differences between pulses q the same mean and variance: τ = t t for all integers q. Ψ is merely the qth k k+q q passage−times density function g averaged over the sta- (V˜(0)+τ)2 q g (τ V˜(0)) exp . (6) tionary probability distribution of “initial” states V(0) 1 | ∝ − 2D V˜(0) | | Ψq(τ)= gq(τ V(0)) V(0). (4) FromtheFPTDFonecouldnowcomputehigherpassage h | i times density functions g in principle by convolution Since g can be computed from the first passage times n q density function (FPTDF), we will firstfocus onthe lat- τ ter. gn(τ|V˜(0))=Z g1(t|V˜(0))gn−1(τ −t|t)dt. (7) 0 2 However,thisconvolutionresistsananalyticaltreatment We have to point out that the behavior of P(τ) de- due to the Markovian character of our process and the pendsonthespecifickindofvoltageincrease. Foravolt- form of Eqs. (6,7). The former manifests itself in the age increase with (t t )β with 0 < β < 2, we find last fact that the first passage time directly gives the initial P(τ) τβ−2 for ou−r model. Substituting the reflect- ∝ distance between the threshold and the voltage for the ing boundaries by a potential of the form a/x+bx2 for next passage problem. The reason for this is that the example, we still find power law tails in the ISI distribu- norm of the kth reset of the voltage V˜ equals the time tion. This is exactly what measurements show for corti- differencebetweenthekthandthe(k 1)thpulsedueto calneurons[25],especiallyinthemesencephalicreticular the linear increase of the voltage V˜ w−ith slope 1. Hence formation [19]. Hence, our model seems to be especially the reset depends on the last passage time only and the well-suited to explain the occurrence of 1/fα noise in wholestochasticpointprocessistotallydescribedbythe that formation. In general, IaF models can not explain FPTDF plus its Markovianproperty. thehighinterspikeintervalvariabilityexhibitedbycorti- calneurons[26]. A fluctuating thresholdasdescribedby . our model, however, can solve this long-standing issue. −3 10 ForrenewalprocessestheISIdistributionfunctionand the FPTDF are identical and completely describe the 10−4 process. This is the case for a standard IaF neuron with resetofthe voltageto the originafter eachbarriercross- S(f)10−5 ing whichis describedbythe FPTDFin Eq. (5). Such a modelcanneither explaina 1/fα signalnor a power-law −6 inverse Gaussian f −1.02 decay of the ISI distribution. In contrast to our model, 10 Usherandco-workersshowedthatfractalbehaviormight Gaussian be a consequence of the global activity dynamics of a model −7 networkofIaFneurons[27]. Duetoexperimentallimita- 10 10−5 10−4 10−3 10−2 10−1 100 tions, however, the link between 1/fα single-unit power f / f spectraandmacroscopicactivitydynamicsremains,asof Nyq. now, a conjecture. Another attempt to explain the phe- FIG.3. Powerspectrumofthetimesignalgeneratedbythe nomenonof1/fαnoiseisbasedonfractalandfractal-rate linearversionofourmodelandbythetwoapproximationsfor stochastic point processes [28]. Certain types of these thesame parameters as in Fig. 2, (b). processes properly characterize the statistical properties foundindifferentexperiments. However,inmanycasesit These properties enable us to interpret the stochastic is notcleara priori why the realsystemshouldgenerate process generated by the Gaussian approximation (6) as such a process. This is especially true for the clustering a random walk of the inter-spike-intervals (ISI) τ : k Poissonprocess which was applied to explain 1/fα noise in the mesencephalic reticular formation [14]. In con- τ =τ + Dτ η . (8) k+1 k k k p trasttothat,ourmodelprovidesasimple explanationof Here, η denotes the white noise source. Note the spe- 1/fα noise in integrate-and-fire systems in general and, k cial kind of multiplicative noise which distinguishes our hence, applies to neurons in particular. The crucial as- model from the one in [23]. The variance of the step sumption is a fluctuating threshold. Such a behavior is length is proportional to the current “position”. Hence, related to models presented in [29] and was already con- the origin is a fixed point of the random walk. To omit sideredin[25]toexplainthepowerlawdecayofneuronal this difficulty andinspiritofouroriginalmodel, wecon- ISI distribution functions. However,the stochastic point sidertherandomwalkτ tobeconfinedbytworeflecting process was still assumed to be renewal and could not k boundaries, i.e., 0 < (C V ) < τ < (C V ). The explain a 1/fα behavior of the power spectral density u 0 k l 0 − − power spectrum of the pulse-train generated by such a function. Consequently, the second crucial assumption random walk with t = t + τ shows a clear 1/fα of our model is the Markovian character of the process k+1 k k behavior with the same exponent as for our model (see in agreement with Refs. [22,23]. To directly verify our Fig. 3). This is also true if the new inter-spike-interval model for neurons in the mesencephalic reticular forma- is chosen from an inverse Gaussian distribution given in tion, one should measure the evolution of the reset from Eq. (5) with an initial condition depending on the last the threshold potential. interval. Hence, these approximations seem to be jus- Allthreemaincharacteristicsofourmodel,i.e.,a1/fα tified. This is further confirmed by the stationary ISI spectrum, a power-lawdecayof the ISI distribution, and distribution. For the random walk, the ISI distribution a Markovian dependence of the ISI, can not only be ob- function P(τ) is proportional to τ−1 [24]. Simulations served in nerve cells but also in other systems. For ex- show that this is in excellent agreement with our model ample, the pulse signal, defined by associating the oc- and the inverse Gaussian approximation. currence of pulses of unit height with earthquakesin the 3 Mojaveregion[30],shows1/fαfluctuationswithα=1.3 [8] J.M.BergerandB.B.Mandelbrot,IBMJ.Res.Dev.7, andhasanISI distributionwith exponent 1.1(Fig. 4). 224 (1963). Moreover, there is evidence that the ISI d−o not obey a [9] T. Musha and H. Higuchi, Jpn. J. Appl. Phys. 15, 1271 renewal process. Rather a Markovian dynamics can be (1976). [10] A.NakaharaandT.Isoda,Phys.Rev.E55,4264(1997). found [31]. These findings imply that the effective evo- [11] M. C. Teich, IEEE Trans. Biomed. Eng. 36, 150 (1989); lution of the ISI is similar to the one generated by our M. C. Teich, in Single Neuron Computation, ed. by T. model and caricatured by Eq. (8). McKenna,J.Davis,andS.F.Zornetzer(AcademicPress, San Diego, 1992), p. 589. .. −1 . [12] M. C. Teich, R. G. Turkott, and R. M. Siegel, IEEE 10 f −1.3 1100−01 1/τ1.1 Eng. Med. Biol. Mag. 15, 79 (1996); M. C. Teich, C. 10−2 τ)P(1100−−32 HOepnt.egShoacn.,ASm..BA. 1L4o,w5e2n9, (T1.99O7)z.aki, and E. Kaplan, J. 10−4 [13] M. Yamamoto et al.,Brain Research 366, 279 (1986). 10−2 10−1 100 101 102 [14] F. Gru¨neis et al., Biol. Cybern. 60, 161 (1989); F. S(f)10−3 τ/hours Gru¨neis, M. Nakao, and M. Yamamoto, ibid. 62, 407 (1990); F. Gru¨neis et al.,ibid. 68, 193 (1993). [15] T. Kodama et al., Brain Research 487, 26 (1989). 10−4 [16] R. Azouzand C. M. Gray, J. Neurosci. 19, 2209 (1999). [17] D.Heck,S.Rotter,andA.Aertsen,inBrainTheory,ed. by A. Aertsen (Elsevier, Amsterdam, 1993), p. 241. 10−5 [18] M.C.TeichandS.B.Lowen,IEEEEng.Med.Biol.Mag. 10−7 10−6 10−5 10−4 13, 197 (1994), and refs. therein. f/Hz [19] M.YamamotoandH.Nakahama,J.Neurophysiology49, 1182 (1983). FIG. 4. Power spectrum of the earthquake signal. Inset: [20] Theslopeofthevoltageincreaseisnotrelevantandisset ISI distribution. to1.WeassumeCl V0 Cu Cland0 D Cu Cl − ≪ − ≪ ≪ − to exclude trivial behavior. To summarize, we have presented a simple stochastic [21] H. G. Schuster, Deterministic Chaos, (VCH, Weinheim, model which is able to transform integrated white noise 1995), p.93. with a 1/f2 power spectrum into noise with 1/f tail. [22] J. Davidsen and H. G. Schuster, Phys. Rev. E 62, 6111 The basic mechanism is similar to what is expected to (2000). describe the dynamics of single neurons and, thus, our [23] B. Kaulakys and T. Meˇskauskas, Phys. Rev. E 58, 7013 modelcanexplainthebehaviorofneuronsinthecentral- (1998). nervous-system [14,19]. The analysis of earthquake data [24] The random walk with multiplicative noise can be ap- proximatedbyarandomwalkwithconstantvarianceand evensuggeststhattheeffectivedynamicsoftheISIisthe non-uniform waiting times between steps which are pro- same for many systems exhibiting 1/fα noise. −1 portional to τ . This explains the power-law decay of We thank A. Aertsen and C. Goltz for useful discus- k P(τ).Notethatdeviationsfromthisbehavioroccurclose sions. J. D. would like to thank the Land Schleswig- to the boundaries. Holstein, Germany, for financial support. [25] M. E. Wise, in Statistical Distributions in Scientific Work,ed.byC.Taillie,G.P.Patil,andB.A.Baldessari (D. Reidel, Hingham, MA, 1981), Vol. 6, p. 211. [26] W.R.SoftkyandC.Koch,NeuralComp. 4,643(1992), J. Neurosci. 13, 334 (1993). [27] M. Usher,M. Stemmler, and Z. Olami, Phys. Rev.Lett. 74, 326 (1995). ∗ Email address: [email protected] [28] S.Thurneretal.,Fractals5,565(1997),andrefs.therein. [1] W. H.Press, Comments Astrophys. 7, 103 (1978). [29] P. Reimann and P. H¨anggi, in Stochastic Dynamics, ed. [2] S. Kogan, Electronic noise and fluctuations in solids byL. Schimansky-Geier(Springer,Berlin, 1997), p.127; (Cambridge University Press, Cambridge, 1996). A. Fulin´ski, Phys.Rev.Lett. 79, 4926 (1997). [3] R. F. Voss and J. Clarke, Nature (London) 258, 317 [30] The investigated time signal contains all earthquakes (1975); J. Acoust. Soc. Am. 63, 258 (1978). above magnitude 1.8 in the area of 34◦0′ lat. N to [4] D. L. Gilden, T. Thornton, M. W. Mallon, Science 267, 35◦45′ lat.Nand115◦45′ lon.Wto118◦15′ lon.Wfrom 1837 (1995). 01/01/84 till12/31/00 whichcanbefoundintheSouth- [5] H. Yoshinaga, S. Miyazima, and S. Mitake, Physica A ern California Seismographic Network Catalogues. 280, 582 (2000). [31] J. Davidsen, C. Goltz, and H. G. Schuster, to be pub- [6] S. M. Bezrukov and M. Winterhalter, Phys. Rev. Lett. lished. 85, 202 (2000). [7] T. Field, U. Harder, and P. Harrison, cs.PF/0107001 (2001) and refs. herein. 4