Macroscopic bursting in physiological networks: node or network property?

Activity pattern modalities of neuronal ensembles are determined by node properties as well as network structure. For many purposes, it is of interest to be able to relate activity patterns to either node properties or to network properties (or to a combination of both). When in physiological neural networks we observe bursting on a coarse-grained time and space scale, a proper decision on whether bursts are the consequence of individual neurons with an inherent bursting property or whether we are dealing with a genuine network effect has generally not been possible because of the noise in these systems. Here, by linking different orders of time and space scales, we provide a simple coarse-grained criterion for deciding this question.


Introduction
Neuronal bursting activity is a ubiquitous physiological state described by a precipitate train of spikes followed by a quiescent period. The phenomenon can be observed on the node or on the network level, and it can be produced by neurons bursting by their own virtue ('inherent bursting') or by neurons that individually respond with regular spiking but exhibit a bursting behavior when embedded as a node into a network (possibly conditional on a particular input to the network that drives the latter into a particular functional mode). Bursting has several distinct functional roles. Within the neocortex's layer IV, bursting activity emerges as a collective phenomenon [1], enabling the amplification of weak thalamic input into the network. In neuronal embryonic cultures, we generally observe bursting activity after a few days of implementation, when the neuronal network starts to develop its structure. In this case, bursting is the fingerprint of the search by the network for its optimal configuration, which is indicated by an avalanche structure of the firing events (see figure 1), the size of which has a power-law characteristic [2][3][4][5]. The most obvious instance of a functional role of bursting is the neuronal contact to muscles, where bursting is required for muscle contraction [6]. Besides that, frequency features of bursting neurons have also been shown to be important due to their resonant properties, which allow them to transmit reliably selective information between neuronal circuits [7,8]. Not least, bursting or synchronized activity hallmarks a number of important conditions in human health, as it has been shown to be closely related, e.g., to the emergence of epilepsy [9] and migraine [10] and to pacemaker function [11].
The experimental situation that we focus on in our investigation is a coarse-grained, macroscopic one, where we have no microscopic access to individual neuronal spiking or where we have too many elements to deal with on an individual basis. We will show that the two main alternatives that lead to bursting produce different effects on the network level (that may have physiological relevance), and we will provide a mixed qualitative and quantitative analysis of the difference between the two situations. As bursting is an intrinsically nonlinear phenomenon, we compare two extreme cases of nonlinear models: the weak nonlinear coupling of linear phase oscillators (Kuramoto case [12,13]) and the coupling of intrinsically nonlinear oscillators ('Rulkov neurons' [14]). Because of their distinguished position among the extant models of neuronal dynamics, both have been abundantly used for modeling neurons and neuronal networks. Authoritative surveys and examples of their potential use for these tasks are provided, e.g., in [15,16] for the Kuramoto model and in [17,18] for the Rulkov model. Due to their prominence, it is, unfortunately, impossible to do justice to the vast literature available. As Rulkov's model recovers essentially all behaviors of neuronal firing observed in physiology (even regular firing), we may see it as the generic nonlinear model covering the whole of the modeling space, beyond simple phase oscillators onwards to the strongly nonlinear behavior of bursting neurons.
Bursting is characterized in all cases by two time scales: a fast one responsible for individual spiking and a slow one responsible for bursting. From a macroscopic and large time-scale point of view, a regular interburst interval between two successive bursts can be considered as the correspondence to one complete oscillation of a regular neuron; i.e., the burst is considered as a single event. The interburst frequency will then be defined by the number of bursts per unit time. In this sense, a phase θ can be associated with the angular position of an equivalent rotating oscillator, and the interburst angular frequency ω is the average of the angular frequencies evaluated over some time series. When coupled, inherently bursting neurons are able to show regular firing and phase synchronization [14]. Based on this, the emergence of neuronal phase synchronization can be seen as equivalent to the mechanism of phase synchronization for coupled Kuramoto oscillators. Regarding phase dynamics, an explicit mapping between the Rulkov model within a given regime of firing, and the Kuramoto model was recently developed [19]. We will see that, surprisingly, this no longer holds if we, instead, consider frequency as the observable; the frequency ω of coupled bursting neurons and the frequency of Kuramoto oscillators depend distinctively on the coupling strength between the neurons.
In contrast to synchronized coupled Kuramoto oscillators, where the coarse-grained frequency (that often displays burst-like characteristics) does not vary with coupling strength (figure 2(a)), for synchronized coupled inherently bursting Rulkov neurons, the mean interburst frequency (MIF) decreases if the coupling strength ε is increased. This unexpected response of inherently bursting neurons is corroborated by physiological observations of coupled pyloric dilator neurons from the lobster stomatogastric ganglion. (For these neurons, it is well known that when they are coupled with artificial electrical synapses, their MIF changes with the coupling strength [20].) One aim of this paper is to explain this unexpected behavior and to exhibit that for the modeling of certain collective neuronal phenomena, the choice of phase oscillators for node dynamics may be insufficient. To work these points out, we first closely follow and then extend to some point the original analysis provided by Rulkov [14].  ω ω hosting a unimodal symmetric frequency distribution.

Frequency dependence of synchronizing Kuramoto neurons
We first recollect how the emergent macroscopic frequency depends on the coupling for the prominent Kuramoto oscillator model of neurons (i.e., weakly coupled linear phase oscillators). It is well known [13] that increased coupling among the oscillators leads to the emergence of synchronized phase behavior. For our investigations, the frequency of oscillation i ω of the oscillator i θ , i = 1,…, N, where N is the network size, will be chosen randomly from a symmetric and unimodal probability density function g ( ) ω , and the coupling that we shall consider is global: (˙i θ is the temporal phase evolution of the ith neuron, and ε is the global coupling strength.) For the simplest case of the coupling of two oscillators (for the earliest experiments regarding biological neurons see [21,22]), we may perform the transformation and consider the condition to frequency synchronization˙0 φ = . Frequency synchronization will be achieved for for all values of the coupling strength c ε ε > [23,24]. Although for a larger number of oscillators the analytics become increasingly difficult [25], it is known that the frequency average of the coupled Kuramoto oscillators continues to be the average frequency of the uncoupled Kuramoto oscillators, where this property is independent of network topology and network size. This invariance of the mean frequency of Kuramoto oscillators is, however, for most modeling of physiological systems, an unrealistic limitation.

Burst-frequency dependence of Rulkov neurons
Most physiological ensembles of neurons that undergo a synchronization process also contain neurons with inherent bursting behavior. (As we will see, neurons may also change from regular to bursting behavior, depending on physiological conditions.) A generic and convenient framework of neurons with different physiological responses is given by Rulkov's two-dimensional (2D) map [14] based on a fast x-variable and a slow y-variable, n n n 1 2 . Through a decomposition of the fast from the slow time-scale, it is possible to arrive from the 2D system at a 1D system, where the slow variable is replaced by the parameter γ as n n 1 2 This dimensional reduction simplifies the bifurcation analysis [14]. As we shall see, the behavior of the 2D system is well described by the simplified system. The gray region in figure 4 shows the accessible parameter space for the 2D Rulkov model. The reduced Rulkov model hosts three important dynamical features: a saddlenode bifurcation, a flip bifurcation, and a crisis [26]. The saddle-node bifurcation, which consists in a collision between a stable and an unstable fixed point, occurs when the following four conditions are satisfied: [27]. From equation (4) for F x ( , ) 0 0 γ , the α and γ that satisfy the four conditions are given by the relation  [28]. In this case α and γ must obey the relation which is represented in figure 4 by a full red line. A crisis is a sudden change in the chaotic attractor (here: its disappearance). This happens if the maximum of function F(x) is mapped into a stable fixed point [29]. This condition is satisfied if the relation  holds, which is represented in figure 4 by the full green line. From the bifurcations of the reduced Rulkov model, we may understand how the individual dynamics separates into different regimes. Whereas for 2.0 α ≲ the xand y trajectories -will be attracted to a fixed point (region I in figure 4), for 2.0 2.58 α ≲ ≲ the y-variable will oscillate between two saddle-node bifurcation points, and the x-variable will be confined to a periodic motion (region II in figure 4). For 2.58 4.0 α ≲ < we have coexistence of the saddle-node bifurcation and of the flip bifurcation points, which produces the triangle bursting. In this case, when the y-variable reaches the largest saddle-node bifurcation point ( max γ ), the x-variable starts to oscillate rapidly as the y-variable decreases towards the flip bifurcation point ( fp γ ), and the amplitude of the oscillation decreases. When the fast oscillations disappear, the bursting terminates, and the trajectory is attracted to the fixed point. After the orbit reaches the stable fixed point, the y-variable slowly increases towards max γ , which completes the loop (region III in figure 4). The square bursting happens if 4.0 4.62 α < ≲ : in this case the y-variable grows up to max γ . At this point, a chaotic attractor emerges, and the y-variable decreases towards the crisis point ( cs γ ). At this point the chaotic attractor vanishes, and, consequently, bursting is terminated. After that, the orbit is attracted by the stable fixed point, and the y-variable increases its value towards max γ again (region IV in figure 4). Already in Rulkov's original paper [14], a hint can be found that in some area of the parameter space, ω might be basically inversely proportional to the distance fp max γ γ | − |. However, no indication was given as to whether this holds generally, or for what network topologies this would be the case. We infer from figure 4 that if the argument holds, then the strength of ε may have a noticeable impact on the burst frequency ω. To pinpoint this, we numerically simulated different network topologies (for diffusive coupling, small-world networks, and scale-free networks) that all exhibit essentially identical behavior. In this paper, though, we restrict ourselves to the case of globally coupled networks. A burst starts whenever y n max γ = .
Defining the oscillation period as the time between the beginning of two successive bursts, we obtain a corresponding phase (φ) as k n n n n n n n As is shown in figure 5(a), in the triangle-bursting regime, the distance between the flip bifurcation point ( fp γ ) and max γ increases linearly in α. ω, being inversely proportional to the distance between these two points, fp max 1 ω γ γ ∝ | − | − , thus decreases (figures 5(b) and (c)). If we increase α within the square-bursting regime, the distance between the crisis bifurcation point cs γ and max γ decreases. The mean interburst frequency thus increases until, towards the end of the square-bursting regime, the proximity to another regime of neuronal activity at 8 3 3 4.62 α = ≈ , and the small distance between the bifurcation points causes other types of modulations that for the interburst frequency become relevant, forcing MIF to decay again.
We proceed from the single-neuron case to the collective behavior by investigating the behavior of a globally coupled set of Rulkov neurons where ε is the coupling strength, N is the network size, and i = 1,…, N. Also for synchronized Rulkov neurons, a reduction of the 2D model to a 1D model is possible, so that We may now use the same arguments that we have used for isolated Rulkov neurons. The critical points for global coupling corresponding to equations (5)-(7) are given by (see [26] for more details) 18 (1 ) 2 (2 6(1 ) ) 3(1 )  (dashed lines), there is a displacement of the bifurcation lines and, as consequence, a change in the MIF. It is now evident that increased coupling increases the distance between the bifurcation points ( figure 6(a)). For triangle bursting, MIF decreases almost linearly with the coupling and decreases for square bursting with an approximately quadratic dependence ( figure 6(b)). While the dependence between frequency and bifurcation point distances decreases in an approximately linear manner for triangle bursting, for square bursting the decrease is nonlinear ( figure 6(c)).
These results demonstrate that, indeed, inside a given neuronal activity regime (triangle bursting, square bursting), the increase of the distance between the bifurcation points in consequence of an increase of ε (or of α), decreases the interburst frequency ω. This dependence is observed at the neuronal level ( figure 5(c)), as well as if the neurons are coupled ( figure 6(c)). If the increase (or decrease) of the coupling strength is too big, there will be additional, non-continuous changes in neuronal activity. Abrupt changes can be seen upon a change between activity regimes, e.g., when, upon the increase of α, we change from triangle bursting into square bursting ( figure 5(b)), or in dependence of the coupling strength.
The difference between the simple dependence on the coupling strength ε by the Kuramoto model versus the essentially strictly monotonous dependence interrupted by events of bifurcations exhibited by Rulkov neurons should therefore be seen as the hallmark of a preponderance of inherently versus non-inherently bursting neurons.

Conclusions
Despite the differences between the Rulkov and Kuramoto neuron models, in the periodic regime 2.0 2.58 α ≲ ≲ the coupling effect in the spike frequency can be neglected, and the Kuramoto model can be used to model phase and frequency synchronization. In the bursting regime of neurons, 2.58 α ≳ , the differences between the models matter, as for inherently bursting neurons, consistently a decrease of MIF with the coupling strength was observed for global, small-world, and scale-free topologies of sizes varying from N = 100 up to N = 10000. Related observations made earlier for different systems ( [30] for diffusive coupling, [31] for smallworld networks, and [32] for scale-free networks) corroborate the interpretation that our observation deals with a general feature of inherently bursting neurons. We have focused on the MIF of Rulkov's model of bursting neurons. This should, however, not be considered a particular modeling case, but rather as a generic framework that reflects dynamical properties essential for both regular and bursting behavior. Neurons that can become inherently bursting may therefore be the more typical case than neurons that do not offer this possibility. The dependence of the MIF on the coupling exhibited by these neurons (figure 6) then would (at least in matters of frequency) prohibit a modeling by simple phase oscillators, which still is the predominant approach. The apparent independence of the observed phenomenon with respect to network topology makes it a rare example of a nontrivial invariant of network topology.
In view of the exhibited different origins that can underlie bursting, we suggest that our observation will be helpful on different levels of physiological coarse-graining. Particularly on higher levels of abstraction of the representation of physiology (i.e., working with higher levels of modularity), it will be to a lesser degree evident whether for modeling subsystems, a phase oscillator model suffices, or whether a bursting model is necessary. Our results provide a simple experimentally accessible indicator for deciding this question. It is conceivable that for neuronal networks with well-controlled architecture, experimental procedures (e.g., the addition of serotonin or related substances) could be developed to mimic an increase or decrease of the coupling strength ε among the nodes of a biological neural network, obtaining in this way the desired information about the dominant nature of the nodes.
Traditionally, most emphasis in the analysis of complex physiological networks has been dedicated to the microscale and the macroscale. Our results, however, put a warning sign against too straightforward an extrapolation from a microscopic level to the macroscopic scale: mesoscale effects, such as the pattern of bursting investigated here from this angle, can generate quite unexpected effects that can cause these extrapolations to fail. This implies that increased efforts on the mesoscale will be necessary to better understand physiological neuronal systems that span different levels of hierarchical organization.