A Mean-Field Description of Bursting Dynamics in Spiking Neural Networks with Short-Term Adaptation

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.

The aim of this article is to provide and validate mean-field descriptions of collective bursting emerging from the dynamic interaction of short-term adaptation mechanisms and recurrent excitation in populations of coupled spiking neurons. Hence, our work closes the gap between investigations of population bursting in spiking neural networks and neural mass models. First efforts into this direction were made for the special case of SFA in a network of coupled linear integrate-and-fire neurons, employing the Fokker-Planck formalism and an adiabatic approximation given long SFA timescales (Gigante, Mattia, & Giudice, 2007). Analyzing this mean-field description, Gigante et al. (2007) were able to identify different types of collective bursting. In this article, we show under which conditions bursting can emerge as a collective phenomenon from different short-term adaptation (STA) mechanisms within a population of coupled spiking neurons. To this end, we consider multiplicative and additive adaptation as generalized descriptions of the vast number of STA mechanisms that have been reported to affect neural excitability. We incorporate these STA mechanisms in a globally coupled population of quadratic integrate-and-fire (QIF) neurons and derive mean-field descriptions of the macroscopic dynamics following the approach by Montbrió, Pazó, and Roxin (2015). Using bifurcation analysis, we identify states of collective bursting as well as the boundary conditions for such bursting to occur. Our results show how changes in either the average input or the strength of short-term adaptation within a neural population can drive the population into and out of collective bursting. We find that short-term adaptation gives rise to bistable regimes of concurrent bursting and nonbursting states, and we show how different short-term adaptation mechanisms can have very similar effects on the population dynamics. Finally, we perform a finite size analysis in which we examine the effects of network size and coupling probability of the spiking neural network on the correspondence between this microscopic network description and the mean-field model.

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.

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).

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).

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.

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.

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.

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).

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.

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.

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.

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.

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).

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.

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.