Mean-field dynamics of a population of stochastic map neurons
Igor Franovic, Oleg V. Maslennikov, Iva Bacic, Vladimir I. Nekorkin
MMean-field dynamics of a population of stochastic map neurons
Igor Franovi´c, ∗ Oleg V. Maslennikov, † Iva Baˇci´c, and Vladimir I. Nekorkin Scientific Computing Laboratory, Center for the Study of Complex Systems,Institute of Physics Belgrade, University of Belgrade, Pregrevica 118, 11080 Belgrade, Serbia Institute of Applied Physics of the Russian Academy of Sciences,46 Ulyanov Street, 603950 Nizhny Novgorod, Russia (Dated: November 8, 2018)We analyze the emergent regimes and the stimulus-response relationship of a population of noisymap neurons by means of a mean-field model, derived within the framework of cumulant approachcomplemented by the Gaussian closure hypothesis. It is demonstrated that the mean-field modelcan qualitatively account for stability and bifurcations of the exact system, capturing all the genericforms of collective behavior, including macroscopic excitability, subthreshold oscillations, periodicor chaotic spiking and chaotic bursting dynamics. Apart from qualitative analogies, we find asubstantial quantitative agreement between the exact and the approximate system, as reflectedin matching of the parameter domains admitting the different dynamical regimes, as well as thecharacteristic properties of the associated time series. The effective model is further shown toreproduce with sufficient accuracy the phase response curves of the exact system and the assembly’sresponse to external stimulation of finite amplitude and duration.
Cortical connectivity patterns exhibit a hierarchicalmodular organization from the microscopic level of in-teracting neurons, via cortical columns and other typesof mesoscopic circuitry, up to fibers projecting betweenthe distributed brain areas [1–3]. Architecture of neuralassemblies comprising the anatomical modules closely re-flects the functional specialization over different modali-ties [4], such that the macroscopic dynamics of neuronalpopulations, as well as the interplay of the associatedcollective modes, underpin various stages of informationprocessing and higher cognitive functions [5]. In termsof evoked activity in the cortex, an ample example con-cerns the sensory regions, where many neurons displaysimilar responses to a given stimulus, which strongly in-dicates that the information content is primarily encodedby the collective assembly response [6]. With regard toself-organized dynamics, the different forms of synchro-nization between spiking or bursting neurons give rise tomacroscopic oscillations known to span several orders ofmagnitude in frequency range [7, 8]. Such oscillationsare deemed to facilitate coordinated rhythmic tasks andprovide the dynamical background behind perception, at-tention, memory consolidation, motor planning or sleep[7–9]. Nonetheless, interaction between multiple rhythmshas been indicated as critical for merging the activitiesof distant brain regions through cross-frequency couplingand other mechanisms [10, 11]. For these reasons, gain-ing a deeper understanding of macroscopic behavior ofneural populations has become an outstanding issue inneuroscience, focused primarily on the mechanisms guid-ing the emergent collective dynamics and the fashionin which assemblies respond to various external stimuli.Conceptually, such an approach is reminiscent of a com-mon paradigm in nonlinear dynamics, where systems of ∗ [email protected] † [email protected] coupled oscillators exhibiting a collective mode are typi-cally treated as macroscopic oscillators, which may thenbe subjected to an external drive or can be influenced bycollective rhythms from afferent populations [12].In the present paper, we systematically analyze theemergent dynamics and the stimulus-response relation-ship of a population of stochastic map neurons using amean-field ( M F ) approach. The considered map neuronscan exhibit a variety of regimes, including excitability,subthreshold oscillations, regular and chaotic spiking orbursting, as well as mixed spiking-bursting oscillations[13–16]. Despite that the collective motion of spiking orbursting neurons subjected to noise has been extensivelystudied using different models of discrete local dynamics,such as Rulkov [17–24] or Izhikevich neuron maps [25, 26],a
M F theory for a population of stochastic map neuronsis obtained here for the first time. Nevertheless, in caseof continuous time systems, the
M F approach has beena standard analytical tool for treating diverse problemsin neuroscience and other fields [27–32]. Note that ourderivation of the effective model relies on Gaussian ap-proximation, which is introduced within the frameworkof a Gaussian closure hypothesis [33–39].The particular set of issues we address consists in es-tablishing whether and how the
M F model can be usedto ( i ) qualitatively analyze the network stability and bi-furcations of the exact system associated to emergence ofgeneric macroscopic regimes; ( ii ) provide adequate quan-titative predictions in terms of bifurcation thresholds,and the average interspike intervals or bursting cyclesof the exact system, as well as ( iii ) accurately anticipatethe population’s response to different forms of externalstimuli. Within this context, it will be examined whetherthe effective model is capable of reproducing the proper-ties of noise-activated, noise-induced and noise-perturbedmodes of collective behavior.The paper is organized as follows. In Section I, wemake an overview of the local map dynamics and in- a r X i v : . [ n li n . C D ] M a r troduce the population model. Section II outlines theingredients most relevant for the derivation of the M F system, with the remaining technical details left for theAppendix. In Section III, the qualitative and quantita-tive agreement between the dynamics of the exact andthe
M F model is illustrated by the appropriate bifurca-tion diagrams, as well as by comparing the characteris-tic features of the associated regimes. Section IV con-cerns the assembly’s stimulus-response relationship, firstinvestigating the analogy between the respective phase-response curves (
P RCs ) of the exact system and the ef-fective model in spiking and bursting regimes, and thenconsidering the extent to which the
M F model repro-duces the population’s response to rectangular pulses offinite amplitude and duration. In Section V, we providea summary of our main results.
I. MAP NEURON DYNAMICS AND THEPOPULATION MODEL
The dynamics of an isolated neuron conforms to a mapmodel first introduced in [40, 41], which is given by x n +1 = x n + G ( x n ) − βH ( x n − d ) − y n , (1) y n +1 = y n + (cid:15) ( x n − J ) , where n denotes the iteration step. The variable x n qual-itatively accounts for the membrane potential, whereasthe recovery variable y n , whose rate of change is set bya small parameter (cid:15) = 10 − , mimics the behavior of ion-gating channels. The parameters a, β and d modify theprofile of the ensuing oscillations, while J crucially in-fluences the neural excitability, viz. the transitions fromsilence to active regimes.The x n evolution features two nonlinear terms, one be-ing a FitzHugh-Nagumo-like cubic nonlinearity G ( x n ) = x n ( x n − a )(1 − x n ), which is complemented by a disconti-nuity term − βH ( x n − d ), where H stands for the Heavi-side step function. The parameters a = 0 . d = 0 . (cid:15) = 0) a Lorenz-type map within certain parameterdomains [41, 42], which endows the model with the abilityto generate chaotic spike or burst oscillations, otherwiselacking in the Fitzhugh-Nagumo type of systems.Under variation of J and β , the map (1) may repro-duce a rich repertoire of generic regimes displayed by thereal neurons, as demonstrated in Fig. 1. In particular,the main frame shows amplitudes of the corresponding x time series for the given ( J, β ), while the remaining sub-figures illustrate the characteristic waveforms pertainingto excitable regime (region I ), subthreshold oscillations( II ), regular ( III ) or chaotic spiking ( IV ), chaotic burst-ing ( V ), as well as the mixed chaotic spike-burst activity( V I ). Some of the indicated boundaries, such as thoseinvolving domains
IV, V and
V I should be understoodas tentative, since the associated transitions are smoothand therefore difficult to discern. The detailed phase plane analysis concerning the rel-evant unstable invariant curves and the mechanismsunderlying transitions between the different dynamicalregimes can be found in [43]. Here we briefly mentionthat under increasing J , the equilibrium loses stabilityvia the Neimarck-Sacker bifurcation, which gives rise tosubthreshold oscillations. Note that the latter may beconsidered as an excitable state, in a sense that a strongenough perturbation can elicit genuine spike, though thephase point does not relax to the equilibrium, but ratherto a closed invariant curve.Adopting model (1) for local dynamics, we focus ona population of N stochastic neurons coupled in the all-to-all fashion via electrical synapses (diffusive couplings).Each neuron receives input from the units within the as-sembly, and is further influenced by synaptic noise fromthe embedding environment. The population activity isthen described by the following system x i,n +1 = x i,n + G ( x i,n ) − βH ( x i,n − d ) − y i,n + I syni,n , (2) y i,n +1 = y i,n + (cid:15) ( x i,n − J ) ,I syni,n = I coupi,n + I randi,n = cN N (cid:88) j =1 ,j (cid:54) = i ( x j,n − x i,n ) + σξ i,n , where i specifies the particular neuron. The synapticcurrents I syni,n comprise two types of terms. The diffusivecouplings I coupi,n are characterized by the strength c , whichis assumed to be uniform over the network and is setto c = 1 in the remainder of the paper. The randominputs I randi,n involve uncorrelated white noise ( E [ ξ i,n ] =0 , E [ ξ i,n ξ j,n (cid:48) ] = δ ij δ ( n − n (cid:48) )) of intensity σ .Confined to a single unit, the stochastic componentmay influence its dynamics either by perturbing the de-terministic oscillatory regimes, or by inducing oscillationsin the excitable regime, cf. Fig. 2(b). The onset of noise-induced spiking or bursting within the parameter domainwhere the fixed point is deterministically stable (domain I in Fig. 1) corresponds to a phenomenon of stochasticbifurcation [38, 44–47]. The latter are typically describedphenomenologically, in a sense that certain time-averagedquantities, such as the asymptotic probability distribu-tions of relevant variables or the associated power spec-tra, exhibit a qualitative change under variation of noiseintensity. For instance, in continuous-time systems, it hasbeen shown that the stochastic Hopf bifurcation from astochastically stable fixed point to a stochastically stablelimit cycle is accompanied by the loss of Gaussian prop-erty for the asymptotic distributions of the appropriatevariables [48]. At variance with standard deterministicbifurcations, where one clearly observes a critical valueof the control parameter, the change of systems behaviorin noise-induced transitions is gradual [38]. Note thatnoise can also play an important part in the ( J, β ) region II where the deterministic map shows subthreshold os-cillations. Here noise can give rise to a form of dynamicsreminiscent of mixed-mode oscillations, cf. Fig. 2(c). J A x n x n x n x n x n I II IIIIV V VII II IIIIVV VI VI x n FIG. 1. (Color online) Dynamical regimes exhibited by model (1). The heat map refers to variation of the amplitude ofoscillations A of the x time series in the J − β plane. The waveforms shown in subfigures I − V I illustrate the different formsof neuron’s behavior, including excitability ( I ), subthreshold oscillations ( II ), regular spiking ( III ), chaotic bursting ( IV ),chaotic spiking ( V ), as well as the mixed spike-burst activity ( V I ). The dots in the heat map indicate the particular (
J, β )values where the representative waveforms are obtained. -0.2 0.2 0.4 0.60-0.020.020.040.060.080 y n x n x = d x = J x n x n (a) (b)(c) y=G(x) FIG. 2. (Color online) Impact of noise on a single map neu-ron in the excitable regime. (a) indicates the mechanismbehind noise-induced spiking. The data are obtained for J = 0 . , β = 0 . , σ = 0 . x = J intersects theinvariant curve y = G ( x ) below the curve’s minimum. (b)shows the x n series corresponding to noise-induced bursting( J = 0 . , β = 0 . , σ = 0 . J = 0 . , β = 0 . , σ = 0 . So far, models similar to (2) have been applied to ad-dress a number of problems associated to collective phe-nomena in networks of coupled neurons, including syn-chronization of electrically coupled units with spike-burstactivity [49, 50], pattern formation in complex networkswith modular architecture [13, 14, 51], transient clusteractivity in evolving dynamical networks [16], as well asthe basin stability of synchronization regimes in small-world networks [15]. Within this paper, the collectivemotion will be described in terms of the global variables X n = N (cid:80) Ni =1 x i,n and Y n = N (cid:80) Ni =1 y i,n . II. DERIVATION OF THE MEAN-FIELDMODEL
Considering a
M F approximation, our main goal liesin deriving a reduced low-dimensional deterministic set of nonlinear difference equations whose dynamics is qual-itatively analogous to the collective motion of the originalsystem (2) comprised of 2 N coupled stochastic maps. Inparticular, the M F model should be able to generateall the regimes exhibited by the exact system, qualita-tively reproducing the bifurcations that the latter under-goes. Also, applying the effective model, one should becapable of inferring with sufficient accuracy the parame-ter domains which admit the different collective states ofthe exact system, with the corresponding time series ex-hibiting similar characteristic quantitative features. Re-garding the explicit effects of noise, the
M F model isexpected to account for the onset or suppression of dif-ferent types of collective modes associated to macroscopicspiking or bursting activity, which are mediated by syn-chronization or desynchronization of individual neurondynamics, respectively. The synchronization processesmay be influenced by noise in a variety of ways, includ-ing the scenarios where noise acts as a perturbation tomainly deterministic (and chaotic) local oscillations, orthe ones where noise plays a facilitatory role, in a sensethat the collective mode emerges via synchronization ofnoise-induced local dynamics.Given that we consider a system of discrete-time equa-tions, one cannot adopt the usual method of derivingthe
M F model via Fokker-Planck formalism [39]. Nev-ertheless, an analytically tractable
M F model may stillbe built by focusing on the evolution of cumulants [33–35, 38], whereby the full density of states is factorized intoa series of marginal densities. The advantage of such anapproach is that the simplifying approximations aimedat truncating the underlying cumulant series can be in-troduced in a controlled fashion. Such approximations,stated in a form of closure hypothesis [33], are requireddue to nonlinearity of the original system, which causesthe dynamics of cumulants of the given order to be cou-pled to those of the higher order.In our case, the derivation of the effective model in-corporates an explicit Gaussian closure hypothesis [33–35, 38], by which all the cumulants above second orderare assumed to vanish. The collective dynamics is thendescribed by a set of five variables (the first- and second-order cumulants), including(i) the means, given by m x,n = lim N →∞ N (cid:80) Ni =1 x i,n ≡(cid:104) x i,n (cid:105) , m y,n = lim N →∞ N (cid:80) Ni =1 y i,n ≡ (cid:104) y i,n (cid:105) ;(ii) the variances, defined as S x,n = (cid:104) x i,n (cid:105) − (cid:104) x i,n (cid:105) = (cid:104) x i,n (cid:105) − m x,n and S y,n = (cid:104) y i,n (cid:105) − (cid:104) y i,n (cid:105) = (cid:104) y i,n (cid:105) − m y,n ;(iii) the covariance U n = (cid:104) x i,n y i,n (cid:105) − m x,n m y,n . The expressions for higher order moments (cid:104) x ki,n (cid:105) in termsof the first- and second-order cumulants [52], such as (cid:104) x i (cid:105) = m x + 3 m x S x (3) (cid:104) x i (cid:105) = m x + 6 m x S x + 3 S x (cid:104) x i y i (cid:105) = m y S x + m y m x + 2 m x U (cid:104) x i y i (cid:105) = 3 S x U + 3 S x m x m y + 3 m x U + m y m x (cid:104) x i (cid:105) = m x + 15 m x S x + 10 m x S x (cid:104) x i (cid:105) = m x + 15 S x + 15 m x S x + 45 m x S x , can be derived using the closure hypothesis.The Gaussian approximation effectively amounts to anassumption that the relationlim N →∞ N N (cid:88) i =1 x ki,n ≈ E [ x ki,n ] , (4)holds, whereby E refers to expectation value obtained byaveraging over an ensemble of different stochastic realiza-tions. In other words, one supposes that the local vari-ables are independent and are drawn from a normal dis-tribution N ( m x , S x ). We do not know a priori whethersuch an assumption is fulfilled or not, but can only judgeon its validity by verifying the correctness of the predic-tions on the population dynamics provided by the M F model. Also note that the effective model concerns theassembly dynamics in the thermodynamic limit N → ∞ .The stochastic terms in this case can be neglected, as onemay show them to contribute to finite size effects whichscale as 1 /N . This means that the influence of noise inour M F model is felt only via the noise intensity, whichassumes the role of an additional bifurcation parameter.Let us now illustrate the main technical points requiredfor the derivation of the
M F model. Our focus will liewith a couple of relevant examples, whereas the remain-ing details are provided in the Appendix. We begin byconsidering the dynamics of the global variable m x , whichis given by m x,n +1 = m x,n − m y,n + (cid:104) G ( x i,n ) (cid:105) − β (cid:104) H ( x j,n − d ) (cid:105) (5) It is easy to see that there is no contribution from thecoupling term. As far as the third term on the r.h.s. ofEq. (5) is concerned, using Eq. (3), one arrives at (cid:104) G ( x i ) (cid:105) = (cid:104)− x i +(1+ a ) x i − ax i (cid:105) = G ( m x )+ S x (1+ a − m x ) . (6)In the last expression, we have dropped the time index forsimplicity and have introduced the shorthand notation G ( m x ) ≡ − m x + (1 + a )( m x + S x ).The key problem is how to treat the final term inthe r.h.s. of Eq. (5). Our approach consists in re-placing the assembly average by the expectation value( (cid:104) H ( x i − d ) (cid:105) ≈ E [ H ( x i − d )]), obtained by assuming thatthe local variables at an arbitrary time moment are nor-mally distributed according to P ( x i ) ∼ N ( m x , S x ). Theexpectation may then be evaluated as E [ − β (cid:104) H ( x i − d ) (cid:105) ] = (cid:90) dx (cid:90) dx ... (cid:90) dx N × (7)( − βN (cid:88) i H ( x i − d )) p ( x , x , ..., x N ) = − β (cid:90) ∞−∞ dx H ( x − d ) p ( x ) = − β (cid:90) ∞ d √ πS x e − ( x − mx )22 Sx = − β − Erf (cid:20) d − m x √ S x (cid:21) ) , with the error function Erf ( x ) = √ π (cid:82) x e − t dt . In theabove calculation, we have explicitly used the assumptionon the independence of distributions of local variables atany given moment of time.In a similar fashion, one may consider the S x dynam-ics, which constitutes the most demanding part of thederivation. In particular, proceeding from the S x defini-tion, we obtain S x,n +1 = (cid:104) x i,n +1 (cid:105) − (cid:104) x i,n +1 (cid:105) = (cid:104) [(1 − c ) x i,n + G ( x i,n )(8) − βH ( x i,n − d ) − y i,n + ξ i,n + cm x,n ] (cid:105)− ( m x,n − m y,n + G ( m x,n ) + S x,n (1 + a − m x,n ) − β (cid:104) H ( x i,n − d ) (cid:105) ) . As an illustration, let us evaluate one of the terms con-taining an average over the threshold function: − βE [ (cid:104) G ( x i ) H ( x i − d ) (cid:105)− (cid:104) G ( x i ) (cid:105)(cid:104) H ( x i − d ) (cid:105) ] = (9) − β (cid:20)(cid:90) dx G ( x ) H ( x − d ) p ( x ) − (cid:90) dx H ( x − d ) p ( x ) [ G ( m x ) + S x (1 + a − m x )] (cid:21) ≈ − β (cid:20)(cid:90) dx ( G ( m x ) + G (cid:48) ( m x )( x − m x ) + 12 G (cid:48)(cid:48) ( m x ) × ( x − m x ) ) H ( x − d ) p ( x ) − (cid:90) dx H ( x − d ) p ( x ) × [ G ( m x ) + S x (1 + a − m x )]] = ... = − β [(1 + a )( m x + d ) − a − m x d ] (cid:114) S x π exp (cid:20) − ( d − m x ) S x (cid:21) . Again, the time indexes have been suppressed to simplifythe notation.Leaving the remaining elements of the derivation forthe Appendix, we now state the final equations of the
M F model in the thermodynamic limit m x,n +1 = m x,n − m y,n + G ( m x,n ) + S x,n (1 + a − m x,n )(10) − β − Erf (cid:34) d − m x,n (cid:112) S x,n (cid:35) ) m y,n +1 = m y,n + (cid:15) ( m x,n − J ) S x,n +1 = (1 − c ) S x,n + S y,n + σ − − c ) U n + S x,n ( − m x,n + 2(1 + a ) m x,n − a ) − − c ) × (3 m x,n S x,n + 3 S x,n − a ) m x,n S x,n + aS x,n )+ 2(3 S x,n U n + 3 m x,n U n − a ) m x,n U n ) − β [(1 + a )( m x,n + d ) − a − dm x,n ] (cid:114) S x,n π × exp (cid:20) − ( d − m x,n ) S x,n (cid:21) − β (1 − c ) (cid:114) S x,n π × exp (cid:20) − ( d − m x,n ) S x,n (cid:21) + S x,n (cid:2) m x,n − a ) × m x,n + 2(1 + a ) + 6 a (cid:3) + 15 S x,n S y,n +1 = S y,n + (cid:15) S x,n + 2 (cid:15)U n U n +1 = U n − ( a + c + (cid:15) ) U n + (cid:15) (1 − c − a ) S x,n − S y,n − ( U n + (cid:15)S x,n )(3 S x,n + 3 m x,n − a ) m x,n ) − β(cid:15) (cid:114) S x,n π exp (cid:20) − ( d − m x,n ) S x,n (cid:21) . III. ANALYSIS OF STABILITY ANDBIFURCATIONS
In this section, our goal is to demonstrate the quali-tative and quantitative analogies between the dynamicsof the exact system and the
M F model. To this end,we first examine the succession of macroscopic regimesin the J − β parameter plane for σ fixed at an interme-diate value σ = 0 . J is relevant for the systems excitability,viz. the transitions from silent to active regimes, while β influences the waveforms of the active states (spiking,bursting, or mixed spike-bursting activity). The assem-bly is found to exhibit the collective modes which qual-itatively correspond to the dynamics of a single unit il-lustrated in plates III − V I of Fig. 1. The heat maps inthe left column of Fig. 3 provide a comparison betweenthe oscillation amplitudes A of the global variable X (toprow) and the M F variable m x (bottom row) for the given( J, β ). The right column indicates how well are matchedthe average interspike interval (or the average burstingcycle) T of the exact system with the corresponding char- J A 0.05 0.10 0.1500.10.20.30.40.50 J T0.05 0.10 0.1500.10.20.30.40.50 J A 0.05 0.10 0.1500.10.20.30.40.50 J T (a) (b)(c) (d) FIG. 3. (Color online) Heat maps in (a) and (b) show thedependencies A ( J, β ) and T ( J, β ) obtained by stochastic av-eraging for a network of N = 100 neurons, respectively. (c)and (d) illustrate the analogous results for the MF model.The noise intensity in all instances is σ = 0 . acteristics of the dynamics of the M F model (11). In thegiven instances, exact system comprises an assembly of N = 100 neurons, having obtained A by averaging over asufficiently long time series, whereas T is determined bytaking average over an ensemble of 20 different stochasticrealizations. With regard to T , we have selected a con-venient threshold θ = 0 .
2, which allows a clear detectionof individual spikes, and also enables one to unambigu-ously discern the initiation stage of bursts, as requiredfor calculating the length of the bursting cycle.Let us begin the analysis by focusing on the domainof J values where the exact system exhibits the stochas-tically stable equilibrium, while the M F model has astable stationary state. The stochastic stability physi-cally implies that fluctuations around the deterministicfixed point are typically of the order of noise, thoughsome rare spikes may still be evoked. For J sufficientlyclose to the region admitting the subthreshold oscilla-tions, the population manifests macroscopic excitability.To properly illustrate this feature, we have analyzed theassembly dynamics in the limit σ = 0, cf. Fig. 4. In par-ticular, figures 4(a) and 4(b) show the maximum X and m x values reached in the corresponding time series ob-tained for sets of different initial conditions ( X , Y ) and( m x, , m y, ), respectively. The comparison between thetwo plots clearly corroborates that the boundary definingthe domain of spiking response is appropriately antici-pated by the M F model. An important remark is thatfor the given J , the assembly may exhibit different formsof macroscopic excitability, generating a single spike ora burst of spikes, as dependent on the value of β . Thisis demonstrated by the time series in figures 4(c) and4(d). The former refers to a one-spike response in case -0.4 -0.2 0.2 0.40-0.3-0.2-0.10.10.20.30 X Y A -0.4 -0.2 0.2 0.40-0.3-0.2-0.10.10.20.30 m x, m y , A X , m x n X m x X , m x n Xm x (a) (b)(c) (d) FIG. 4. (Color online) Macroscopic excitability feature. In (a)and (b) are shown the maximum values of X and m x reachedwithin the time series of the exact and the MF system,starting from the analogous initial conditions ( X , Y ) and( m x, , m y, ), respectively. The parameters are J = 0 . , β =0 .
4. (c) illustrates the case where a strong enough pertur-bation elicits a single-spike response ( J = 0 . , β = 0 . J = 0 . , β = 0 . MF model (dotted line) is indistinguishable fromthat of the exact system (dashed line). of β = 0 .
4. For smaller β , one observes responses com-prising two or more closely packed spikes, with Fig. 4(d)illustrating a three-spike burst encountered for β = 0 . M F model are exactly matched in the limit σ = 0.Next we address the noise-influenced transitions fromsilence to active regimes observed under increasing J . Todo so, in Figure 5(a) we have plotted the change of thefiring (spiking or bursting) frequency R for an assemblyconsisting of N = 100 neurons. The average frequencyis determined by considering an ensemble of 20 differentstochastic realizations, having σ fixed to the moderatevalue from Fig. 4. The results from simulations of thefull system (2) are compared against the data obtainedfor the M F model. In this context, two points should bestressed. First, for moderate σ , note that the firing fre-quencies of the M F model lie in close agreement to thoseof the exact system. As a second point, one finds thatsuch quantitative agreement extends to different formsof collective behavior, viz. it holds for different types oftransitions from silent to active regimes. As already indi-cated, the waveforms pertaining the active states dependon β , such that the associated transitions are mediatedby the distinct synchronization processes. For instance,at β = 0, synchronization involves time series of singleunits that conform to spiking activity of type III fromFigure 1, which are quite resilient to impact of noise. Onthe other hand, for β = 0 . β = 0 .
4, the individualunits exhibit chaotic bursting or spiking activity, respec- tively, such that the underlying synchronization processmay be more susceptible to stochastic effects. The typical X time series illustrating the different collective modesare compared to the corresponding m x series in figures5(b)-(e). The top (bottom) row concerns the data for theexact system ( M F model).In order to investigate more closely the influence ofnoise for J interval in vicinity of the transition fromsilence to active regimes, we examine how the profilesof R ( J ) curves change under increasing σ . The resultsshown in Fig. 6 refer to β = 0 . N = 100 neurons. As expected, the transitionappears quite sharp for moderate noise σ = 0 . σ , e. g. σ = 0 . M F model for σ = 0 . σ , the M F model fails to reproduce the be-havior of the exact system in vicinity of threshold J , ina sense that it overestimates the maximal R value, aswell as the actual critical J characterizing the transition.Viewed from another angle, one may infer that for suf-ficiently large σ and J below the threshold given by the M F model, the latter fails to capture the impact of syn-chronization processes taking place between the noise-induced oscillations of individual units. This especiallyrefers to J interval where the spikes or bursts (dependingon the given β ) are superimposed on the background ofsubthreshold oscillations. An example of such a discrep-ancy between the behavior of the exact and the effectivesystem is provided in Fig. 7, cf. Fig. 7(a) and Fig.7(c). Also, for strong σ and J values above the tran-sition, the firing frequencies anticipated by the effectivemodel are typically higher than those of the exact system(not shown). Within this region, the stochastic effectssuppress synchronization between the chaotic oscillationsof single neurons, thereby reducing the corresponding R value. This is not accounted for with sufficient accu-racy by the M F system. Note that such suppression ofsynchronization is reflected in the corresponding X se-ries by the spike (burst) ”skipping” mechanism, wherethe large-amplitude oscillations are occasionally replacedwith subthreshold oscillations. For the associated J and σ values, such a phenomenon is absent in the dynamics ofthe effective model, cf. Fig. 7(b) and Fig. 7(d). In bothof the scenarios illustrated in Fig. 7, the reason for hav-ing the M F model fail lies in that the Gaussian approx-imation behind it breaks down due to large stochasticfluctuations.The fashion in which the validity of the effective modelspredictions deteriorates with increasing σ is made moreexplicit in Fig. 8, which shows the A ( J, σ ) and T ( J, σ )dependencies for the exact and the approximate systemat fixed β = 0 .
4. The considered size of the network is N = 100. Comparison between the respective A (leftcolumn) and T plots (right column) suggests that therange of σ values where the M F approximation appliesis contingent on J . For instance, in the J region be-low the deterministic threshold, one may estimate this R J b = ex b = mf b = ex b = mf b = ex b = mf b = , ex b = , mf X t m x n X t m x n (a) (b) ( c )( d ) (e) FIG. 5. (Color online) (a) shows a family of R ( J ) curves over β for a network of size N = 100 under fixed σ = 0 . MF model, whereby the symbols × , + , ∗ , (cid:63) correspond to cases β = 0 , . , . . X series associated to the spiking and the bursting collective modes. The considerednetwork is made up of N = 100 neurons, with the parameters set to J = 0 . , β = 0 . , σ = 0 .
001 in (b), and J = 0 . , β =0 . , σ = 0 .
001 in (c). In (d) and (e) are provided the m x series obtained for parameters from (b) and (c). range by noting that the effective bifurcation diagram inFig. 8(a) indicates that noise-induced macroscopic os-cillations emerge for σ ≈ . σ ≈ .
003 within the given J re-gion. Nevertheless, for J above the deterministic thresh-old, the validity of the M F model appears to dependrather strongly on particular J , with the σ values wherethe Gaussian approximation effectively fails spanning therange σ ∈ (0 . , . N = 100to those obtained for the effective system. Nevertheless,within Section II, it has already been emphasized thatthe M F model, deterministic in character, refers to thesystems behavior in the thermodynamic limit N → ∞ ,whereas the explicitly stochastic terms could only be in- R J s = 0 . 0 0 1 s = 0 . 0 1 s = 0 . 0 2 s = 0 . 0 5 m f , s = 0 . 0 0 1 FIG. 6. (Color online) Family of R ( J ) curves over σ obtainedfor a network of N = 100 neurons under fixed β = 0 .
2. Thedifferent symbols correspond to cases σ = 0 .
001 (squares), σ = 0 .
01 (circles), σ = 0 .
02 (triangles) and σ = 0 .
05 (dia-monds). The crosses connected by the dashed line highlightthe R ( J ) curve for the MF model at σ = 0 . corporated as finite-size effects. This makes it relevantto examine how the behavior of the exact system withinthe J domain around deterministic threshold changes forlarge and fixed σ under increasing N . To this end, wehave plotted in Fig. 9 the R ( J ) curves calculated for N = 100 (squares), N = 500 (circles) and N = 1500(diamonds) at fixed β = 0 . , σ = 0 .
05. The curve for N = 100 evinces that the given σ value is quite large ina sense of being sufficient to induce collective oscillationswithin the excitable regime. Apart from the dependen-cies for the full system, we also show the R ( J ) curve as-sociated to the M F model (dashed line with crosses). Aninteresting point regarding the latter is that the J thresh-old for the emergence of the collective mode is shifted to-ward a larger value compared to the case σ ≈ .
01. Whilethe given transition itself appears quite sharp, the curvescorresponding to the exact system approach it with in-creasing N , both in terms of the J threshold and the R values above the transition. This corroborates that the( J, σ ) domain where the Gaussian approximation behindthe
M F model fails expectedly reduces with the increas-ing system size. X n m x n X n m x n ( a ) (b)(c) (d) FIG. 7. Noise-induced phenomena within the J interval invicinity of the deterministic threshold. X series in (a) showsthe noise-induced spike-bursting activity on top of subthresh-old oscillations ( J = 0 . , β = 0 . , σ = 0 . X variable ( J = 0 . , β = 0 . , σ = 0 . m x series corresponding to parameter sets from(a) and (b), respectively. J A J T J A J T (a) (b)(c) (d) FIG. 8. (Color online) (a) and (b) respectively refer to A ( J, σ )and T ( J, σ ) dependencies for the network of N = 100 neuronsunder fixed β = 0 .
4. The results in (a) are obtained by aver-aging over a sufficiently long time series, whereas data in (b)derive from averaging over an ensemble of 20 different stochas-tic realizations. In (c) and (d) are provided the A ( J, σ ) and T ( J, σ ) dependencies determined by numerical simulations ofthe MF model. IV. RESPONSE TO EXTERNAL STIMULI
The aim of this section is to investigate the extent towhich the
M F model can be used to predict the stimulus-response relationship of an assembly exhibiting differentmacroscopic regimes, including the excitable state, aswell as the spiking and bursting collective modes. Let R J N = 1 0 0 N = 5 0 0 N = 1 5 0 0 m f
FIG. 9. (Color online) R ( J ) dependencies for increasing N under fixed ( β, σ ) = (0 . , . N = 100, N = 500 and N = 1500,respectively. The results predicted by the MF model are in-dicated by crosses connected via dashed line. us first focus on the two latter instances and examinethe sensitivity of a population to an external pulse per-turbation within the framework of phase resetting theory[53–56]. In order to compare the behavior of the exactsystem and the effective model, we determine the corre-sponding phase resetting curves ( P RC s), which describethe phase shift ∆ φ , induced by the perturbation, in termsof the phase φ p when the perturbation is applied. Theconsidered stimulus has a form of a short pulse current I p = a p H ( n − n i ) H ( n − n f ), whose magnitude a p andwidth ∆ = n i − n f are small compared to the amplitudeand duration of the spiking (or bursting) cycle T , re-spectively. In case of the exact system, the same pulsecurrent is delivered to each neuron i , adding the term I p to x i dynamics, whereas in the effective model, stimula-tion is administered via the m x variable. The phase ϕ p is defined in reference to T by ϕ p = n p /T . The associ-ated phase difference following the reset is calculated as∆ ϕ = 1 − T /T , where T denotes the duration of theperturbed spiking or bursting cycle.The P RC s characterizing the assembly response in thespiking regime are provided in Fig. 10(a) and Fig. 10(b),whereby the former is obtained under the action of anexcitatory ( a p > a p < ∆ φ exmf X ∆ φ exmf ∆ φ exmf X ∆ φ exmf φ p φ p φ p φ p (a) (b)(c) (d) φφ FIG. 10. (Color online) Assembly phase resetting. (a)and (b) show the
P RC s for a population in spiking regime( J = 0 . , β = 0) under excitatory ( a = 0 . a = − . N = 500) are indicated by the solid line,whereas the data for the MF model are denoted by circles.The bottom row illustrates the P RC s for an assembly ex-hibiting macroscopic bursting ( J = 0 . , β = 0 . a = 0 . a = − . π . ex X n m x n mf I p n X n ex m x n mf I p n X n ex m x n mf I p n (a)(b)(c) (d)(e)(f) (g)(h)(i) FIG. 11. Stimulus-response relationship in the excitable regime ( J = 0 . MF model), whereas the bottom row shows the profile of the external stimulation. In the panels (a)-(c), thesystem parameters are β = 0 . , σ = 0, while the perturbation is characterized by a p = 0 . , ∆ = 200. Panels (d)-(f) concern theresponse of an assembly ( β = 0 . , σ = 0 . a p = 0 . , ∆ = 200. Panels (g)-(i) illustrate theresponse of a population ( β = 0 . , σ = 0 . a p = 0 . , ∆ = 50. The considered networkis of size N = 500. in both instances, the results derived from the effectivemodel, denoted by circles, show excellent agreement withthe data for the exact system (solid lines). In qualitativeterms, one observes that an excitatory stimulation mayadvance the phase of the spiking cycle if it arrives suffi-ciently close to the spike, but still before the sharp risingstage. Nevertheless, an excitatory perturbation whichacts during the spike or within the effective refractoryperiod has a suppression effect, reflected in delaying ofthe next spike. At variance with the excitatory stimula-tion, the inhibitory pulse postpones the next firing timeif it is introduced within the interval close to the risingstage of spike.The P RC s determined for an assembly exhibiting col-lective bursting show qualitatively analogous effects tothose described so far, see Fig. 10(c) and Fig. 10(d).This especially refers to impact of perturbation deliv-ered sufficiently close to a moment of burst initiation.An apparent difference compared to Fig. 10(a) and Fig.10(b) emerges during the bursting stage itself, where theassociated
P RC s expectedly exhibit strong fluctuations.Apart from that, one finds an interesting effect that boththe excitatory and the inhibitory stimulation have a fa-cilitatory role, i. e. cause phase advancement during therelaxation stage of the bursting cycle.For a population in the excitable state, we considerscenarios where the system is influenced by a rectangularpulse perturbation of finite magnitude and duration, in asense that the latter are comparable to corresponding fea-tures of typical spiking or bursting cycles. Note that theselected J value J = 0 .
02 lies sufficiently away from theinterval admitting the subthreshold oscillations. Again,our objective is to determine whether the
M F model cor-rectly anticipates the response of the exact system, nowin presence of small to moderate noise. Some of the il-lustrative examples concerning the stimulus-response re-lationship under the finite perturbation are provided inFig. 11. The top and the middle row refer to X and cor-responding m x time series, respectively, while the bot-tom row shows the profile of the applied stimulus. We find that in the absence of noise or for sufficiently small σ , the effective model reproduces the evoked behavior ofthe full system quite accurately. This also refers to somehighly complex forms of responses, as corroborated inFig. 11(a)-(c), which concern relatively large a p and ∆.Under increasing σ , the ability of the M F model to pre-dict the dynamics of the exact system gradually reduces,but in a fashion that involves a nontrivial dependenceon β . In particular, for smaller β ≈ .
1, which wouldfacilitate macroscopic spiking mode for supercritical J ,it turns out that the dynamics of the M F model lies inclose agreement to the one of the exact system even formoderate noise σ = 0 . β , such an analogy between the responses ofthe exact and the M F system is lost, see Fig. 11(g)-(i). Naturally, the validity of the predictions given bythe
M F model deteriorates if the stimulation amplitude a p and the duration ∆ are large, especially in presenceof non-negligible noise. V. SUMMARY AND DISCUSSION
We have developed a
M F approach in order to system-atically analyze the emergent dynamics and the input-output relationship of a population of stochastic mapneurons. The reduced low-dimensional model has beenderived within the framework of Gaussian approxima-tion, formally introduced in a form of a closure hypoth-esis. In physical terms, such an approximation suggeststhat the local variables at an arbitrary moment of timeare independent and conform to a normal distributioncentered about the assembly mean and characterized bythe associated assembly variance. Validity of such an ap-proximation cannot be established a priori , but has beensystematically verified by numerically corroborating thatthe
M F model reproduces the behavior of the exact sys-tem with sufficient accuracy.In particular, we have first demonstrated that the ef-fective model can qualitatively capture all the bifurca-0tions of the exact system leading to the onset of differentgeneric regimes of collective behavior. As far as the quan-titative agreement is concerned, we have established sub-stantial matching between the parameter domains admit-ting the respective dynamical regimes for the exact andthe approximate system. Moreover, the typical featuresof the associated regimes, such as the average interspikeinterval or the average bursting cycle, exhibit analogouschanges with parameter variation, and in many parame-ter domains display numerically similar values.An important issue has been to explicitly examine howthe effects of noise are reflected in the behavior of the
M F model. For the noise-perturbed activity, where thesufficiently small noise weakly influences the determinis-tic attractors of the system, the obtained results indicatethat the Gaussian approximation holds. Nevertheless,the physical picture changes in case of noise-induced col-lective behavior. In particular, for different scenarios ofstochastic bifurcations, typically corresponding to transi-tions from subthreshold oscillations, which involve gener-alized excitability feature, to spiking or bursting regimes,the exact system undergoes a gradual (smooth) changeof collective dynamics, whereas the
M F model exhibits astandard deterministic bifurcation with a sharp bifurca-tion threshold. In such instances, the collective variablesof exact system manifest large fluctuations, which explic-itly violate the Gaussian approximation behind the effec-tive model. Note that the loss of Gaussianity propertyfor asymptotic distribution of relevant variables, whichaccompanies the described stochastic bifurcations, doesnot imply per se that our Gaussian approximation fails inthe supercritical state. This point is evinced by the factthat the dynamics of the effective model shows qualita-tively and quantitatively similar features to those of theexact system if the considered parameters lie sufficientlyabove the stochastic bifurcation. In fact, the Gaussianapproximation applied in the derivation of the
M F modelbreaks down only in vicinity of such transitions, wherethe finite-size effects neglected in Eq. (11) become mostprominent. We have numerically verified the prevalenceof finite-size effects in these parameter domains, show-ing that the change of the appropriate order parameter,such as the spiking frequency, becomes sharper as thesize of the neural assembly is increased. Nevertheless,the validity of Gaussian approximation is regained oncethe system is sufficiently above the bifurcation.Apart from considering asymptotic dynamics, we haveverified that the
M F model is capable of capturing thestimulus-response features of the exact system. For shortpulse-like perturbations, it has been found that the ap-proximate system reproduces the
P RC s of the exact sys-tem for both the spiking and bursting regimes of collec-tive activity with high accuracy. Substantial analogieshave also been observed in case of macroscopic excitableregime for scenarios where the assembly is stimulated byrectangular pulse perturbations of finite amplitude andduration.Having developed a viable
M F approach, the present research has set the stage for a more systematic explo-ration of collective dynamics of assemblies of map neu-rons by analytical means. We believe that the introducedtechniques can be successfully applied for treating theemergent behavior of populations in case of chemicallyand delay-coupled neurons [13]. Moreover, the methodmay likely be used to explore the effects of parameter in-homogeneity, as well as to study the impact of complexnetwork topologies [13, 15]. Our ultimate goal will be toextend the
M F approach to account for collective behav-ior of interacting populations of map neurons [13, 14].
ACKNOWLEDGMENTS
This work is supported by the Ministry of Education,Science and Technological Development of Republic ofSerbia under project No. 171017, and by the RussianFoundation for Basic Research under project No. 15-02-04245.
VI. APPENDIX
In the following, we provide the remaining details con-cerning the calculation of the S x dynamics, which isthe most complex part of the derivation of the effectivemodel. Following some algebra, Eq. (9) can be trans-formed to S x,n +1 = (1 − c ) S x,n + S y,n + σ − − c ) U n (11)+ ( (cid:104) G ( x i,n ) (cid:105) − (cid:104) G ( x i,n ) (cid:105) ) (cid:124) (cid:123)(cid:122) (cid:125) V ar ( G ( x i,n )) +2(1 − c )( (cid:104) x i,n G ( x i,n ) (cid:105)− m x,n (cid:104) G ( x i,n ) (cid:105) ) − (cid:104) y i,n G ( x i,n ) (cid:105) − m y,n (cid:104) G ( x i,n ) (cid:105) ) − β (1 − c ) [ (cid:104) x i,n H ( x i,n − d ) (cid:105) − m x,n (cid:104) H ( x i,n − d ) (cid:105) ] − β ( (cid:104) G ( x i,n ) H ( x i,n − d ) (cid:105) − (cid:104) G ( x i,n ) (cid:105)(cid:104) H ( x i,n − d ) (cid:105) )+ β ( (cid:104) H ( x i,n − d ) (cid:105) − (cid:104) H ( x i,n − d ) (cid:105) ) (cid:124) (cid:123)(cid:122) (cid:125) V ar ( H ( x i,n − d ))) . The partial results required for completing the calcu-lation are given by (cid:104) x i G ( x i ) (cid:105) − m x (cid:104) G ( x i ) (cid:105) = G (cid:48) ( m x ) S x − S x (12) (cid:104) y i G ( x i ) (cid:105) − m y (cid:104) G ( x i ) (cid:105) = − S x U xy − m x U xy + 2(1 + a ) m x U xy , where G (cid:48) ( m x ) ≡ − m x + 2(1 + a ) m x − a . Note that thetime indexes have been omitted for simplicity. After sometedious work, it may also be shown that the expressionfor variance V ar ( G ( x i )) reads V ar ( G ( x i )) = G (cid:48) ( m x ) S x + S x (cid:2) m x − a ) m x (13)+ 2(1 + a ) + 6 a (cid:3) + 15 S x . − β (1 − c ) [ (cid:104) x i H ( x i − d ) (cid:105) − (cid:104) x i (cid:105)(cid:104) H ( x i − d ) (cid:105) ] = − β × (14)(1 − c ) (cid:34)(cid:90) dx dx ...dx N N (cid:88) i x i H ( x i − d ) p ( x , ..., x N ) − m x (cid:90) dx dx ...dx N N (cid:88) i H ( x i − d ) p ( x , ..., x N ) (cid:35) = ... = − β (1 − c ) (cid:20)(cid:90) dx ( x − m x ) H ( x − d ) p ( x ) (cid:21) = − β (1 − c ) (cid:114) S x π exp (cid:20) − ( d − m x ) S x (cid:21) . Note that the second term containing the threshold func-tion has been evaluated in the main text, cf. Eq. (10).Finally, let us address the term β V ar ( H ( x i − d )),which can be estimated by considering theassociated expectation β V ar ( H ( x i − d )) ≈ β (cid:2) (cid:104) H ( x i − d ) (cid:105) − (cid:104) H ( x i − d ) (cid:105) (cid:3) . Applying the tech- nique introduced in Sec. II, we obtain E [ β H ( x i − d ) ] = β (cid:90) dx (cid:90) dx ... (cid:90) dx N × (15)( 1 N (cid:88) i (cid:88) j H ( x i − d ) H ( x j − d )) p ( x , x , ..., x N )= β N N (cid:90) dx H ( x − d ) p ( x ) (cid:124) (cid:123)(cid:122) (cid:125) N cases where i = j + β N N ( N − (cid:90) dx (cid:90) dx H ( x − d ) H ( x − d ) p ( x ) p ( x ) (cid:124) (cid:123)(cid:122) (cid:125) N ( N − cases where i (cid:54) = j = β N [1 − Erf [ d − m x √ S x ]] + β N N ( N − − Erf [ d − m x √ S x ]] . Given that β (cid:104) H ( x i − d ) (cid:105) = β (cid:104) − Erf [ d − m x √ S x ] (cid:105) , onearrives at β V ar ( H ( x i − d )) = β N (1 − Erf [ d − m x √ S x ])([1+ Erf [ d − m x √ S x ]) . (16)This shows that the variance of the threshold functionultimately contributes to a finite-size effect which can beneglected in the thermodynamic limit. [1] C. Zhou, L. Zemanova, G. Zamora, C. Hilgetag, and J.Kurths, Phys. Rev. Lett. , 238103 (2006).[2] E. Bullmore, and O. Sporns, Nat. Rev. Neurosci. , 186(2009).[3] D. Meunier, R. Lambiotte, and E. Bullmore, Front. Neu-rosci. , 200 (2010).[4] O. Sporns, D. Chialvo, M. Kaiser, and C. C. Hilgetag,Trends Cogn. Sci. , 418 (2004).[5] Dynamic Coordination in the Brain , edited by C. von derMalsburg, W. A. Phillips, and W. Singer, (MIT Press,Cambridge, 2010).[6] B. B. Averbeck, P. E. Latham, and A. Pouget, Nat. Rev.Neurosci. , 358 (2006).[7] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Mar-tinerie, Nat. Rev. Neurosci. , 229 (2001).[8] G. Buzs´aki, Rhythms of the Brain , (Oxford UniversityPress, Oxford, 2009).[9] J. L. P. Velazquez, and R. Wennberg,
Coordinated Activ-ity in the Brain: Measurements and Relevance to BrainFunction and Behavior , (Springer, New York, 2009).[10] R. T. Canolty et al., Science , 1626 (2006).[11] J. E. Lisman, and O. Jensen, Neuron , 1002 (2013).[12] Y. Baibolatov, M. Rosenblum, Z. Z. Zhanabaev, M. Kyz-garina, and A. Pikovsky, Phys. Rev. E , 046211 (2009).[13] O. V. Maslennikov, and V. I. Nekorkin, Phys. Rev. E ,012901 (2014).[14] O. V. Maslennikov, D. V. Kasatkin, N. F. Rulkov, andV. I. Nekorkin, Phys. Rev. E , 042907 (2013).[15] O. V. Maslennikov, V. I. Nekorkin, and J. Kurths, Phys.Rev. E , 042803 (2015). [16] O. V. Maslennikov, and V. I. Nekorkin, Commun. Non-linear Sci. Numer. Simul. , 10 (2015).[17] N. F. Rulkov, Phys. Rev. E , 041922 (2002).[18] N. F. Rulkov, I. Timofeev, and M. Bazhenov, J. Comput.Neurosci. , 203 (2004).[19] D. Q. Wei, and X. S. Luo, Europhys. Lett. , 68004(2007).[20] Q.Y. Wang, Z. Duan, M. Perc, and G. Chen, Europhys.Lett. , 50008 (2008).[21] C. A. S. Batista, A. M. Batista, J. A. C. de Pontes, R. L.Viana, and S. R. Lopes, Phys. Rev. E , 016218 (2007).[22] B. Ibarz, J. M. Casado, and M. A. F. Sanju´an, Phys.Rep. , 1 (2011).[23] I. Franovi´c, and V. Miljkovi´c, Europhys. Lett. , 68007(2010).[24] I. Franovi´c, and V. Miljkovi´c, Commun. Nonlinear Sci.Numer. Simul. , 623 (2011).[25] E. M. Izhikevich, Neural Comput. , 245 (2006).[26] E. M. Izhikevich, and G. M. Edelman, Proc. Natl. Acad.Sci. USA , 3593 (2008).[27] S. E. Folias, and P. C. Bressloff, Phys. Rev. Lett. ,208107 (2005).[28] C. R. Laing, W. C. Troy, B. Gutkin, and G. B. Ermen-trout, SIAM J. Appl. Math. , 62 (2002).[29] P. C. Bressloff, Phys. Rev. E , 051903 (2010).[30] M. A. Buice, and J. D. Cowan, Phys. Rev. E , 051919(2007).[31] N. Brunel, and V. Hakim, Neural Comput. , 1621(1999).[32] H. Hasegawa, Phys. Rev. E , 041903 (2003). [33] B. Lindner, J. Garcia-Ojalvo, A. Neiman, and L.Schimansky-Geier, Phys. Rep. , 321 (2004).[34] I. Franovi´c, K. Todorovi´c, N. Vasovi´c, and N. Buri´c,Phys. Rev. E , 022926 (2014).[35] I. Franovi´c, K. Todorovi´c, N. Vasovi´c, and N. Buri´c,Phys. Rev. E , 012922 (2013).[36] I. Franovi´c, K. Todorovi´c, N. Vasovi´c, and N. Buri´c,Chaos , 033147 (2012).[37] V. Klinshov, and I. Franovi´c, Phys. Rev. E , 062813(2015).[38] M. A. Zaks, X. Sailer, L. Schimansky-Geier, and A. B.Neiman, Chaos , 026117 (2005).[39] B. Sonnenschein, M. A. Zaks, A. B. Neiman, and L.Schimansky-Geier, Eur. Phys. J. Special Topics ,2517 (2013).[40] V. I. Nekorkin, and L. V. Vdovin, Izv. Vyssh. Uchebn.Zaved. Prikladn. Nelinejn. Din , 36 (2007).[41] M. Courbage, V. I. Nekorkin, and L. V. Vdovin, Chaos , 043109 (2007).[42] O. V. Maslennikov, and V. I. Nekorkin, Chaos , 023129(2013).[43] O. V. Maslennikov, and V. I. Nekorkin, ”Map-BasedApproach to Problems of Spiking Neural Network Dy-namics” in Nonlinear Dynamics and Complexity , V.Afraimovich, A. C. J. Luo, and X. Fu (Editors) (SpringerInternational Publishing Switzerland, 2014).[44] L. Arnold,
Random Dynamical Systems , (SpringerVerlag,Berlin, 1999). [45] J. A. Acebr´on, A. R. Bulsara, and W.-J. Rappel, Phys.Rev. E , 026202 (2004).[46] M. Gaudreault, J. M. Berbert, and J. Vi˜nals, Phys. Rev.E , 011903 (2011).[47] P. Kaluza, C. Strege, and H. Meyer-Ortmanns, Phys.Rev. E 82, 036104 (2010).[48] S. Tanabe and K. Pakdaman, Phys. Rev. E , 031911(2001).[49] V. I. Nekorkin, and O. V. Maslennikov, Radiophys.Quantum Electron. (Engl. Transl.) , 56 (2011).[50] M. Courbage, O. V. Maslennikov, and V. I. Nekorkin,Chaos Soliton. Fract. , 645 (2012).[51] O. V. Maslennikov, and V. I. Nekorkin, Radiophys.Quantum Electron. (Engl. Transl.) , 198 (2012).[52] C. W. Gardiner, Handbook of Stochastic Methods forPhysics, Chemistry and the Natural Sciences , 3rd ed.(Springer-Verlag, Berlin, 2004).[53]
Phase Response Curves in Neuroscience: Theory, Exper-iment, and Analysis , edited by N. W. Schultheiss, A. A.Prinz, and R. J. Butera, (Springer, New York, 2012).[54] P. A. Tass,
Phase Resetting in Medicine and Biology:Stochastic Modeling and Data Analysis (Springer, Berlin,Heidelberg, 2007).[55] E. M. Izhikevich,
Dynamical Systems in Neuroscience:The Geometry of Excitability and Bursting (MIT Press,Cambridge, 2007), Chap. 10.[56] C. C. Canavier, Scholarpedia1