Mean-field dynamics of a population of stochastic map neurons

TL;DR

This paper develops a mean-field model for a population of stochastic map neurons, accurately capturing diverse collective behaviors and responses to stimuli, with strong qualitative and quantitative agreement with the exact system.

Contribution

It introduces a cumulant-based mean-field approach that effectively models complex neuronal dynamics and responses, extending understanding of stochastic neural populations.

Findings

Model captures stability and bifurcations of the exact system.

Qualitative and quantitative agreement with the original dynamics.

Reproduces phase response curves and stimulus responses accurately.

Abstract

We analyze the emergent regimes and the stimulus-response relationship of a population of noisy map neurons by means of a mean-field model, derived within the framework of cumulant approach complemented by the Gaussian closure hypothesis. It is demonstrated that the mean-field model can qualitatively account for stability and bifurcations of the exact system, capturing all the generic forms of collective behavior, including macroscopic excitability, subthreshold oscillations, periodic or chaotic spiking and chaotic bursting dynamics. Apart from qualitative analogies, we find a substantial quantitative agreement between the exact and the approximate system, as reflected in matching of the parameter domains admitting the different dynamical regimes, as well as the characteristic properties of the associated time series. The effective model is further shown to reproduce with sufficient accuracy the phase response curves of the exact system and the assembly's response to external stimulation of finite amplitude and duration.

Figures

Equations

$$x_{n+1}$$

$$y_{n+1}$$

$$x_{i,n+1}$$

$$y_{i,n+1}$$

$$I_{i,n}^{syn}$$

$$⟨x_i^3⟩$$

$$⟨x_i^4⟩$$

$$⟨x_i^2{y}_i⟩$$

$$⟨x_i^3{y}_i⟩$$

$$⟨x_i^5⟩$$

$$⟨x_i^6⟩$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{i=1}^N{x}_{i,n}^k\approx E[x_{i,n}^k],$$

$$m_{x,n+1}=m_{x,n}-m_{y,n}+⟨G(x_{i,n})⟩-\beta ⟨H(x_{j,n}-d)⟩$$

$$⟨G(x_i)⟩=⟨-x_i^3+(1+a)x_i^2-a{x}_i⟩=G(m_x)+S_x(1+a-3m_x).$$

$$E[-\beta ⟨H(x_i-d)⟩]=\int d{x}_1\int d{x}_2...\int d{x}_N\times$$

$$(-\frac{\beta }{N}\sum _{i}H(x_i-d))p(x_1,x_2,...,x_N)=$$

$$-\beta {\int }_{-\infty }^{\infty }d{x}_{1}H(x_1-d)p(x_1)=-\beta {\int }_d^{\infty }\frac{1}{\sqrt{2\pi S_x}}{e}^{-\frac{(x_1-m_x{)}^2}{2S_x}}=$$

$$-\frac{\beta }{2}(1-Erf\left[\frac{d-m_x}{\sqrt{2S_x}}\right]),$$

$$S_{x,n+1}$$

$$-\beta H(x_{i,n}-d)-y_{i,n}+{\xi }_{i,n}+c{m}_{x,n}{]}^2⟩$$

$$-(m_{x,n}-m_{y,n}+G(m_{x,n})+S_{x,n}(1+a-3m_{x,n})$$

$$-\beta ⟨H(x_{i,n}-d)⟩{)}^2.$$

$$-2\beta E[⟨G(x_i)H(x_i-d)⟩-⟨G(x_i)⟩⟨H(x_i-d)⟩]=$$

$$-2\beta [\int d{x}_{1}G(x_1)H(x_1-d)p(x_1)$$

$$-\int d{x}_{1}H(x_1-d)p(x_1)\left[G(m_x)+S_x(1+a-3m_x)\right]]$$

$$\approx -2\beta [\int d{x}_1(G(m_x)+G^{\prime }(m_x)(x_1-m_x)+\frac{1}{2}{G}^{\prime \prime }(m_x)\times$$

$$(x_1-m_x{)}^2)H(x_1-d)p(x_1)-\int d{x}_{1}H(x_1-d)p(x_1)\times$$

$$\left[G(m_x)+S_x(1+a-3m_x)\right]]=...=$$

$$-2\beta \left[(1+a)(m_x+d)-a-3m_{x}d\right]\sqrt{\frac{S_x}{2\pi }}\exp \left[-\frac{(d-m_x{)}^2}{2S_x}\right].$$

$$m_{x,n+1}$$

$$-\frac{\beta }{2}(1-Erf\left[\frac{d-m_{x,n}}{\sqrt{2S_{x,n}}}\right])$$

$$m_{y,n+1}$$

$$S_{x,n+1}$$

$$+S_{x,n}(-3m_{x,n}^2+2(1+a)m_{x,n}-a{)}^2-2(1-c)\times$$

$$(3m_{x,n}^2{S}_{x,n}+3S_{x,n}^2-2(1+a)m_{x,n}{S}_{x,n}+a{S}_{x,n})$$

$$+2(3S_{x,n}{U}_n+3m_{x,n}^2{U}_n-2(1+a)m_{x,n}{U}_n)$$

$$-2\beta \left[(1+a)(m_{x,n}+d)-a-3d{m}_{x,n}\right]\sqrt{\frac{S_{x,n}}{2\pi }}\times$$

$$\exp \left[-\frac{(d-m_{x,n}{)}^2}{2S_{x,n}}\right]-2\beta (1-c)\sqrt{\frac{S_{x,n}}{2\pi }}\times$$

$$\exp \left[-\frac{(d-m_{x,n}{)}^2}{2S_{x,n}}\right]+S_{x,n}^2[36m_{x,n}^2-24(1+a)\times$$

$$m_{x,n}+2(1+a{)}^2+6a]+15S_{x,n}^3$$

$$S_{y,n+1}$$

$$U_{n+1}$$

$$-S_{y,n}-(U_n+ϵS_{x,n})(3S_{x,n}+3m_{x,n}^2$$

$$-2(1+a)m_{x,n})-\beta ϵ\sqrt{\frac{S_{x,n}}{2\pi }}exp\left[-\frac{(d-m_{x,n}{)}^2}{2S_{x,n}}\right].$$

$$S_{x,n+1}$$

$$+\underset{Var(G(x_{i,n}))}{\underbrace{(⟨G(x_{i,n}{)}^2⟩-⟨G(x_{i,n}){⟩}^2)}}+2(1-c)(⟨x_{i,n}G(x_{i,n})⟩$$

$$-m_{x,n}⟨G(x_{i,n})⟩)-2(⟨y_{i,n}G(x_{i,n})⟩-m_{y,n}⟨G(x_{i,n})⟩)$$

$$-2\beta (1-c)\left[⟨x_{i,n}H(x_{i,n}-d)⟩-m_{x,n}⟨H(x_{i,n}-d)⟩\right]$$

$$-2\beta (⟨G(x_{i,n})H(x_{i,n}-d)⟩-⟨G(x_{i,n})⟩⟨H(x_{i,n}-d)⟩)$$

Full text

Cortical connectivity patterns exhibit a hierarchical modular organization from the microscopic level of interacting neurons, via cortical columns and other types of mesoscopic circuitry, up to fibers projecting between the distributed brain areas ZZHK06 ; BS09 ; MLB10 . Architecture of neural assemblies comprising the anatomical modules closely reflects the functional specialization over different modalities SCKH04 , such that the macroscopic dynamics of neuronal populations, as well as the interplay of the associated collective modes, underpin various stages of information processing and higher cognitive functions MPS10 . In terms of evoked activity in the cortex, an ample example concerns the sensory regions, where many neurons display similar responses to a given stimulus, which strongly indicates that the information content is primarily encoded by the collective assembly response ALP06 . With regard to self-organized dynamics, the different forms of synchronization between spiking or bursting neurons give rise to macroscopic oscillations known to span several orders of magnitude in frequency range VLRM01 ; B09 . Such oscillations are deemed to facilitate coordinated rhythmic tasks and provide the dynamical background behind perception, attention, memory consolidation, motor planning or sleep VLRM01 ; B09 ; VW09 . Nonetheless, interaction between multiple rhythms has been indicated as critical for merging the activities of distant brain regions through cross-frequency coupling and other mechanisms CS06 ; LJ13 . For these reasons, gaining a deeper understanding of macroscopic behavior of neural populations has become an outstanding issue in neuroscience, focused primarily on the mechanisms guiding the emergent collective dynamics and the fashion in which assemblies respond to various external stimuli. Conceptually, such an approach is reminiscent of a common paradigm in nonlinear dynamics, where systems of coupled oscillators exhibiting a collective mode are typically treated as macroscopic oscillators, which may then be subjected to an external drive or can be influenced by collective rhythms from afferent populations BRZK09 .

In the present paper, we systematically analyze the emergent dynamics and the stimulus-response relationship of a population of stochastic map neurons using a mean-field ($MF$) approach. The considered map neurons can exhibit a variety of regimes, including excitability, subthreshold oscillations, regular and chaotic spiking or bursting, as well as mixed spiking-bursting oscillations MN14 ; MKRN13 ; MNK15 ; MN15 . Despite that the collective motion of spiking or bursting neurons subjected to noise has been extensively studied using different models of discrete local dynamics, such as Rulkov R02 ; RTB04 ; WL07 ; WDPC08 ; BPVL07 ; ICS11 ; FM10 ; FM11 or Izhikevich neuron maps I06 ; IE08 , a $MF$ theory for a population of stochastic map neurons is obtained here for the first time. Nevertheless, in case of continuous time systems, the $MF$ approach has been a standard analytical tool for treating diverse problems in neuroscience and other fields FB05 ; LTGE02 ; B10 ; BC07 ; BH99 ; H03 . Note that our derivation of the effective model relies on Gaussian approximation, which is introduced within the framework of a Gaussian closure hypothesis LGNS04 ; FTVB14 ; FTVB13 ; FTVB12 ; KF15 ; ZSSN05 ; SZNS13 .

The particular set of issues we address consists in establishing whether and how the $MF$ model can be used to $(i)$ qualitatively analyze the network stability and bifurcations of the exact system associated to emergence of generic macroscopic regimes; $(ii)$ provide adequate quantitative predictions in terms of bifurcation thresholds, and the average interspike intervals or bursting cycles of the exact system, as well as $(iii)$ accurately anticipate the population’s response to different forms of external stimuli. Within this context, it will be examined whether the effective model is capable of reproducing the properties of noise-activated, noise-induced and noise-perturbed modes of collective behavior.

The paper is organized as follows. In Section I, we make an overview of the local map dynamics and introduce the population model. Section II outlines the ingredients most relevant for the derivation of the $MF$ system, with the remaining technical details left for the Appendix. In Section III, the qualitative and quantitative agreement between the dynamics of the exact and the $MF$ model is illustrated by the appropriate bifurcation diagrams, as well as by comparing the characteristic features of the associated regimes. Section IV concerns the assembly’s stimulus-response relationship, first investigating the analogy between the respective phase-response curves ($PRCs$) of the exact system and the effective model in spiking and bursting regimes, and then considering the extent to which the $MF$ model reproduces the population’s response to rectangular pulses of finite amplitude and duration. In Section V, we provide a summary of our main results.

I Map neuron dynamics and the population model

The dynamics of an isolated neuron conforms to a map model first introduced in NV07 ; CNV07 , which is given by

[TABLE]

where $n$ denotes the iteration step. The variable $x_{n}$ qualitatively accounts for the membrane potential, whereas the recovery variable $y_{n}$, whose rate of change is set by a small parameter $\epsilon=10^{-2}$, mimics the behavior of ion-gating channels. The parameters $a,\beta$ and $d$ modify the profile of the ensuing oscillations, while $J$ crucially influences the neural excitability, viz. the transitions from silence to active regimes.

The $x_{n}$ evolution features two nonlinear terms, one being a FitzHugh-Nagumo-like cubic nonlinearity $G(x_{n})=x_{n}(x_{n}-a)(1-x_{n})$, which is complemented by a discontinuity term $-\beta H(x_{n}-d)$, where $H$ stands for the Heaviside step function. The parameters $a=0.1$ and $d=0.45$ are kept fixed throughout the paper. The impact of discontinuity consists in making the fast subsystem (Eq. (1) with $\epsilon=0$) a Lorenz-type map within certain parameter domains CNV07 ; MN13 , which endows the model with the ability to generate chaotic spike or burst oscillations, otherwise lacking in the Fitzhugh-Nagumo type of systems.

Under variation of $J$ and $\beta$, the map (1) may reproduce a rich repertoire of generic regimes displayed by the real neurons, as demonstrated in Fig. 1. In particular, the main frame shows amplitudes of the corresponding $x$ time series for the given $(J,\beta)$, while the remaining subfigures illustrate the characteristic waveforms pertaining to excitable regime (region $I$), subthreshold oscillations ($II$), regular ($III$) or chaotic spiking ($IV$), chaotic bursting ($V$), as well as the mixed chaotic spike-burst activity ($VI$). Some of the indicated boundaries, such as those involving domains $IV,V$ and $VI$ should be understood as tentative, since the associated transitions are smooth and therefore difficult to discern.

The detailed phase plane analysis concerning the relevant unstable invariant curves and the mechanisms underlying transitions between the different dynamical regimes can be found in MN14b . Here we briefly mention that under increasing $J$, the equilibrium loses stability via the Neimarck-Sacker bifurcation, which gives rise to subthreshold oscillations. Note that the latter may be considered as an excitable state, in a sense that a strong enough perturbation can elicit genuine spike, though the phase point does not relax to the equilibrium, but rather to a closed invariant curve.

Adopting model (1) for local dynamics, we focus on a 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 assembly, and is further influenced by synaptic noise from the embedding environment. The population activity is then described by the following system

[TABLE]

where $i$ specifies the particular neuron. The synaptic currents $I_{i,n}^{syn}$ comprise two types of terms. The diffusive couplings $I_{i,n}^{coup}$ are characterized by the strength $c$, which is assumed to be uniform over the network and is set to $c=1$ in the remainder of the paper. The random inputs $I_{i,n}^{rand}$ involve uncorrelated white noise ($E[\xi_{i,n}]=0,E[\xi_{i,n}\xi_{j,n^{\prime}}]=\delta_{ij}\delta(n-n^{\prime})$) of intensity $\sigma$.

Confined to a single unit, the stochastic component may influence its dynamics either by perturbing the deterministic oscillatory regimes, or by inducing oscillations in the excitable regime, cf. Fig. 2(b). The onset of noise-induced spiking or bursting within the parameter domain where the fixed point is deterministically stable (domain $I$ in Fig. 1) corresponds to a phenomenon of stochastic bifurcation A99 ; AB04 ; GBV11 ; KSM10 ; ZSSN05 . The latter are typically described phenomenologically, in a sense that certain time-averaged quantities, such as the asymptotic probability distributions of relevant variables or the associated power spectra, exhibit a qualitative change under variation of noise intensity. For instance, in continuous-time systems, it has been shown that the stochastic Hopf bifurcation from a stochastically stable fixed point to a stochastically stable limit cycle is accompanied by the loss of Gaussian property for the asymptotic distributions of the appropriate variables TP01 . At variance with standard deterministic bifurcations, where one clearly observes a critical value of the control parameter, the change of system’s behavior in noise-induced transitions is gradual ZSSN05 . Note that noise can also play an important part in the $(J,\beta)$ region $II$ where the deterministic map shows subthreshold oscillations. Here noise can give rise to a form of dynamics reminiscent of mixed-mode oscillations, cf. Fig. 2(c).

So far, models similar to (2) have been applied to address a number of problems associated to collective phenomena in networks of coupled neurons, including synchronization of electrically coupled units with spike-burst activity NM11 ; CMN12 , pattern formation in complex networks with modular architecture MKRN13 ; MN14 ; MN12 , transient cluster activity in evolving dynamical networks MN15 , as well as the basin stability of synchronization regimes in small-world networks MNK15 . Within this paper, the collective motion will be described in terms of the global variables $X_{n}=\frac{1}{N}\sum_{i=1}^{N}x_{i,n}$ and $Y_{n}=\frac{1}{N}\sum_{i=1}^{N}y_{i,n}$.

II Derivation of the mean-field model

Considering a $MF$ approximation, our main goal lies in deriving a reduced low-dimensional deterministic set of nonlinear difference equations whose dynamics is qualitatively analogous to the collective motion of the original system (2) comprised of $2N$ coupled stochastic maps. In particular, the $MF$ model should be able to generate all the regimes exhibited by the exact system, qualitatively reproducing the bifurcations that the latter undergoes. Also, applying the effective model, one should be capable of inferring with sufficient accuracy the parameter domains which admit the different collective states of the exact system, with the corresponding time series exhibiting similar characteristic quantitative features. Regarding the explicit effects of noise, the $MF$ model is expected to account for the onset or suppression of different types of collective modes associated to macroscopic spiking or bursting activity, which are mediated by synchronization or desynchronization of individual neuron dynamics, respectively. The synchronization processes may be influenced by noise in a variety of ways, including the scenarios where noise acts as a perturbation to mainly deterministic (and chaotic) local oscillations, or the ones where noise plays a facilitatory role, in a sense that the collective mode emerges via synchronization of noise-induced local dynamics.

Given that we consider a system of discrete-time equations, one cannot adopt the usual method of deriving the $MF$ model via Fokker-Planck formalism SZNS13 . Nevertheless, an analytically tractable $MF$ model may still be built by focusing on the evolution of cumulants LGNS04 ; FTVB14 ; FTVB13 ; ZSSN05 , whereby the full density of states is factorized into a series of marginal densities. The advantage of such an approach is that the simplifying approximations aimed at truncating the underlying cumulant series can be introduced in a controlled fashion. Such approximations, stated in a form of closure hypothesis LGNS04 , are required due to nonlinearity of the original system, which causes the dynamics of cumulants of the given order to be coupled to those of the higher order.

In our case, the derivation of the effective model incorporates an explicit Gaussian closure hypothesis LGNS04 ; FTVB14 ; FTVB13 ; ZSSN05 , by which all the cumulants above second order are assumed to vanish. The collective dynamics is then described by a set of five variables (the first- and second-order cumulants), including

• (i)

the means, given by $m_{x,n}=\lim\limits_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^{N}x_{i,n}\equiv\langle x_{i,n}\rangle$, $m_{y,n}=\lim\limits_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^{N}y_{i,n}\equiv\langle y_{i,n}\rangle$;

• (ii)

the variances, defined as $S_{x,n}=\langle x_{i,n}^{2}\rangle-\langle x_{i,n}\rangle^{2}=\langle x_{i,n}^{2}\rangle-m_{x,n}^{2}$ and $S_{y,n}=\langle y_{i,n}^{2}\rangle-\langle y_{i,n}\rangle^{2}=\langle y_{i,n}^{2}\rangle-m_{y,n}^{2}$;

• (iii)

the covariance $U_{n}=\langle x_{i,n}y_{i,n}\rangle-m_{x,n}m_{y,n}.$

The expressions for higher order moments $\langle x_{i,n}^{k}\rangle$ in terms of the first- and second-order cumulants G04 , such as

[TABLE]

can be derived using the closure hypothesis.

The Gaussian approximation effectively amounts to an assumption that the relation

[TABLE]

holds, whereby $E$ refers to expectation value obtained by averaging over an ensemble of different stochastic realizations. In other words, one supposes that the local variables are independent and are drawn from a normal distribution $\mathcal{N}(m_{x},S_{x})$. We do not know a priori whether such an assumption is fulfilled or not, but can only judge on its validity by verifying the correctness of the predictions on the population dynamics provided by the $MF$ model. Also note that the effective model concerns the assembly dynamics in the thermodynamic limit $N\rightarrow\infty$. The stochastic terms in this case can be neglected, as one may show them to contribute to finite size effects which scale as $1/N$. This means that the influence of noise in our $MF$ model is felt only via the noise intensity, which assumes the role of an additional bifurcation parameter.

Let us now illustrate the main technical points required for the derivation of the $MF$ model. Our focus will lie with a couple of relevant examples, whereas the remaining details are provided in the Appendix. We begin by considering the dynamics of the global variable $m_{x}$, which is given by

[TABLE]

It is easy to see that there is no contribution from the coupling term. As far as the third term on the r.h.s. of Eq. (5) is concerned, using Eq. (3), one arrives at

[TABLE]

In the last expression, we have dropped the time index for simplicity and have introduced the shorthand notation $G(m_{x})\equiv-m_{x}^{3}+(1+a)(m_{x}^{2}+S_{x})$.

The key problem is how to treat the final term in the r.h.s. of Eq. (5). Our approach consists in replacing the assembly average by the expectation value $(\langle H(x_{i}-d)\rangle\approx E[H(x_{i}-d)])$, obtained by assuming that the local variables at an arbitrary time moment are normally distributed according to $P(x_{i})\sim\mathcal{N}(m_{x},S_{x})$. The expectation may then be evaluated as

[TABLE]

with the error function $Erf(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt$. In the above calculation, we have explicitly used the assumption on the independence of distributions of local variables at any given moment of time.

In a similar fashion, one may consider the $S_{x}$ dynamics, which constitutes the most demanding part of the derivation. In particular, proceeding from the $S_{x}$ definition, we obtain

[TABLE]

As an illustration, let us evaluate one of the terms containing an average over the threshold function:

[TABLE]

Again, the time indexes have been suppressed to simplify the notation.

Leaving the remaining elements of the derivation for the Appendix, we now state the final equations of the $MF$ model in the thermodynamic limit

[TABLE]

III Analysis of stability and bifurcations

In this section, our goal is to demonstrate the qualitative and quantitative analogies between the dynamics of the exact system and the $MF$ model. To this end, we first examine the succession of macroscopic regimes in the $J-\beta$ parameter plane for $\sigma$ fixed at an intermediate value $\sigma=0.002$, see Fig. 3. As in case of a single unit, changing $J$ is relevant for the system’s excitability, viz. the transitions from silent to active regimes, while $\beta$ influences the waveforms of the active states (spiking, bursting, or mixed spike-bursting activity). The assembly is found to exhibit the collective modes which qualitatively correspond to the dynamics of a single unit illustrated in plates $III-VI$ of Fig. 1. The heat maps in the left column of Fig. 3 provide a comparison between the oscillation amplitudes $A$ of the global variable $X$ (top row) and the $MF$ variable $m_{x}$ (bottom row) for the given $(J,\beta)$. The right column indicates how well are matched the average interspike interval (or the average bursting cycle) $T$ of the exact system with the corresponding characteristics of the dynamics of the $MF$ model (II). In the given instances, exact system comprises an assembly of $N=100$ neurons, having obtained $A$ by averaging over a sufficiently long time series, whereas $T$ is determined by taking average over an ensemble of $20$ different stochastic realizations. With regard to $T$, we have selected a convenient threshold $\theta=0.2$, which allows a clear detection of individual spikes, and also enables one to unambiguously discern the initiation stage of bursts, as required for calculating the length of the bursting cycle.

Let us begin the analysis by focusing on the domain of $J$ values where the exact system exhibits the stochastically stable equilibrium, while the $MF$ model has a stable stationary state. The stochastic stability physically implies that fluctuations around the deterministic fixed point are typically of the order of noise, though some rare spikes may still be evoked. For $J$ sufficiently close to the region admitting the subthreshold oscillations, the population manifests macroscopic excitability. To properly illustrate this feature, we have analyzed the assembly dynamics in the limit $\sigma=0$, cf. Fig. 4. In particular, figures 4(a) and 4(b) show the maximum $X$ and $m_{x}$ values reached in the corresponding time series obtained for sets of different initial conditions $(X_{0},Y_{0})$ and $(m_{x,0},m_{y,0})$, respectively. The comparison between the two plots clearly corroborates that the boundary defining the domain of spiking response is appropriately anticipated by the $MF$ model. An important remark is that for the given $J$, the assembly may exhibit different forms of macroscopic excitability, generating a single spike or a burst of spikes, as dependent on the value of $\beta$. This is demonstrated by the time series in figures 4(c) and 4(d). The former refers to a one-spike response in case of $\beta=0.4$. For smaller $\beta$, one observes responses comprising two or more closely packed spikes, with Fig. 4(d) illustrating a three-spike burst encountered for $\beta=0.25$. Note that the time series of the full system and the $MF$ model are exactly matched in the limit $\sigma=0$.

Next we address the noise-influenced transitions from silence to active regimes observed under increasing $J$. To do so, in Figure 5(a) we have plotted the change of the firing (spiking or bursting) frequency $R$ for an assembly consisting of $N=100$ neurons. The average frequency is determined by considering an ensemble of $20$ different stochastic realizations, having $\sigma$ fixed to the moderate value from Fig. 4. The results from simulations of the full system (2) are compared against the data obtained for the $MF$ model. In this context, two points should be stressed. First, for moderate $\sigma$, note that the firing frequencies of the $MF$ model lie in close agreement to those of the exact system. As a second point, one finds that such quantitative agreement extends to different forms of collective behavior, viz. it holds for different types of transitions from silent to active regimes. As already indicated, the waveforms pertaining the active states depend on $\beta$, such that the associated transitions are mediated by the distinct synchronization processes. For instance, at $\beta=0$, synchronization involves time series of single units that conform to spiking activity of type $III$ from Figure 1, which are quite resilient to impact of noise. On the other hand, for $\beta=0.3$ or $\beta=0.4$, the individual units exhibit chaotic bursting or spiking activity, respectively, such that the underlying synchronization process may be more susceptible to stochastic effects. The typical $X$ time series illustrating the different collective modes are compared to the corresponding $m_{x}$ series in figures 5(b)-(e). The top (bottom) row concerns the data for the exact system ($MF$ model).

In order to investigate more closely the influence of noise for $J$ interval in vicinity of the transition from silence to active regimes, we examine how the profiles of $R(J)$ curves change under increasing $\sigma$. The results shown in Fig. 6 refer to $\beta=0.2$ and a population comprised of $N=100$ neurons. As expected, the transition appears quite sharp for moderate noise $\sigma=0.001$, but is considerably flattened for larger $\sigma$, e. g. $\sigma=0.05$. The crosses indicate the firing frequencies predicted by the $MF$ model for $\sigma=0.001$.

For larger $\sigma$, the $MF$ model fails to reproduce the behavior of the exact system in vicinity of threshold $J$, in a sense that it overestimates the maximal $R$ value, as well as the actual critical $J$ characterizing the transition. Viewed from another angle, one may infer that for sufficiently large $\sigma$ and $J$ below the threshold given by the $MF$ model, the latter fails to capture the impact of synchronization processes taking place between the noise-induced oscillations of individual units. This especially refers to $J$ interval where the spikes or bursts (depending on the given $\beta$) are superimposed on the background of subthreshold oscillations. An example of such a discrepancy between the behavior of the exact and the effective system is provided in Fig. 7, cf. Fig. 7(a) and Fig. 7(c). Also, for strong $\sigma$ and $J$ values above the transition, the firing frequencies anticipated by the effective model are typically higher than those of the exact system (not shown). Within this region, the stochastic effects suppress synchronization between the chaotic oscillations of single neurons, thereby reducing the corresponding $R$ value. This is not accounted for with sufficient accuracy by the $MF$ system. Note that such suppression of synchronization is reflected in the corresponding $X$ series by the spike (burst) ”skipping” mechanism, where the large-amplitude oscillations are occasionally replaced with subthreshold oscillations. For the associated $J$ and $\sigma$ values, such a phenomenon is absent in the dynamics of the effective model, cf. Fig. 7(b) and Fig. 7(d). In both of the scenarios illustrated in Fig. 7, the reason for having the $MF$ model fail lies in that the Gaussian approximation behind it breaks down due to large stochastic fluctuations.

The fashion in which the validity of the effective model’s predictions deteriorates with increasing $\sigma$ is made more explicit in Fig. 8, which shows the $A(J,\sigma)$ and $T(J,\sigma)$ dependencies for the exact and the approximate system at fixed $\beta=0.4$. The considered size of the network is $N=100$. Comparison between the respective $A$ (left column) and $T$ plots (right column) suggests that the range of $\sigma$ values where the $MF$ approximation applies is contingent on $J$. For instance, in the $J$ region below the deterministic threshold, one may estimate this range by noting that the effective bifurcation diagram in Fig. 8(a) indicates that noise-induced macroscopic oscillations emerge for $\sigma\approx 0.003$. Since this point is not adequately represented by the effective model, cf. Fig. 8(c), one may state that the Gaussian approximation breaks down around $\sigma\approx 0.003$ within the given $J$ region. Nevertheless, for $J$ above the deterministic threshold, the validity of the $MF$ model appears to depend rather strongly on particular $J$, with the $\sigma$ values where the Gaussian approximation effectively fails spanning the range $\sigma\in(0.002,0.006)$.

So far, we have investigated the impact of noise by comparing the results for the network of size $N=100$ to those obtained for the effective system. Nevertheless, within Section II, it has already been emphasized that the $MF$ model, deterministic in character, refers to the system’s behavior in the thermodynamic limit $N\rightarrow\infty$, whereas the explicitly stochastic terms could only be incorporated as finite-size effects. This makes it relevant to examine how the behavior of the exact system within the $J$ domain around deterministic threshold changes for large and fixed $\sigma$ under increasing $N$. To this end, we have plotted in Fig. 9 the $R(J)$ curves calculated for $N=100$ (squares), $N=500$ (circles) and $N=1500$ (diamonds) at fixed $\beta=0.2,\sigma=0.05$. The curve for $N=100$ evinces that the given $\sigma$ value is quite large in a sense of being sufficient to induce collective oscillations within the excitable regime. Apart from the dependencies for the full system, we also show the $R(J)$ curve associated to the $MF$ model (dashed line with crosses). An interesting point regarding the latter is that the $J$ threshold for the emergence of the collective mode is shifted toward a larger value compared to the case $\sigma\approx 0.01$. While the given transition itself appears quite sharp, the curves corresponding to the exact system approach it with increasing $N$, both in terms of the $J$ threshold and the $R$ values above the transition. This corroborates that the $(J,\sigma)$ domain where the Gaussian approximation behind the $MF$ model fails expectedly reduces with the increasing system size.

IV Response to external stimuli

The aim of this section is to investigate the extent to which the $MF$ model can be used to predict the stimulus-response relationship of an assembly exhibiting different macroscopic regimes, including the excitable state, as well as the spiking and bursting collective modes. Let us first focus on the two latter instances and examine the sensitivity of a population to an external pulse perturbation within the framework of phase resetting theory SPB12 ; T07 ; I07 ; C06 . In order to compare the behavior of the exact system and the effective model, we determine the corresponding phase resetting curves ($PRC$s), which describe the phase shift $\Delta\phi$, induced by the perturbation, in terms of the phase $\phi_{p}$ when the perturbation is applied. The considered 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}$ and width $\Delta=n_{i}-n_{f}$ are small compared to the amplitude and duration of the spiking (or bursting) cycle $T_{0}$, respectively. In case of the exact system, the same pulse current is delivered to each neuron $i$, adding the term $I_{p}$ to $x_{i}$ dynamics, whereas in the effective model, stimulation is administered via the $m_{x}$ variable. The phase $\varphi_{p}$ is defined in reference to $T_{0}$ by $\varphi_{p}=n_{p}/T_{0}$. The associated phase difference following the reset is calculated as $\Delta\varphi=1-T_{1}/T_{0}$, where $T_{1}$ denotes the duration of the perturbed spiking or bursting cycle.

The $PRC$s characterizing the assembly response in the spiking regime are provided in Fig. 10(a) and Fig. 10(b), whereby the former is obtained under the action of an excitatory ($a_{p}>0$), and the latter under the influence of an inhibitory stimulation ($a_{p}<0$). We stress that in both instances, the results derived from the effective model, denoted by circles, show excellent agreement with the data for the exact system (solid lines). In qualitative terms, one observes that an excitatory stimulation may advance the phase of the spiking cycle if it arrives sufficiently close to the spike, but still before the sharp rising stage. Nevertheless, an excitatory perturbation which acts during the spike or within the effective refractory period has a suppression effect, reflected in delaying of the next spike. At variance with the excitatory stimulation, the inhibitory pulse postpones the next firing time if it is introduced within the interval close to the rising stage of spike.

The $PRC$s determined for an assembly exhibiting collective bursting show qualitatively analogous effects to those described so far, see Fig. 10(c) and Fig. 10(d). This especially refers to impact of perturbation delivered 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 the associated $PRC$s expectedly exhibit strong fluctuations. Apart from that, one finds an interesting effect that both the excitatory and the inhibitory stimulation have a facilitatory role, i. e. cause phase advancement during the relaxation stage of the bursting cycle.

For a population in the excitable state, we consider scenarios where the system is influenced by a rectangular pulse perturbation of finite magnitude and duration, in a sense that the latter are comparable to corresponding features of typical spiking or bursting cycles. Note that the selected $J$ value $J=0.02$ lies sufficiently away from the interval admitting the subthreshold oscillations. Again, our objective is to determine whether the $MF$ model correctly anticipates the response of the exact system, now in presence of small to moderate noise. Some of the illustrative examples concerning the stimulus-response relationship under the finite perturbation are provided in Fig. 11. The top and the middle row refer to $X$ and corresponding $m_{x}$ time series, respectively, while the bottom row shows the profile of the applied stimulus. We find that in the absence of noise or for sufficiently small $\sigma$, the effective model reproduces the evoked behavior of the full system quite accurately. This also refers to some highly complex forms of responses, as corroborated in Fig. 11(a)-(c), which concern relatively large $a_{p}$ and $\Delta$. Under increasing $\sigma$, the ability of the $MF$ model to predict the dynamics of the exact system gradually reduces, but in a fashion that involves a nontrivial dependence on $\beta$. In particular, for smaller $\beta\approx 0.1$, which would facilitate macroscopic spiking mode for supercritical $J$, it turns out that the dynamics of the $MF$ model lies in close agreement to the one of the exact system even for moderate noise $\sigma=0.001$, cf. Fig. 11(d)-(f). However, for higher $\beta$, such an analogy between the responses of the exact and the $MF$ system is lost, see Fig. 11(g)-(i). Naturally, the validity of the predictions given by the $MF$ model deteriorates if the stimulation amplitude $a_{p}$ and the duration $\Delta$ are large, especially in presence of non-negligible noise.

V Summary and discussion

We have developed a $MF$ approach in order to systematically analyze the emergent dynamics and the input-output relationship of a population of stochastic map neurons. The reduced low-dimensional model has been derived within the framework of Gaussian approximation, formally introduced in a form of a closure hypothesis. In physical terms, such an approximation suggests that the local variables at an arbitrary moment of time are independent and conform to a normal distribution centered about the assembly mean and characterized by the associated assembly variance. Validity of such an approximation cannot be established a priori, but has been systematically verified by numerically corroborating that the $MF$ model reproduces the behavior of the exact system with sufficient accuracy.

In particular, we have first demonstrated that the effective model can qualitatively capture all the bifurcations of the exact system leading to the onset of different generic regimes of collective behavior. As far as the quantitative agreement is concerned, we have established substantial matching between the parameter domains admitting the respective dynamical regimes for the exact and the approximate system. Moreover, the typical features of the associated regimes, such as the average interspike interval or the average bursting cycle, exhibit analogous changes with parameter variation, and in many parameter domains display numerically similar values.

An important issue has been to explicitly examine how the effects of noise are reflected in the behavior of the $MF$ model. For the noise-perturbed activity, where the sufficiently small noise weakly influences the deterministic attractors of the system, the obtained results indicate that the Gaussian approximation holds. Nevertheless, the physical picture changes in case of noise-induced collective behavior. In particular, for different scenarios of stochastic bifurcations, typically corresponding to transitions from subthreshold oscillations, which involve generalized excitability feature, to spiking or bursting regimes, the exact system undergoes a gradual (smooth) change of collective dynamics, whereas the $MF$ model exhibits a standard deterministic bifurcation with a sharp bifurcation threshold. In such instances, the collective variables of exact system manifest large fluctuations, which explicitly violate the Gaussian approximation behind the effective model. Note that the loss of Gaussianity property for asymptotic distribution of relevant variables, which accompanies the described stochastic bifurcations, does not imply per se that our Gaussian approximation fails in the supercritical state. This point is evinced by the fact that the dynamics of the effective model shows qualitatively and quantitatively similar features to those of the exact system if the considered parameters lie sufficiently above the stochastic bifurcation. In fact, the Gaussian approximation applied in the derivation of the $MF$ model breaks down only in vicinity of such transitions, where the finite-size effects neglected in Eq. (II) become most prominent. We have numerically verified the prevalence of finite-size effects in these parameter domains, showing that the change of the appropriate order parameter, such as the spiking frequency, becomes sharper as the size of the neural assembly is increased. Nevertheless, the validity of Gaussian approximation is regained once the system is sufficiently above the bifurcation.

Apart from considering asymptotic dynamics, we have verified that the $MF$ model is capable of capturing the stimulus-response features of the exact system. For short pulse-like perturbations, it has been found that the approximate system reproduces the $PRC$s of the exact system for both the spiking and bursting regimes of collective activity with high accuracy. Substantial analogies have also been observed in case of macroscopic excitable regime for scenarios where the assembly is stimulated by rectangular pulse perturbations of finite amplitude and duration.

Having developed a viable $MF$ approach, the present research has set the stage for a more systematic exploration of collective dynamics of assemblies of map neurons by analytical means. We believe that the introduced techniques can be successfully applied for treating the emergent behavior of populations in case of chemically and delay-coupled neurons MN14 . Moreover, the method may likely be used to explore the effects of parameter inhomogeneity, as well as to study the impact of complex network topologies MN14 ; MNK15 . Our ultimate goal will be to extend the $MF$ approach to account for collective behavior of interacting populations of map neurons MKRN13 ; MN14 .

Acknowledgements.

This work is supported by the Ministry of Education, Science and Technological Development of Republic of Serbia under project No. $171017$, and by the Russian Foundation for Basic Research under project No. 15-02-04245.

VI Appendix

In the following, we provide the remaining details concerning the calculation of the $S_{x}$ dynamics, which is the most complex part of the derivation of the effective model. Following some algebra, Eq. (II) can be transformed to

[TABLE]

The partial results required for completing the calculation are given by

[TABLE]

where $G^{\prime}(m_{x})\equiv-3m_{x}^{2}+2(1+a)m_{x}-a$. Note that the time indexes have been omitted for simplicity. After some tedious work, it may also be shown that the expression for variance $Var(G(x_{i}))$ reads

[TABLE]

Let us now explicitly calculate the terms containing the threshold function. First we have

[TABLE]

Note that the second term containing the threshold function has been evaluated in the main text, cf. Eq. (II).

Finally, let us address the term $\beta^{2}Var(H(x_{i}-d))$, which can be estimated by considering the associated expectation $\beta^{2}Var(H(x_{i}-d))\approx\beta^{2}\left[\langle H(x_{i}-d)^{2}\rangle-\langle H(x_{i}-d)\rangle^{2}\right]$. Applying the technique introduced in Sec. II, we obtain

[TABLE]

Given that $\beta^{2}\langle H(x_{i}-d)\rangle^{2}=\frac{\beta^{2}}{4}\left[1-Erf[\frac{d-m_{x}}{\sqrt{2S_{x}}}]\right]^{2}$, one arrives at

[TABLE]

This shows that the variance of the threshold function ultimately contributes to a finite-size effect which can be neglected in the thermodynamic limit.