## Abstract

Bursting plays an important role in neural communication. At the population level, macroscopic bursting has been identified in populations of neurons that do not express intrinsic bursting mechanisms. For the analysis of phase transitions between bursting and non-bursting states, mean-field descriptions of macroscopic bursting behavior are a valuable tool. In this article, we derive mean-field descriptions of populations of spiking neurons and examine whether states of collective bursting behavior can arise from short-term adaptation mechanisms. Specifically, we consider synaptic depression and spike-frequency adaptation in networks of quadratic integrate-and-fire neurons. Analyzing the mean-field model via bifurcation analysis, we find that bursting behavior emerges for both types of short-term adaptation. This bursting behavior can coexist with steady-state behavior, providing a bistable regime that allows for transient switches between synchronized and nonsynchronized states of population dynamics. For all of these findings, we demonstrate a close correspondence between the spiking neural network and the mean-field model. Although the mean-field model has been derived under the assumptions of an infinite population size and all-to-all coupling inside the population, we show that this correspondence holds even for small, sparsely coupled networks. In summary, we provide mechanistic descriptions of phase transitions between bursting and steady-state population dynamics, which play important roles in both healthy neural communication and neurological disorders.

## 1  Introduction

The brain, composed of billions of single cells, has been demonstrated to possess a hierarchical, modular organization, indicative of a complex dynamical system (Ballard, 2015). Within this hierarchy, populations of neurons form functional entities, the states of which are defined by the collective dynamics of the population rather than by the activities of each single cell. Mean-field descriptions of the macroscopic dynamics of such populations are a valuable tool for the mathematical analysis of collective phenomena, as well as for computational models of multiple coupled populations of neurons. Population bursting is a particular mode of collective behavior that plays a major role in both healthy and pathological neural dynamics. At the single neuron level, bursting is characterized by the neuron firing a train of spikes, followed by a period of quiescence (Izhikevich, 2000). This behavior has been suggested to result from adaptive mechanisms introducing a slow timescale that enables dynamic regimes of bursting and controls the burst period (Izhikevich, 2000; Dhamala, Jirsa, & Ding, 2004). Mathematical descriptions of such adaptation mechanisms have been developed accordingly at the level of single cells. Importantly, bursting has also been reported in populations of cells without intrinsic bursting mechanisms (Izhikevich, 2000; Marder & Thirumalai, 2002; Zeldenrust, Wadman, & Englitz, 2018). In such cases, bursting can be conceived as a property of the collective dynamic interactions within the population, henceforth referred to as emergent bursting.

In healthy neural communication, emergent bursting activity may allow for a more reliable information transmission via chemical synapses (Lisman, 1997). This can be explained by the synchronized activity of the population during the burst, which stabilizes neural information transmission against different types of noise (Hahn, Ponce-Alvarez, Deco, Aertsen, & Kumar, 2019). Increased bursting, however, activity has been found in various neurological diseases, such as epilepsy and Parkinson's disease, and can act disruptively on neural communication if it exceeds certain levels of occurrence (Connors, 1984; Lobb, 2014). The mechanisms behind emergent bursting are not well understood, since most of the computational literature on bursting focuses on single cells (Guckenheimer, Harris-Warrick, Peck, & Willms, 1997; Izhikevich, 2000). Typical approaches to model bursting at the population level either use coupled circuits of excitatory and inhibitory populations (Kudela, Franaszczuk, & Bergey, 2003; Zeldenrust, Wadman, & Englitz, 2018) or include an explicit bursting mechanism such as the action of a neuromodulator (Marder & Thirumalai, 2002), feedforward inhibition (Zeldenrust & Wadman, 2013), or spike-frequency adaptation (SFA) (van Vreeswijk & Hansel, 2001). For example, van Vreeswijk and Hansel (2001) demonstrated in a network of coupled, excitatory leaky integrate-and-fire neurons that SFA can lead to the emergence of network bursting.

Importantly, these approaches investigate spiking neural networks, where the macroscopic state variables have to be inferred from single cell activities. However, a direct mathematical description of the macroscopic dynamics would be beneficial for both mathematical analyses of emergent bursting and studies on multiple coupled bursting populations. This is evident from the number of neural mass models that have been applied to study phase transitions between asynchronous and bursting states in coupled neural populations ((Suffczynski, Kalitzin, & Lopes Da Silva, 2004); (Breakspear et al., 2006; Spiegler et al., 2011; Chen et al., 2014); (Müller, van Albada, Kim, & Robinson, 2017)). These models describe macroscopic state variables such as average firing rates inside a population. Importantly, they have not been derived from a spiking neural network but were designed to resemble experimentally observed macroscopic features of neural behavior, such as the input-output relationships of a population or spectral features of population activity (Wilson & Cowan, 1972; Jansen & Rit, 1995; Robinson, Rennie, & Wright, 1997). Due to the continuous nature of such macroscopic state variables, neural mass models allow applying various methods from dynamical systems theory that could not be applied to spiking neural networks—for example, directly linking changes in the model structure to phase transitions. However, as the employed neural mass models are of a phenomenological nature, the mechanistic link to the underlying spiking neurons remains unclear in this type of work. While suited to investigate the emergence of bursting behavior within circuits of coupled neural populations, these models offer limited insight into the emergence of bursting within a population based on its intrinsic dynamics.

We find that emergent bursting can be found for all network sizes and coupling probabilities we considered. Furthermore, we demonstrate that the mean-field model generalizes well to networks of biologically plausible size and coupling probability. This makes our work applicable to a broad range of neurodynamic scenarios in which the impact of changes in population input or short-term adaptation strength on the behavior of spiking neuron populations is of interest.

## 2  Model Definition

The QIF neuron is the canonical form of type 1 neurons and has previously been used in combination with linear adaptation as a basis for models of bursting cells (Izhikevich, 2000). The evolution equation of the membrane potential $Vi$ of a single QIF neuron $i$ is given by
$τV˙i=Vi2+ηi+I(t)+Jsτ,$
(2.1)
$s=1N∑j=1N∑k∖tjk
(2.2)
with background current $ηi$, synaptic strength $J$, evolution time constant $τ$, extrinsic input $I(t)$, and synaptic activation $s$. The integral kernel $a(t-t')$ represents synaptic dynamics, for example, in the case of exponential synapses $a(t)=e-t/τs/τs$ with synaptic timescale $τs$. In the limit $τs→0$, $s$ represents instantaneous synaptic coupling in an all-to-all coupled network of $N$ neurons. A neuron $i$ emits its $k$th spike at time $tik$ when it reaches $Vθ$ upon which $Vi$ is reset to $Vr$. Recently, Montbrió et al. (2015) have shown that there exists an exact mean-field description for the macroscopic dynamics of the population dynamics given by equations 2.1 and 2.2 in the limit of $N→∞$ and $Vθ=-Vr→∞$. These can be expressed in terms of the evolution of the average firing rate $r$ and average membrane potential $v$:
$τr˙=Δπτ+2rv,$
(2.3)
$τv˙=v2+η¯+I(t)+Jrτ-(πrτ)2.$
(2.4)
Here, $η¯$ and $Δ$ are the center and half-width at half maximum of a Lorentzian distribution over the single neuron parameters $ηi$. Furthermore, synaptic conversion is assumed to be instantaneous ($τs→0$) throughout this article, which allows us to set $s=r$. However, we note that the adaptation mechanisms we consider here may lead to scenarios where $s≠r$ even in this limit. In these cases, we will state the functional relationship between $s$ and $r$ explicitly. For simplicity, we set the membrane time constant to $τ=1$ and define all other time-dependent variables used throughout this article in units of $τ$. While Figure 1 provides a visual representation of the relationship between the microscopic variables given by equations 2.1 and 2.2, and the macroscopic variables governed by equations 2.3 and 2.4, the mathematical details of this relationship are summarized in the following paragraphs.
Figure 1:

Relationship between microscopic and macroscopic representation of the QIF network. Microscopic network simulations were performed with $N=10,000$ neurons. Spikes are shown for 50 randomly selected neurons. Model parameters: $Δ=2$, $J=15Δ$, $η¯=-4Δ$, $I(t)=2.5$ if $10.0≤t<30$ and $I(t)=0$ otherwise.

Figure 1:

Relationship between microscopic and macroscopic representation of the QIF network. Microscopic network simulations were performed with $N=10,000$ neurons. Spikes are shown for 50 randomly selected neurons. Model parameters: $Δ=2$, $J=15Δ$, $η¯=-4Δ$, $I(t)=2.5$ if $10.0≤t<30$ and $I(t)=0$ otherwise.

As Montbrió et al. (2015) show the variables $V$ of the microscopic system always follow a Lorentzian distribution,
$ρ(V|η,t)=1πx(η,t)[V-y(η,t)]2+x(η,t)2,$
(2.5)
where $x(η,t)$ is associated with the mean firing rate $r(η,t)$ via $x(η,t)=πr(η,t)$, and $y(η,t)=v(η,t)$ is the membrane potential averaged over $V$ for specific values of $η$. The microscopic system follows the continuity equation,
$∂tρ+∂V[(V2+η+Js+I)ρ]=0,$
(2.6)
and by inserting equation 2.5 into 2.6, one obtains the dynamics of $x(η,t)$ and $y(η,t)$ in complex form,
$∂tw(η,t)=i[-w(η,t)2+η+Js(t)+I(t)],$
(2.7)
with $w(η,t)=x(η,t)+iy(η,t)$. Equation 2.7 can then be reduced to equations 2.3 and 2.4 using the Lorentzian distribution,
$g(η)=1πΔ(η-η¯)2+Δ2,$
(2.8)
given in complex form by
$πr˙+iv˙=iη¯-iΔ+Jr-(πr+iv)2.$
(2.9)
There exists a regime in which two stable states of the system defined by equations 2.3 and 2.4 can coexist, a fixed point representing low firing activity and a focus representing high firing activity. To introduce bursting to this system, a mechanism is required that alternatingly switches between those two states.

## 3  Synaptic Depression

The first short-term adaptation mechanism we consider is synaptic depression (SD), a multiplicative downscaling of the synaptic efficacy. Neurobiologically, this can represent various mechanisms such as postsynaptic receptor desensitization, alterations in the density of postsynaptic receptors, or resource depletion at the synapse (Jones & Westbrook, 1996; Zucker & Regehr, 2002).

### 3.1  Mathematical Definition of SD

To introduce SD to our system, we change equation 2.2 as follows:
$si=(1-Ai)1N∑j=1N∑k∖tjk
(3.1)
This adds a dependency of the postsynaptic activation $si$ on an adaptation variable $Ai$ pertaining to the $i$th postsynaptic neuron. The latter follows the microscopic evolution equations,
$Ai=∫0∞αN∑j=1N∑k∖tjk
(3.2)
$=αGA*r,$
(3.3)
$GA(t)=τA-2texp(-t/τA),$
(3.4)
which express adaptation as a convolution with an alpha kernel with rate $α$ and timescale $τA$. This choice accounts for (1) the delay between post-synaptic activation and peak adaptation, (2) the slow decay of the adaptation to baseline, and (3) the exponential shape of the rise and decay of adaptation that have been reported in experimental studies (Jones & Westbrook, 1996; Chung, Li, & Nelson, 2002; Hanson & Jaeger, 2002; Zucker & Regehr, 2002). The alpha kernel results from the biexponential kernel when the rise and decay time constants of the biexponential dynamics are identical (Robinson et al., 1997). This relationship allows application of our model to scenarios where the adaptation dynamics can be described by up to two different timescales, and it allows for an easy extension of our adaptation mechanism to differentiate between rise and decay times of adaptation. Importantly, this kind of adaptive mechanism has no effect on an isolated neuron, because it affects only the susceptibility to synaptic input, which is effectively zero shortly after it spiked (due to refractoriness). Therefore, any changes to the network behavior caused by SD have to emerge from the interaction between the network units. Furthermore, the adaptation is coupled to the mean-field firing rate of the population, that is, each spike in the network triggers postsynaptic adaptation at all network units. This becomes clear from equation 3.2, where $Ai$ is driven by the mean firing rate $r$ of the network and thus behaves like a global variable $A$ that can be written as equation 3.3, with $*$ denoting the convolution operation. Therefore, the synaptic input into each neuron also behaves like a global variable $s=(1-A)r$ with adaptation dynamics given by equations 3.3 and 3.4. Solving the convolution integral in equation 3.2 leads to the following evolution equations for the macroscopic adaptation dynamics:
$τAA˙=B,$
(3.5)
$τAB˙=-2B-A+αr.$
(3.6)
Hence, to introduce SD at the macroscopic level, we can change equations 2.3 and 2.4 to
$τr˙=Δπτ+2rv,$
(3.7)
$τv˙=v2+η¯+I(t)+Jrτ(1-A)-(πrτ)2.$
(3.8)
We demonstrate the accuracy of this mean-field description below.

### 3.2  Effects of SD

Next, we performed numerical bifurcation analysis of the four-dimensional system defined by equations 3.5 to 3.8. For different values of $α$, we initialized the model at a low activity state and continued the model in $η¯$. To this end, and for all other parameter continuations reported in this article, we used the software package AUTO-07p (Doedel et al., 2007). The adaptation timescale was chosen as $τA=10τ$, corresponding to slow adaptation relative to the evolution of the average membrane potential and firing rate. In accordance with the analysis of Montbrió et al. (2015), we found two fold bifurcations for $α=0$, defining the borders of a bistable regime in $η¯$ in which a stable node (representing low firing activity) and a stable focus (representing high firing activity) are separated by a saddle. For an increasing adaptation rate $α$, we identified a parameter regime in which two subcritical Andronov-Hopf bifurcations occur (see Figure 2A). The unstable limit cycle emerging from the branch of nodes (lower branch) exists only in a very narrow parameter range and gets annihilated quickly via a homoclinic bifurcation with the saddle. Thus, it is not of further interest for our examination of emergent bursting. As shown in Figure 2B, the unstable limit cycles emerging from the branch of foci (upper branch) via a subcritical Andronov-Hopf bifurcation undergo a fold bifurcation, marking the birth of a second, stable limit cycle. Continuation of the stable limit cycle in $η¯$ led to a second fold bifurcation. This fold bifurcation occurred at a value of $η¯$ close to the one where the stationary states of the system undergo a saddle-node bifurcation. The period of this unstable limit cycle grew rapidly toward infinity, terminating at a homoclinic bifurcation close to the fold bifurcation. Since this model behavior can, again, only be observed in a very small parameter range in our model, it is of limited relevance for the analysis of population bursting, which is why we omit a detailed analysis of this homoclinic bifurcation. However, this scenario has been reported in neural models with slow-fast dynamics before and was analyzed in detail in De Maesschalck and Wechselberger (2015).

Figure 2:

Bursting due to synaptic depression. (A) Bifurcation diagram of fixed points in $η¯$ for various values of $α$. Stable (unstable) fixed points are marked by solid (dotted) lines. (B) Subcritical Hopf bifurctions give rise to limit cycles representing bursting behavior. Minimum and maximum of the limit cycle are depicted in green. Vertical lines indicate initialization points used for panels C to E. (C) Transient bursting induced by excitation. (D) Transient bursting induced by inhibition. (E) In the bistable regime, excitatory and inhibitory stimuli switch the system between sustained bursting and sustained regular firing. Microscopic network simulations were performed with $N=10,000$ neurons. Spikes are shown for 50 randomly selected neurons. Model parameters: $Δ=2$, $J=15Δ$, $τ=1$, $τA=10$, $α=0.05$.

Figure 2:

Bursting due to synaptic depression. (A) Bifurcation diagram of fixed points in $η¯$ for various values of $α$. Stable (unstable) fixed points are marked by solid (dotted) lines. (B) Subcritical Hopf bifurctions give rise to limit cycles representing bursting behavior. Minimum and maximum of the limit cycle are depicted in green. Vertical lines indicate initialization points used for panels C to E. (C) Transient bursting induced by excitation. (D) Transient bursting induced by inhibition. (E) In the bistable regime, excitatory and inhibitory stimuli switch the system between sustained bursting and sustained regular firing. Microscopic network simulations were performed with $N=10,000$ neurons. Spikes are shown for 50 randomly selected neurons. Model parameters: $Δ=2$, $J=15Δ$, $τ=1$, $τA=10$, $α=0.05$.

As shown in Figure 2B, the stable regime of the bursting limit cycle can coexist with the high-activity focus and hence permits various transitions between bursting and steady-state behavior. Figures 2C and 2D demonstrate that the stable bursting state can be transiently entered from a either low-activity state through excitation (see Figure 2C) or a high-activity state through inhibition (see Figure 2D). Furthermore, the bistable regime allows for hysteresis—switching between limit cycle and focus equilibrium through transient excitatory and inhibitory inputs (see Figure 2E). In neural communication, this regime is particularly relevant, since it allows for quick transitions between highly different firing modes via transient inputs and introduces a form of network memory. However, it is also of interest for pathological neural dynamics such as observed in epilepsy, which have been proposed to reflect switching between a healthy state of low neural synchrony and a coexisting pathological, synchronous state (Suffczynski, Kalitzin, & Lopes Da Silva, 2004; Takeshita, Sato, & Bahar, 2007). Importantly, Figures 2C to 2E show a close correspondence between numerical simulations of the mean-field model and the spiking neural network.

### 3.3  Limit Cycle Characteristics

For a better understanding of the bistable regime, we mapped out the basins of attraction with respect to the state variables $A$, $r$, and $v$ of the model given by equations 3.5, 3.7, and 3.8, respectively. Figure 3 visualizes different trajectories of this three-dimensional projection of the system when initialized at different points near the unstable limit cycle that separates the stable limit cycle (green) and the stable focus (purple). It is the unstable (saddle) limit cycle and its stable manifold (i.e., all points in state space from which the system converges onto the unstable limit cycle) that act as a separatrix between the stable focus and the stable limit cycle. This separating behavior is visible where the unstable manifold is orthogonal to the 3D-projection in Figure 3, especially along the left part of the unstable limit cycle.

Figure 3:

Coexistence of bursting and steady-state behavior. Left column: three-dimensional projection (onto $r$, $v$, and $A$) of the four-dimensional state-space representation of the system dynamics. For these parameters, the stable limit cycle (bold green curve, representing bursting behavior) coexists with a stable focus (purple dot) and an unstable limit cycle (black dashed curve). Thin curves mark trajectories with different initial conditions in the basin of attraction of the limit cycle (green) or the focus (purple). Right column: Two sample time series that have been initiated in either the basin of attraction of the stable limit cycle ($r0=1.8$, $v0=1.0$, $A0=0.4$, $B0=0.01$) or the basin of attraction of the stable focus ($r0=0.75$, $v0=-0.4$, $A0=0.36$, $B0=0.0$). Model parameters: $Δ=2$, $J=15Δ$, $τ=1.0$, $τA=10$, $α=0.05$, $η¯=-4.6$.

Figure 3:

Coexistence of bursting and steady-state behavior. Left column: three-dimensional projection (onto $r$, $v$, and $A$) of the four-dimensional state-space representation of the system dynamics. For these parameters, the stable limit cycle (bold green curve, representing bursting behavior) coexists with a stable focus (purple dot) and an unstable limit cycle (black dashed curve). Thin curves mark trajectories with different initial conditions in the basin of attraction of the limit cycle (green) or the focus (purple). Right column: Two sample time series that have been initiated in either the basin of attraction of the stable limit cycle ($r0=1.8$, $v0=1.0$, $A0=0.4$, $B0=0.01$) or the basin of attraction of the stable focus ($r0=0.75$, $v0=-0.4$, $A0=0.36$, $B0=0.0$). Model parameters: $Δ=2$, $J=15Δ$, $τ=1.0$, $τA=10$, $α=0.05$, $η¯=-4.6$.

In a next step, we performed a two-parameter continuation of the subcritical Andronov-Hopf bifurcation in $η¯$ and $α$ to examine the dependence of population bursting on the interplay between network excitation and adaptation rate. Figure 4 shows that dynamic regimes of stable bursting can be found for a range of the two parameters, which is bounded by the fold of limit cycle bifurcations. For $η¯$, the parameter range in which the limit cycle exists corresponds to most of the cells in the population being in an excitable regime and has been reported for a number of models using QIF neurons (Montbrió et al., 2015; Ratas & Pyragas, 2016; Schmidt, Avitabile, Montbrió, & Roxin, 2018). Within this range, the interburst frequency scales with the input strength, which makes the latter an interesting parameter for tuning this model to reflect experimentally reported interburst frequencies. In summary, we identified SD as a potential mechanism for bursting to occur in networks of globally coupled spiking neurons. We demonstrated that this bursting mechanism could be transiently switched on and off via transient input currents and found that the interburst frequency can be tuned via the input strength.

Figure 4:

Existence and period of bursting. Lines indicate two-parameter continuations of codimension 1 bifurcations in the ($η¯,α$) plane. The color-coded region shows the interburst period of the stable limit cycle (depicted in units of $τ$). Other model parameters: $Δ=2$, $J=15Δ$, $τ=1.0$, $τA=10$.

Figure 4:

Existence and period of bursting. Lines indicate two-parameter continuations of codimension 1 bifurcations in the ($η¯,α$) plane. The color-coded region shows the interburst period of the stable limit cycle (depicted in units of $τ$). Other model parameters: $Δ=2$, $J=15Δ$, $τ=1.0$, $τA=10$.

In this section, we examine how well our results for multiplicative adaptation generalize to additive short-term adaptation mechanisms, known as spike-frequency adaptation (SFA). SFA differs from the above described adaptation mechanism in two aspects: (1) it affects the presynaptic activity instead of the postsynaptic efficacy, and (2) it acts additively instead of multiplicatively (van Vreeswijk & Hansel, 2001; Gigante et al., 2007).

### 4.1  Mathematical Defition of SFA

SFA is a homeostatic mechanism that acts at the single cell level via spike-triggered balancing currents (Guckenheimer et al., 1997; Benda & Herz, 2003). As such, SFA is an adaptive mechanism driven by the firing rate of a single cell rather than the firing rate of the whole network. Therefore, we introduce neuron-specific adaptation variables $Ai$ and $Bi$,
$τAA˙i=Bi,$
(4.1)
$τAB˙i=-2Bi-Ai+αri,$
(4.2)
with neuron-specific spiking activity $ri$ given by
$ri=∑k∖tik
(4.3)
Adding the adaptation variable $Ai$ to 2.1, we receive the following evolution equation for the membrane potential of the single neuron:
$τV˙i=Vi2+ηi+I(t)-Ai+Jsτ,$
(4.4)
where we can set $s=r$, since the synaptic input is not directly affected by adaptation anymore. If $τA≫τ$, then we may assume that the adaption variable $Ai$ changes very slowly in comparison to $Vi$, and that the Lorentzian ansatz, equation 2.5, holds. Effectively, the variable $Ai$ can be regarded as constant and be absorbed into $ηi$ in this limit (for a similar approach, see Gigante et al. (2007)). We note here that the Lorentzian ansatz equation 2.5 is independent of the distribution of ${ηi}$. The neuron-specific spiking activity $ri$ is associated with $x(η,t)$ via $ri=x(ηi,t)/π$. Hence, the microscopic dynamics are given by
$∂tw=i[η+Jr+I-w2-αGA*ℜ[w]/π].$
(4.5)
If $g(η)$ follows the Lorentzian distribution, then equation 4.5 results in
$πr˙+iv˙=i[η¯-iΔ+Jr-(πr+iv)2-αGA*r].$
(4.6)
Thus, the presynaptic additive model is identical to a postsynaptic additive model at the macroscopic scale if $τA≫τ$. At this scale, postsynaptic SFA is straightforward to realize by adding $A$ to the right-hand side of equation 2.3, leading to the macroscopic evolution equations
$τr˙=Δπτ+2rv,$
(4.7)
$τv˙=v2+η¯+I(t)-A+Jrτ-(πrτ)2,$
(4.8)
where $s=r$ and $A$ is still given by equations 3.5 and 3.6. Note that equations 4.7 and 4.8 are equivalent to the complex differential form given by equation 4.6.

### 4.2  Effects of SFA

Using the system defined by equations 3.5, 3.6, 4.7, and 4.8, we repeated the parameter continuation in $η¯$ for different values of $α$. This was done to examine whether the results we obtained for SD-induced population bursting would translate to an additive adaptation mechanism. As can be seen in Figure 5A, we found results strikingly similar to the ones we found for SD. For sufficiently strong levels of SFA (parameterized via $α$), we found a subcritical Andronov-Hopf bifurcation marking the birth of a bursting limit cycle. Furthermore, we again found a bistable region in which the bursting limit cycle coexists with the stable focus, separated by an unstable limit cycle. These regimes could be well traversed via transient inputs, as depicted in Figure 5B. Driving the microscopic model given by equations 2.2 and 4.1 to 4.4 with the same transient inputs, we found that the spiking dynamics were still attracted to the low-dimensional manifold described by the macroscopic system. This shows that even with $τA=10τ$, the condition $τA≫τ$ is sufficiently satisfied for the macroscopic description of the population bursting dynamics to be valid.

Figure 5:

Bursting due to SFA. (A) Similar to SD, SFA changes the fixed-point structure and stability of the system via Hopf bifurcations. The limit cycle minima and maxima are visualized in green. Stable (unstable) equilibria are marked by solid (dotted) lines. The dashed vertical line marks the initialization point used for panel B. (B) In the bistable regime, positive (red) and negative (blue) stimuli switch the system between sustained bursting and sustained regular firing. Model parameters: $Δ=2$, $J=15Δ$, $τ=1$, $τA=10$, $α=1.0$.

Figure 5:

Bursting due to SFA. (A) Similar to SD, SFA changes the fixed-point structure and stability of the system via Hopf bifurcations. The limit cycle minima and maxima are visualized in green. Stable (unstable) equilibria are marked by solid (dotted) lines. The dashed vertical line marks the initialization point used for panel B. (B) In the bistable regime, positive (red) and negative (blue) stimuli switch the system between sustained bursting and sustained regular firing. Model parameters: $Δ=2$, $J=15Δ$, $τ=1$, $τA=10$, $α=1.0$.

## 5  Finite Size Effects

In the previous sections, we examined bursting behavior in a mean-field model of a network of QIF neurons. This model has been derived in the thermodynamic limit ($N→∞$) and under the constraint of all-to-all coupling. We demonstrated that the mean-field model is an accurate representation of the macroscopic behavior of an all-to-all coupled QIF neuron network of size $N=104$. In this section, we investigate how well our findings generalize to more realistic network architectures. Specifically, we ask how sensitive our findings are with respect to the number of cells and the connection probabilities inside the QIF network. For this purpose, we constructed QIF networks with varying numbers of QIF neurons $N$ and varying coupling probabilities $p$ and compared their dynamic behavior against the mean-field model given by equations 3.5 to 3.8. We initialized each model in a stable bursting regime ($η¯=-5.5$, $Δ=2$, $J=15Δ$, $τA=10$, $α=0.05$) and performed numerical simulations over a time interval of $T=1000$. For comparison with the mean-field model, we calculated the average frequency of bursts in time intervals of $T=100$, as well as the peak amplitude of the bursting. The results of this procedure are visualized in Figure 6.

Figure 6:

Finite size effects. (A) Difference between mean-field model and spiking neural network in interburst frequency and maximum bursting amplitude for different network sizes ($N$) and coupling probabilities ($p$) of the spiking neural network. (B) Sample time series of the mean-field model and spiking neural network for specific $N$ and $p$. Model parameters: $Δ=2$, $η¯=-5.5$, $J=15Δ$, $τ=1$, $τA=10$, $α=0.05$.

Figure 6:

Finite size effects. (A) Difference between mean-field model and spiking neural network in interburst frequency and maximum bursting amplitude for different network sizes ($N$) and coupling probabilities ($p$) of the spiking neural network. (B) Sample time series of the mean-field model and spiking neural network for specific $N$ and $p$. Model parameters: $Δ=2$, $η¯=-5.5$, $J=15Δ$, $τ=1$, $τA=10$, $α=0.05$.

Most important, we found stable bursting behavior in each of the microscopic networks we investigated, even in the smallest, most sparsely connected one ($N=1000$, $p=0.01$). Furthermore, we observed that the bursting frequency scales negatively with coupling probability, whereas increases of the network size improved the correspondence with the mean-field model (see Figure 6A). This led to the interesting behavior that the correspondence between small networks (e.g., with $N=1000$ neurons) and the mean-field model had an optimum at intermediate coupling probabilities (e.g., $p=0.1$). Regarding the bursting amplitude, we observed that both a small network size and a small coupling probability led to a reduced amplitude. Investigating the simulated time series more closely, we found that this was due to a decreased within-burst neural synchronization (see Figure 6B).

In summary, we found that our mean-field model generalizes well to networks with realistic cell counts and coupling probabilities, even though it was derived for the special case of an infinitely large network with all-to-all coupling. Most important, our finding that states of population bursting can emerge from the interaction of recurrent, excitatory coupling and short-term plasticity has been reproduced in each microscopic network we examined. Discrepancies between the mean-field model and the microscopic network were found for small, sparsely coupled networks in the interburst frequency and within-burst neural synchrony. Interpreting those quantities in the mean-field model must thus be handled cautiously in cases where the underlying microscopic network may be of small size or very sparsely coupled.

## 6  Discussion

In this work, we have examined the dynamic impact of two different short-term adaptation mechanisms on the collective behavior of a globally coupled QIF population: synaptic depression and spike-frequency adaptation. For both mechanisms, we derived and validated mean-field descriptions of the macroscopic dynamics via the approach described in Montbrió et al. (2015). Using bifurcation analysis, we identified and characterized regimes of collective bursting that emerged given a sufficiently strong adaptation rate. These bursting regimes could coexist with nonbursting regimes, allowing for dynamic phase transitions between bursting and steady-state behavior via transient inputs.

### 6.1  Population Bursting in Health and Disease

In neural communication, bistable regimes are particularly relevant, since they allow for quick transitions between highly different firing modes via transient inputs and hence can be used to describe short-term or working memory as a dynamic property of the collective behavior of neural populations (Kunze, Peterson, Haueisen, & Knösche, 2017; Schmidt et al., 2018). However, these bistable bursting regimes are also of interest for pathological neural dynamics such as observed in epilepsy and Parkinson's disease, which have been proposed to reflect switching between an asynchronous, healthy state and a coexisting synchronous, pathological state with increased bursting behavior (Connors, 1984; Suffczynski et al., 2004; Takeshita et al., 2007; Lobb, 2014). In Parkinson's disease, for example, changes in the input strength and the synaptic short-term plasticity have been reported together with a strong increase of the bursting activity inside key neural populations known to be affected by the disease (located inside the basal ganglia) (Magill, Bolam, & Bevan, 2001; Hanson & Jaeger, 2002; Wichmann & Soares, 2006; Baufreton & Bevan, 2008). Based on these findings, it has been suggested that alterations of synaptic short-term depression might have a critical influence on the recurrent coupling inside the basal ganglia populations, thus leading to more synchronized, bursting neural activity. Since, in our model, the existence of bursting depends directly on strength of population input and synaptic depression, it could serve as a mechanistic link between those observations and be embedded in larger network models of Parkinsonian neural dynamics. Further support for this mechanistic link can be drawn from in vitro studies explaining bursting in neural populations via an interaction between alterations in the average population input and postsynaptic homeostatic plasticity (Zierenberg, Wilting, & Priesemann, 2018). In summary, this renders our mean-field model applicable to a broad range of neurodynamic scenarios that involve generation of bursting behavior in neural populations. As discussed above, a prime example would be the excitatory-inhibitory interactions between subthalamic nucleus and globus pallidus, which have been suggested as a generator of Parkinsonian oscillations inside the basal ganglia and are known to show altered synaptic depression patterns in Parkinson's disease (Magill et al., 2001; Hanson & Jaeger, 2002; Wichmann & Soares, 2006; Baufreton & Bevan, 2008).

### 6.2  Neurophysiological Interpretation of the Mean-Field Model

Our model provides a description of the emergence of synchronous, bursting neural dynamics in recurrently connected populations of spiking neurons that could either arise from spike-frequency adaptation (additive) or postsynaptic efficacy reduction (multiplicative). Experimental work suggests that the former is a result of different balancing currents triggered at a single cell after it generated a spike (Fuhrmann, Markram, & Tsodyks, 2002; Benda & Herz, 2003). The latter has been linked to various mechanisms such as receptor desensitization (Jones & Westbrook, 1996; Wong, Graham, Billups, & Forsythe, 2003), receptor density reduction (Turrigiano, 2008; Pozo & Goda, 2010), or resource depletion at glial cells involved in synaptic transmission (Virkar, Shew, Restrepo, & Ott, 2016; Huang, Chang, Chen, Lai, & Chan, 2017). Even though these adaptation mechanisms can express tremendously different timescales, ranging from a few hundred milliseconds (e.g., spike-frequency adaptation; Fuhrmann et al., 2002) to days (e.g., postsynaptic receptor density reduction; Pozo & Goda, 2010), our mean-field descriptions remain applicable. However, note that our model of synaptic depression cannot express vesicle depletion at the presynaptic site, as introduced for single cell models in Tsodyks, Pawelzik, and Markram (1998). For this type of plasticity, it is not clear whether the Lorentzian ansatz can be applied, making a derivation of the mean-field equations more difficult. Apart from this restriction, our mean-field model is independent of the particular form of adaptation that is used. That is, the convolution with an alpha kernel we employed as a second-order approximation of the dynamics of the adaptation variable $A$ could be replaced with any other description. This allows future studies the examination of the influence of specific short-term adaptation characteristics on population dynamics. Finally, we tested the validity of our mean-field model against two of its core assumptions: infinitely large neural populations and all-to-all coupling inside the population. We were able to demonstrate that short-term depression led to population bursting even for very small and sparsely coupled networks. While the bursting frequency and amplitude did scale with the network size and coupling probability, the qualitative behavior of the model as predicted by our bifurcation analysis remained the same. Hence, we conclude that our model is well capable of describing phase transitions between steady-state and bursting behavior in spiking neural networks of realistic size and coupling density.

## Acknowledgments

We thank E. Montbrió and H. G. E. Meijer for valuable comments and discussions. H.S. was supported by the German Research Foundation (DFG (KN 588/7-1) awarded to T.R.K., via Priority Program 2041, “Computational Connectomics”). R.G. was supported by the Max Planck Society and is currently funded by the Studienstiftung des Deutschen Volkes.

## References

Ballard
,
D. H.
(
2015
).
Brain computation as hierarchical abstraction
.
Cambridge, MA
:
MIT Press
.
Baufreton
,
J.
, &
Bevan
,
M. D.
(
2008
).
D2-like dopamine receptor-mediated modulation of activity-dependent plasticity at GABAergic synapses in the subthalamic nucleus.
Journal of Physiology
,
586
(8)
,
2121
2142
. doi:10.1113/jphysiol.2008.151118
Benda
,
J.
, &
Herz
,
A. V. M.
(
2003
).
A universal model for spike-frequency adaptation.
Neural Computation
,
15
(11)
,
2523
2564
. doi:10.1162/089976603322385063
Breakspear
,
M.
,
Roberts
,
J. A.
,
Terry
,
J. R.
,
Rodrigues
,
S.
,
Mahant
,
N.
, &
Robinson
,
P. A.
(
2006
).
A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis.
Cerebral Cortex
,
16
(9)
,
1296
1313
. doi:10.1093/cercor/bhj072
Chen
,
M.
,
Guo
,
D.
,
Wang
,
T.
,
Jing
,
W.
,
Xia
,
Y.
,
Xu
,
P.
, &
Yao
,
D.
(
2014
).
Bidirectional control of absence seizures by the basal ganglia: A computational evidence.
PLOS Computational Biology
,
10
(3)
. doi:10.1371/journal.pcbi.1003495
Chung
,
S.
,
Li
,
X.
, &
Nelson
,
S. B.
(
2002
).
Short-term depression at thalamo-cortical synapses contributes to rapid adaptation of cortical sensory responses in vivo.
Neuron
,
34
(3)
,
437
446
. doi:10.1016/S0896-6273(02)00659-1
Connors
,
B. W.
(
1984
).
Initiation of synchronized neuronal bursting in neo- cortex.
Nature
,
310
(
5979
),
685
. doi:10.1038/310685a0
De Maesschalck
,
P.
, &
Wechselberger
,
M.
(
2015
).
Neural excitability and singular bifurcations.
Journal of Mathematical Neuroscience
,
5
(1)
,
16
. doi:10.1186/s13408-015-0029-2
Dhamala
,
M.
,
Jirsa
,
V. K.
, &
Ding
,
M.
(
2004
).
Transitions to synchrony in coupled bursting neurons.
Physical Review Letters
,
92
(2)
,
028101
. doi:10.1103/PhysRevLett.92.028101
Doedel
,
E. J.
,
Fairgrieve
,
T. F.
,
Sandstede
,
B.
,
Champneys
,
A. R.
,
Kuznetsov
,
Y. A.
, &
Wang
,
X.
(
2007
).
Auto-07p: Continuation and bifurcation software for ordinary differential equations
(Technical Report)
.
Montreal
:
Department of Computer Science, Concordia University
.
Fuhrmann
,
G.
,
Markram
,
H.
, &
Tsodyks
,
M.
(
2002
).
Spike frequency adaptation and neocortical rhythms.
Journal of Neurophysiology
,
88
(2)
,
761
770
. doi:10.1152/jn.2002.88.2.761
Gigante
,
G.
,
Mattia
,
M.
, &
Giudice
,
P. D.
(
2007
).
Diverse population-bursting modes of adapting spiking neurons.
Physical Review Letters
,
98
(14)
,
148101
. doi:10.1103/PhysRevLett.98.148101
Guckenheimer
,
J.
,
Harris-Warrick
,
R.
,
Peck
,
J.
, &
Willms
,
A.
(
1997
).
Bifurcation, bursting, and spike frequency adaptation.
Journal of Computational Neuroscience
,
4
(3)
,
257
277
. doi:10.1023/A:1008871803040
Hahn
,
G.
,
Ponce-Alvarez
,
A.
,
Deco
,
G.
,
Aertsen
,
A.
, &
Kumar
,
A.
(2019)
.
Portraits of communication in neuronal networks.
Nature Reviews Neuroscience
,
20
(2)
,
117
. doi:10.1038/s41583-018-0094-0
Hanson
,
J. E.
, &
Jaeger
,
D.
(2002)
.
Short-term plasticity shapes the response to simulated normal and Parkinsonian input patterns in the globus pallidus.
Journal of Neuroscience
,
22
(12)
,
5164
5172
. doi:10.1523/JNEUROSCI.22-12-05164.2002
Huang
,
Y.-T.
,
Chang
,
Y.-L.
,
Chen
,
C.-C.
,
Lai
,
P.-Y.
, &
Chan
,
C. K.
(
2017
).
Positive feedback and synchronized bursts in neuronal cultures.
PLOS One
,
12
(11)
,
e0187276
. doi:10.1371/journal.pone.0187276
Izhikevich
,
E. M.
(2000)
.
Neural excitability, spiking and bursting.
International Journal of Bifurcation and Chaos
,
10
(6)
,
1171
1266
. doi:10.1142/S0218127400000840
Jansen
,
B. H.
, &
Rit
,
V. G.
(1995)
.
Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns.
Biological Cybernetics
,
73
(4)
,
357
366
. https://doi.org/10.1007/BF00199471 doi:10.1007/BF00199471
Jones
,
M. V.
, &
Westbrook
,
G. L.
(1996)
.
The impact of receptor desensitization on fast synaptic transmission.
Trends in Neurosciences
,
19
(3)
,
96
101
. doi:10.1016/S0166-2236(96)80037-3
Kudela
,
P.
,
Franaszczuk
,
P. J.
, &
Bergey
,
G. K.
(2003)
.
Changing excitation and inhibition in simulated neural networks: Effects on induced bursting behavior.
Biological Cybernetics
,
88
(4)
,
276
285
. doi:10.1007/s00422-002-0381-7
Kunze
,
T.
,
Peterson
,
A. D. H.
,
Haueisen
,
J.
, &
Knösche
,
T. R.
(2017)
.
A model of individualized canonical microcircuits supporting cognitive operations.
PLOS One
,
12
(12)
. doi:10.1371/journal.pone.0188003
Lisman
,
J. E.
(1997)
.
Bursts as a unit of neural information: Making unreliable synapses reliable.
Trends in Neurosciences
,
20
(1)
,
38
43
. doi:10.1016/S0166-2236(96)10070-9
Lobb
,
C. J.
(2014)
.
Abnormal bursting as a pathophysiological mechanism in Parkinson's disease.
Basal Ganglia
,
3
(4)
,
187
195
. doi:10.1016/j.baga.2013.11.002
Magill
,
P. J.
,
Bolam
,
J. P.
, &
Bevan
,
M. D.
(2001)
.
Dopamine regulates the impact of the cerebral cortex on the subthalamic nucleus–globus pallidus network.
Neuroscience
,
106
(2)
,
313
330
. doi:10.1016/S0306-4522(01)00281-0
Marder
,
E.
, &
Thirumalai
,
V.
(2002)
.
Cellular, synaptic and network effects of neuromodulation.
Neural Networks
,
15
(4)
,
479
493
. doi:10.1016/S0893-6080(02)00043-6
Montbrió
,
E.
,
Pazó
,
D.
, &
Roxin
,
A.
(
2015
).
Macroscopic description for networks of spiking neurons.
Physical Review X
,
5
(2)
,
021028
. doi:10.1103/Phys-RevX.5.021028
Müller
,
E. J.
,
,
S. J.
,
Kim
,
J. W.
, &
Robinson
,
P. A.
(
2017
).
Unified neural field theory of brain dynamics underlying oscillations in Parkinson's disease and generalized epilepsies.
Journal of Theoretical Biology
,
428
,
132
146
. doi:10.1016/j.jtbi.2017.06.016
Pozo
,
K.
, &
Goda
,
Y.
(
2010
).
Unraveling mechanisms of homeostatic synaptic plasticity.
Neuron
,
66
(3)
,
337
351
. doi:10.1016/j.neuron.2010.04.028
Ratas
,
I.
, &
Pyragas
,
K.
(
2016
).
Macroscopic self-oscillations and aging transition in a network of synaptically coupled quadratic integrate-and-fire neurons.
Physical Review E
,
94
(3)
,
032215
. doi:10.1103/PhysRevE.94.032215
Robinson
,
P. A.
,
Rennie
,
C. J.
, &
Wright
,
J. J.
(1997)
.
Propagation and stability of waves of electrical activity in the cerebral cortex.
Physical Review E
,
56
(1)
,
826
840
. doi:10.1103/PhysRevE.56.826
Schmidt
,
H.
,
Avitabile
,
D.
,
Montbrió
,
E.
, &
Roxin
,
A.
(2018)
.
Network mechanisms underlying the role of oscillations in cognitive tasks.
PLOS Computational Biology
,
14
(9)
,
e1006430
. doi:10.1371/journal.pcbi.1006430
Spiegler
,
A.
,
Knösche
,
T. R.
,
Schwab
,
K.
,
Haueisen
,
J.
, &
Atay
,
F. M.
(2011)
.
Modeling brain resonance phenomena using a neural mass model.
PLOS Computational Biology
,
7
(12)
,
e1002298
. doi:10.1371/journal.pcbi.1002298
Suffczynski
,
P.
,
Kalitzin
,
S.
, &
Lopes Da Silva
,
F. H.
(2004)
.
Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network.
Neuroscience
,
126
(2)
,
467
484
. doi:10.1016/j.neuroscience.2004.03.014
Takeshita
,
D.
,
Sato
,
Y. D.
, &
Bahar
,
S.
(2007)
.
Transitions between multistable states as a model of epileptic seizure dynamics.
Physical Review E
,
75
(5)
,
051925
. doi:10.1103/PhysRevE.75.051925
Tsodyks
,
M.
,
Pawelzik
,
K.
, &
Markram
,
H.
(1998)
.
Neural networks with dynamic synapses.
Neural Computation
,
10
(4)
,
821
835
. doi:10.1162/089976698300017502
Turrigiano
,
G. G.
(2008)
.
The self-tuning neuron: Synaptic scaling of excitatory synapses.
Cell
,
135
(3)
,
422
435
. doi:10.1016/j.cell.2008.10.008
van Vreeswijk
,
C.
, &
Hansel
,
D.
(2001)
.
Patterns of synchrony in neural networks with spike adaptation.
Neural Computation
,
13
(5)
,
959
992
. doi:10.1162/08997660151134280
Virkar
,
Y. S.
,
Shew
,
W. L.
,
Restrepo
,
J. G.
, &
Ott
,
E.
(2016)
.
Feedback control stabilization of critical dynamics via resource transport on multilayer networks: How glia enable learning dynamics in the brain.
Physical Review E
,
94
(4)
,
042310
. doi:10.1103/PhysRevE.94.042310
Wichmann
,
T.
, &
Soares
,
J.
(
2006
).
Neuronal firing before and after burst discharges in the monkey basal ganglia is predictably patterned in the normal state and altered in Parkinsonism.
Journal of Neurophysiology
,
95
(4)
,
2120
2133
. doi:10.1152/jn.01013.2005
Wilson
,
H. R.
, &
Cowan
,
J. D.
(1972)
.
Excitatory and inhibitory interactions in localized populations of model neurons.
Biophysical Journal
,
12
(1)
,
1
24
. doi:10.1016/S0006-3495(72)86068-5
Wong
,
A. Y. C.
,
Graham
,
B. P.
,
Billups
,
B.
, &
Forsythe
,
I. D.
(
2003
).
Distinguishing between presynaptic and postsynaptic mechanisms of short-term depression during action potential trains.
Journal of Neuroscience
,
23
(12)
,
4868
4877
. doi:10.1523/JNEUROSCI.23-12-04868.2003
Zeldenrust
,
F.
, &
,
W. J.
(
2013
).
Modulation of spike and burst rate in a minimal neuronal circuit with feed-forward inhibition.
Neural Networks
,
40
,
1
17
. doi:10.1016/j.neunet.2012.12.008
Zeldenrust
,
F.
,
,
W. J.
, &
Englitz
,
B.
(
2018
).
Neural coding with bursts—current state and future perspectives.
Frontiers in Computational Neuroscience
,
12
. doi:10.3389/fncom.2018.00048
Zierenberg
,
J.
,
Wilting
,
J.
, &
Priesemann
,
V.
(
2018
).
Homeostatic plasticity and external input shape neural network dynamics.
Physical Review X
,
8
(3)
,
031018
. doi:10.1103/PhysRevX.8.031018
Zucker
,
R. S.
, &
Regehr
,
W. G.
(2002)
.
Short-term synaptic plasticity.
Annual Review of Physiology
,
64
(1)
,
355
405
. doi:10.1146/annurev.physiol.64.092501.11454