Understanding the Generation of Network Bursts by Adaptive Oscillatory Neurons

Experimental and numerical studies have revealed that isolated populations of oscillatory neurons can spontaneously synchronize and generate periodic bursts involving the whole network. Such a behavior has notably been observed for cultured neurons in rodent's cortex or hippocampus. We show here that a sufficient condition for this network bursting is the presence of an excitatory population of oscillatory neurons which displays spike-driven adaptation. We provide an analytic model to analyze network bursts generated by coupled adaptive exponential integrate-and-fire neurons. We show that, for strong synaptic coupling, intrinsically tonic spiking neurons evolve to reach a synchronized intermittent bursting state. The presence of inhibitory neurons or plastic synapses can then modulate this dynamics in many ways but is not necessary for its appearance. Thanks to a simple self-consistent equation, our model gives an intuitive and semi-quantitative tool to understand the bursting behavior. Furthermore, it suggests that after-hyperpolarization currents are sufficient to explain bursting termination. Through a thorough mapping between the theoretical parameters and ion-channel properties, we discuss the biological mechanisms that could be involved and the relevance of the explored parameter-space. Such an insight enables us to propose experimentally-testable predictions regarding how blocking fast, medium or slow after-hyperpolarization channels would affect the firing rate and burst duration, as well as the interburst interval.

Introduction

Network bursting is an intermittent collective behavior that occurs spontaneously in neuronal populations. It is characterized by long quiet periods, with almost no spike emission, punctuated by brief periods of intense spiking activity, where the whole network displays high firing rates—most neurons emit at least 2 closely-packed spikes. This particular pattern is then repeated, with varying regularity, over long time intervals.

Such periodic and synchronized activity has been observed as an emergent phenomenon in large neuronal populations, both in brain regions (Meister et al., 1991; Blankenship and Feller, 2009; Rybak et al., 2014) and unperturbed neuronal cultures (Wagenaar et al., 2006; Stegenga et al., 2008; Penn et al., 2016). It has been investigated as a plausible candidate for rhythmogenesis (Ramirez et al., 2004), but also in various disorders such as epilepsy (Derchansky et al., 2008) or Parkinson's disease.

Recent experiments by Penn et al. (2016), studying dissociated neuronal cultures where the chemical environment was precisely controlled, provide evidence that the majority of hippocampal pyramidal neurons are self-sustaining oscillators. These oscillators spontaneously synchronize to give birth to a very regular network bursting phenomenon.

Starting from these results and others (Ramirez et al., 2004; Suresh et al., 2016), we propose here a detailed understanding of the synchronized network bursting dynamics that explains and reproduces other experimental observations (Sipilä et al., 2006; Masquelier and Deco, 2013; Orlandi et al., 2013) of bursting on a variety of different timescales and with inter-burst intervals (IBIs) ranging from less than 1 s up to several minutes. We focus specifically on the characterization of the synchronized attractor and do not consider the transient synchronization process from an asynchronous to a synchronized phase. Indeed, synchronization of pulse-coupled oscillators is a known asymptotic behavior (Somers and Kopell, 1993; Bottani, 1995), which has been shown to lead to bursting in the presence of adaptation (Van Vreeswijk and Hansel, 2001). This was confirmed in all our simulations, regardless of the precise neuronal parameters, as long as they corresponded to adaptive oscillatory neurons. By oscillatory, we mean that a single neurons will spike periodically if uncoupled and considered independently.

Let us insist on the fact that collective bursting, giving rise to “network bursts,” should not be confused with the individual behavior observed at the cellular level for “bursting” or “chattering” neurons. Though they share similar intervals of rapid firing followed by long quiet periods (Connors and Gutnick, 1990; Sipilä et al., 2006), hence the common name, collective bursting can stem from radically different mechanisms and occur on different timescales (see Supplementary S1.3). Here, population-wide bursts are a specific synchronized behavior emerging from the interaction of oscillating, adaptive-spiking neurons which do not display intrinsic bursting behavior when considered independently but only emit single spikes.

The periodic activity of the intrinsically oscillatory neurons present in culture populations and brain regions is assumed to rely on leak currents which affect their excitability (Suresh et al., 2016). More specifically, the persistent, non-inactivating, sodium current I_{Na, p} (Golomb et al., 2006; Penn et al., 2016) and the H-current I_h (Lüthi and McCormick, 1998) are the prime candidates for this intrinsic depolarization. Adaptation, on the other hand refers to the capacity of a neuron to change—here, more precisely, to lower—its excitability in response to continuous or repeated excitation, such as a step-current in electrophysiological experiments, or the intense synaptic input received from its neighbors during a collective burst. Adaptive neurons indeed display periodic firing with a spiking frequency that progressively slows down from its initial high frequency value. The biophysical processes mediating adaptation are thus distinct for the origin of the rhythmic behavior which they modulate, and several potassium currents are considered for this frequency adaptation, like the muscarinic K⁺ (I_M) current or the Ca²⁺ activated K⁺ currents (I_AHP) (Sah and Louise Faber, 2002; Golomb et al., 2006).

We show here that adaptive spiking is a sufficient condition for network bursting, confirming what was suggested by previous studies (Van Vreeswijk and Hansel, 2001; Masquelier and Deco, 2013; Ferguson et al., 2015), and that intrinsically bursting or chattering neurons are not required. Indeed, we focus on the role of adaptation to explain why, as observed in the experiments, the presence of inhibitory neurons is not necessary to obtain regular collective bursting dynamics. Likewise, though short-term synaptic plasticity might play a role in shaping the dynamics (Gritsun et al., 2010; Masquelier and Deco, 2013), we also demonstrate that it is not required to reproduce characteristic timescales of this dynamics.

Methods

We first describe the models used for the different units composing the system (neurons, synapses and network structure). Based on these, we derive an effective model which remains almost completely tractable, so that most of the properties of the collective dynamics can be predicted analytically. This model is based on successive approximations which were validated by numerical experiments: by dividing the cyclic behavior into several subdomains, we isolate regions where the activity can be solved under different approximations. The final solution is thus composed of the concatenation of these different approximations. We also used these simulations to verify and extend the predictions of our analytic equivalent model.

Neuronal Model

We chose the adaptive Exponential Integrate-and-Fire (aEIF) model (Brette and Gerstner, 2005) because of its compromise between simplicity and biological relevance. The dynamical evolution of a neuron is described by two variables, its membrane potential Ṽ, and a slow adaptation current $\tilde{w}$, which are governed by the following equations:

$\begin{array}{l}\text{if }\tilde{V}\le {\tilde{V}}_{peak}\{\begin{array}{l}{\tilde{C}}_m\frac{d\tilde{V}}{d\tilde{t}}=-{\tilde{g}}_L(\tilde{V}-{\tilde{E}}_L)+{\tilde{g}}_L{\tilde{\Delta }}_T{e}^{\frac{\tilde{V}-{\tilde{V}}_{th}}{{\tilde{\Delta }}_T}}\hfill \\ -\tilde{w}+{\tilde{I}}_e+{\tilde{I}}_s\hfill \\ {\tilde{\tau }}_w\frac{d\tilde{w}}{d\tilde{t}}=\tilde{a}(\tilde{V}-{\tilde{E}}_L)-\tilde{w}\hfill \end{array}\\ \text{else if}\tilde{V}>{\tilde{V}}_{peak},\text{then}\{\begin{array}{ccc}\tilde{V}& \leftarrow & {\tilde{V}}_r\\ \tilde{w}& \leftarrow & \tilde{w}+\tilde{b}\end{array}& (1)\end{array}$

where ${\tilde{C}}_m$ is the membrane capacitance, ${\tilde{g}}_L$ is the leak conductance of the neuron, Ẽ_L is its resting potential, $\tilde{\Delta }T$ affects both the slope and the strength of the spiking current, Ṽ_th is the threshold potential, ${\tilde{\tau }}_w$ is the adaptation timescale, ã gives the strength of subthreshold adaptation, $\tilde{b}$ gives the intensity of the spike-triggered adaptation, and Ṽ_r is the reset potential. Ṽ_peak is the spike cutoff for the model. Ĩ_e is an external current to which the neuron can be submitted.

The main difference of this model compared to the well-known integrate-and-fire model is the presence of the second variable, the current $\tilde{w}$, which modulates the neuronal excitability. The synaptic input received by a neuron is represented by the variable Ĩ_s, which is usually time dependent. The neuronal adaptation can be either subthreshold, through the coupling between Ṽ and $\tilde{w}$ via ã, or spike-driven, from the step increments of size $\tilde{b}$ that $\tilde{w}$ undergoes after a spike.

The exponential spike generation present in the aEIF model is more realistic than the hard threshold of the original Integrate-and-Fire model, which leads to unrealistically fast spiking during bursts. The soft threshold of the Izikevich model (Izhikevich, 2007), which also includes adaptation and could have been a possible choice, is similar to that of the aEIF model and would be analytically more tractable. However, it generates a divergence which is not sharp enough, thus leading to overly long interspikes and induces an undesired influence of the cutoff value (V_peak) on the neuronal dynamics (Touboul, 2009). Despite its non-analytic nature, this feature of the aEIF model was therefore critical to capture the inter-spike dynamics inside bursts.

In this study, and in accordance with the experimental observations for several types of pyramidal neurons, we use only neuronal parameters leading to adapting neurons which exhibit periodic spiking. This state is reached through the persistent current Ĩ_e, which drives their progressive depolarization and makes them spike periodically; setting V_r < V_th ensures that the neurons are not intrinsically bursting, as described in Naud et al. (2008).

Contrary to the resting state, where one stable and one unstable fixed point exist (points where both $\dot{Ṽ}$ and $\dot{\tilde{w}}$ are zero), the periodic activity occurs after these two points disappear through a bifurcation, as described in Brette and Gerstner (2005) and Touboul and Brette (2008), when I_e becomes high enough. In this spiking regime, no fixed point is present in phase space, which allows the neuron to depolarize until V_peak before being reset to V_r, thus following a discontinuous limit cycle.

Illustration of the resting and spiking behaviors can be found on Figures S2, S3, while biologically-relevant values of the parameters used for the aEIF model can be found in Table S1, in the Supplementary Material.

For these parameter sets, we have ${\tilde{\tau }}_w\gg {\tilde{\tau }}_m=\frac{{\tilde{C}}_m}{{\tilde{g}}_L}$, as the typical timescales for the continuous variation of $\tilde{w}$ relate to medium and slow after-hyperpolarization, which occur over hundreds of milliseconds (Sah and Louise Faber, 2002).

During the rest of the study, we use the dimensionless version of the model:

$\begin{array}{ll}\{\begin{array}{l}\dot{V}=-(V-E_L)+e^V-w+I_e+I_s\hfill \\ {\tau }_w\dot{w}=a(V-E_L)-w\hfill \end{array}& (2)\end{array}$

Details for the change of variables can be found in the first section of the Supplementary Material, “Neuronal model and parameters.” From then on, all equations involve only dimensionless variables and parameters.

Synaptic Model

The coupling strength between a pre-synaptic neuron j and a post-synaptic neuron i, such that j → i, is represented by the total charge Q_s transmitted from j to i. This charge is passed dynamically through the ion channels of the synapses, which we represent here by an alpha-shaped post-synaptic current (PSC) (Roth and van Rossum, 2009). If neuron j spikes at time t_j, the triggered PSC is felt by i, after a delay d_ji, and is described by:

$\begin{array}{l}{I}_s(t)=\Theta (t-t_j-d_{ji})I_{ji}(t-t_j-d_{ji})\\ =s_{ji}{I}_0·(t-t_j-d_{ji})\Theta (t-t_j-d_{ji})e^{-\frac{t-t_j-d_{ji}}{{\tau }_s}}.& (3)\end{array}$

Where s_ji is the strength of the synaptic connection from j to i, τ_s is the characteristic synaptic time, Θ(x) is the Heaviside step function, such that Θ(x) = 0 is x ≤ 0 and Θ(x) = 1 if x > 1, and $I_0=\frac{1pA}{{\tilde{g}}_L{\tilde{\Delta }}_T}$ is the unit current which we set in this way to be coherent with the conventions of the NEST simulator (Kunkel et al., 2017). As such, the total charge delivered to i reads:

$\begin{array}{ll}{Q}_{s,ji}={\displaystyle {\int }_0^{\infty }{I}_{ji}}(t)dt=s_{ji}{I}_0{\tau }_s^2.& (4)\end{array}$

Network Models

This study is based on two non-spatial random network models: a fully homogeneous network with fixed in-degree which is useful to introduce the equivalent model, and more heterogeneous Gaussian in-degree networks which are supposed to be representative of connectivity in dissociated cultures (Cohen et al., 2010). Both random networks are generated in the same way by drawing a number k_i (in-degree) of incoming connections originating from randomly chosen other neurons in the population. In the case of fixed in-degree networks, the in-degree k_i is fixed and identical for each neuron. For Gaussian random networks, k_i is drawn for each neuron from a Gaussian distribution with mean value $\overline{k}$ and standard deviation σ_k. Note that the fixed-in-degree networks can be seen as the limit case of the Gaussian ones when the variance goes to zero. The out-degree distributions are binomial and identical in both cases. All networks where generated using the graph-tool or igraph backends of the NNGT library.

All transmissions between neurons in the network are subjected to the same delay d and have the same synaptic strength s. which means that the complete dynamical system describing the network is given, for each neuron i, by:

$\begin{array}{ll}\{\begin{array}{l}{\dot{V}}_i=-(V_i-E_L)+e^{V_i}-w_i+I_e+\underset{I_{syn,i}(t)}{\underbrace{{\displaystyle \sum _{j\to i}{\displaystyle \sum _{t_j}s}}{I}_0·(t-t_j-d)\Theta (t-t_j-d)e^{-\frac{t-t_j-d}{{\tau }_s}}}}\hfill \\ {\tau }_w{\dot{w}}_i=a(V_i-E_L)-w_i\hfill \end{array}& (5)\end{array}$

Where {j → i} is the set of neurons j that are presynaptic neurons for i and {t_j} is the set of spike times for neuron j.

Numerical Simulations

All dynamical simulations were performed using the NEST simulator (Kunkel et al., 2017) with the aeif_psc_alpha model implementation, that we developed, and which corresponds to the Equations (2, 3) presented above. Neurons were set to adaptive spiking using the neuronal and synaptic parameters detailed in Table S1 and were connected using static_synapses, i.e., without plasticity, but including a delay d in the spike transmission. Simulations were started from a population of neurons in an asynchronous random state, with their state variable $\tilde{w}$ following a normal distribution of average value 50 pA and standard deviation 10 pA. The runs were performed on networks containing 1,000–100,000 neurons with an average degree of 100, which is the typical value estimated in mature neuronal cultures (Cohen et al., 2010).

Activity Analysis

For each simulation we computed the average firing rate $\nu =\frac{N_s}{T}$, where N_s is the total number of spike and T is the simulation time. This gives us a characteristic timescale t_ν, which would be the average interspike if the spikes were distributed uniformly. Considering d as the transmission delay of action potentials, bursts are identified as uninterrupted sequences of spikes separated by less than min(t_ν/2, 3d); they must also involve at least 20% of the neurons. This analysis was performed using tools from the NNGT library and extra functions available on our GitHub repository.

Equivalent Analytical Model

We derived an equivalent model that describes the system dynamics and predicts the range over which the characteristic frequencies can vary without the need to simulate the network dynamics.The model focuses on the fully synchronized dynamics, for which all neurons behave almost identically. The rationale of the model is most apparent if we first consider the case of a fixed in-degree network. As illustrated on Figure 1, in this case, once the population is synchronized, all neurons receive the same input, that is the contribution of k simultaneous spikes given by the sum of k PSCs. Here, one neuron behaves exactly as any other neuron, thus, ∀i, j, t V_i(t) = V_j(t) = V(t).

FIGURE 1

Figure 1. Schematic representation of the equivalence between a fixed-in-degree network containing N = 6 synchronous neurons, with in-degree k = 2 and connection strength s, and N isolated neurons, with a self-loop connection of strength ks.

This means that the network of N neurons receiving k inputs of strength s is equivalent to N isolated neurons, each one forming a close loop with one autapse—that is, a self-loop—of strength k × s. This simplification is inexact if all neurons do not have the same number of incoming connections, however, as shown in the Results section, this approximation holds very well for homogeneous Gaussian networks and, through a slight modification of the synaptic dynamics, even the behavior of more heterogeneous Gaussian or scale-free networks can be estimated.

Based on this observation, exact for fixed in-degree networks, we propose a model of bursting dynamics for any synchronized network, where we describe the whole population through the behavior of an equivalent neuron, representative of the “average” dynamics. This neuron is subjected to the “average” input received by neurons in the network, and, under this simplified description, Equation (6) is now the same for every neuron in the network, since they are all approximated by this equivalent neuron. As they all receive the same number of spikes (n_s) emitted at the same times {t_j}, j ∈ [1, …, n_s], and from the same number k of neighbors, we obtain:

$\begin{array}{ll}\{\begin{array}{l}\dot{V}=-(V-E_L)+e^V-w+I_e+\underset{I_{syn}(t)}{\underbrace{{\displaystyle \sum _{j}k}s{I}_0·(t-t_j-d)\Theta (t-t_j-d)e^{-\frac{t-t_j-d}{{\tau }_s}}}}\hfill \\ {\tau }_w\dot{w}=a(V-E_L)-w\hfill \end{array}& (6)\end{array}$

This single dynamical system is then solved through several approximations depending on the network state. A typical approximation in the burst, on the interval [t_i, t_i + d] between the emission of a spike and its arrival, consists in linearizing the exponential term when V < V_th. On this interval, I_syn = 0 and since d ≪ τ_w, w can be considered as constant. This leads to an approximate solution for V(t) that we will call V_l(t), with V_l(t_i) = V_r (see also Equations S2, S3 in Supplementary Material):

$\begin{array}{ll}{V}_l(t)=V_r{e}^{-t}+(E_L+I_e-w)\left(1-e^{-t}\right)\text{for}{t}_i\le t<t_i+d.& (7)\end{array}$

For V ∈ [V_th, V_peak], we cannot solve the equation, but know from simulations that this simply leads the neuron to spike with a typical timescale of τ_m = 1.

From these analytic formula, we can then constrain the final solution through a self-consistent equation. The solution of the self-consistent equation will therefore assure that the spikes of one neuron during a burst sustain the burst itself and drive the subsequent ones (self-loop in the equivalent representation of Figure 1), thus shaping a permanent and self-sustained bursting activity, as observed experimentally. Such a solution gives a complete description of the neuron's dynamical properties in time and allows us to obtain all the characteristics of the bursting dynamics.

This equivalent approach is applied here to three different synaptic models (instantaneous, continuous, and alpha-shaped synapses) leading to three transcendental self-consistent equations; details of mathematical developments can be found in the Supplementary Material. Python tools to solve the self-consistent equations and compute the characteristics of the bursting behavior are available on our GitHub repository; they are based on the scipy implementation of Brent's root-finding method.

Exploration of Parameter Space

Thanks to the fast computation of the equivalent model, we were able to compute the dynamical properties for a large number of parameter sets. These results were normalized and analyzed through a Principal Component Analysis algorithm, using the scikit.learn package, in order to obtain the correlation matrix linking the collective dynamical properties to the precise values of the neuronal parameters.

For each parameter set, we first ensure that there is no stable fixed-point in phase-space and that the model predicts a solution, i.e., the existence of bursts with mathematically coherent properties. Secondly, we assess the biological relevance of the solution by (1) ruling out dynamics for which the voltage decreases to values lower than −120 mV during the giant hyperpolarization following a burst; (2) restricting the maximum value of the slow current w to 1,000 pA; (3) preventing cellular bursting for individual neurons by asserting V_r < V_th—this restricts the neurons to single-spike intrinsic behaviors (Naud et al., 2008).

These constraints limit the number of “valid” parameter sets and make the parameters inter-dependent; this leads to a non-trivial parameter/parameter correlation matrix (Figure S1).

Results

As mentioned in the introduction and discussed in the Supplementary Material, synchronization is highly resilient and we focus here solely on the fully synchronized bursting network. We start from individual neurons which are spiking periodically, a behavior that seems to originate from persistent sodium currents like I_{Na, p} or I_h in neuronal cultures (Penn et al., 2016); it is modeled here by a constant input current I_e. When these neurons are coupled, however, their periodic dynamics is drastically modified as they adopt a collective bursting behavior (Borges et al., 2017).

We describe the attractor characterizing the dynamics of the synchronous bursting state. Our key result details the properties of this attractor and shows how they are linked to both the biological parameters of the neurons and the network topology.

The behavior shows features of a relaxation oscillator (see Figure 4A): the current w slowly decreases during the quiescent phase, then rapidly increases during the bursting phase until it reaches a threshold value w*, which determines the burst termination and the start of a new cycle.

The main characteristic which determines the dynamics is the maximum value of the adaptation current, w* reached at the end of a burst. It depends on the neuronal and network parameters, and qualitatively obeys the following equation (details in subsection 5.3 of the Supplementary Material):

$\begin{array}{ll}{w}^{*}\approx E_L+I_e+\frac{V_r{e}^{-d}+\overline{k}{Q}_s}{1-e^{-d}}+C& (8)\end{array}$

where C is a constant. Since the firing rate during the burst is mostly linked to w*, this equation directly shows that higher coupling ($\overline{k}{Q}_s$), higher excitability (E_L), or higher reset voltage (V_r) will increase the bursting intensity. The effect of the transmission delay d is slightly more complex but roughly decreases bursting intensity when increased.

Taking into account finer effects and spike-driven adaptation then leads to more complete equations delivering additional results about the influence of the remaining parameters. These are considered in more details in the Discussion section.

In the following subsections, we describe and explain the bursting dynamics, then discuss the more detailed, self-consistent versions of Equation (9) (complete derivation of these equations can be found in the Supplementary Material). Finally, we describe how our model accounts for the structural heterogeneity that is present in neuronal cultures.

The Attractor, Inner Structure of a Burst

The synchronous attractor is composed of intermittent bursts of activity, as shown in Figure 2 in the (V, w) phase space. During a cycle, the neuron state variables (V, w) do not follow the attractor at constant velocity: the neuron spends much longer on the recovery path (low V) compared to the bursting period (high V)—see Figure 4 to see this trajectory in time.

FIGURE 2

Figure 2. Attractors for three different networks of 1,000 identical neurons with average degree 100. Fixed-in-degree is represented by the blue solid line (spike positions are represented by empty squares and reset positions by full circles). For Gaussian in-degree networks, the logarithm of the number of states per bin—over 200 simulations with 4 cycles each—was used to compensate the non-constant velocity across the whole attractor. The larger attractor, in green, is associated to σ_k = 4; the smaller one, in purple, is for σ_k = 20; both attractors are delimited by a dashed line (limit of a unique visit per bin). Bin size is approximately 0.05 mV along the V-axis and 1 pA along the w-axis.

This attractor is modified by the presence of heterogeneity in the network's topology—quantified by σ_k for Gaussian in-degree networks—which impacts both its duration and regularity. Indeed, heterogeneity noticeably smooths the average behavior and reduces the number of spikes in a burst which goes down from 6 spikes per burst for the fixed in-degree graph, to 3–5 if σ_k = 4, and is roughly reduced to 2 when σ_k = 20. For the fully synchronized fixed in-degree network, all neurons are responding to the exact same input—they receive spikes from the same number of neighbors—hence they are all equivalent to a single average neuron.

As can be seen on Figure 3 for a fixed in-degree network, synchronized bursting of the population consists of a succession of active periods, called bursts¹, separated by long inactive intervals, which we call recovery periods. As can be seen on the inset, the burst displays a strongly ordered inner structure composed of successive synchronized burst slices, which are consistent sets of spikes stemming from a common input.

FIGURE 3

Figure 3. Spike raster of bursting activity for a fixed 100-in-degree network. Inset provides details on the behavior of the neurons during a single burst, with successive synchronized burst slices separated by longer and longer intervals as the adaptation increases.

This inner structure, based on spike events, helps us define several quantities that characterize the dynamics such as the burst and inter-burst durations. However, information about the spike times alone is not sufficient to provide insights regarding the phenomena involved in the burst initiation or termination. Therefore, we will use the time evolution of the neuron's state variables to perform a phase-plane analysis and investigate possible mechanisms for both the bursting and recovery periods.

Neuronal Trajectory, Assessment of the Theoretical Model

From the simulation, we can record the evolution of V and w during the whole dynamics to reconstruct the trajectory of the neuronal state, both in time and in phase-space. Figure 4A represents the time evolution of the equivalent neuron (see Figure 1) during a bursting dynamics on a regular fixed in-degree network and the comparison with the trajectory predicted by the “alpha” equivalent model (see section 6 of the Supplementary Material, “A more detailed model: alpha-shaped synapses”). The close agreement between these trajectories shows that the theoretical model has a good predictive power. Indeed, the most visible discrepancy between the equivalent model and the simulations concerns the precise spike times, as shown in the inset of Figure 4A; however, though the difference can be significant on the intraburst timescale, it is in fact limited to a few milliseconds, which is negligible compared to the duration of a cycle.

FIGURE 4

Figure 4. (A) Time evolution of parameters V and w for the theoretical model (thick black lines and circled numbers) and for a simulation (thin red curve for V, blue for w) on a fixed in-degree graph with k = 100. Two bursts are represented and the intraburst dynamics is presented in the inset, where the spike times predicted by the equivalent model are marked by black dots. The numbered circles mark the main points of the theoretical dynamics, where the behavior changes, as described in the SI. (B) Plot of the attractors in phase space, both for the theoretical model (thick black curve), and through a simulation (thin brown). The three first periods following a burst are denoted by blue arrows: there is first a sharp decrease of V down to its minimum value as it crosses the V-nullcline (dashes); it is followed by a short interval where the neuronal state moves rapidly toward the “recovery curve” (dotted line), which is then followed until the minimum of the V-nullcline and the bursting sequence. The spike trajectory is cut on the figure (marked by empty squares) and the following reset point is marked by a filled circle, as the voltage is set back to V_r and w is increased by b. Neuronal and synaptic parameters are detailed in Table S1, Set 1. The w-nullcline is outside the range of (B).

The dynamics can be understood most easily when looking at w since its behavior can be seen as relaxation oscillations: after a burst (0), the adaptation variable undergoes a quasi-exponential decrease until it reaches its minimum value w_min—passing through points (1) to (4). At this point, the burst starts and w increases rapidly toward a peak value w*—point (5) on Figure 4, which characterizes the trajectory, and will be derived below. Once this maximum value is reached, the neuron stops spiking, the increase of w stops, then the cycle starts again (see Supplementary animation online).

The evolution of V can then be seen as an interplay between the influence of w, I_e, and the synaptic currents in the active period:

• During the burst, each new spike induces a strong depolarization of the membrane, thus leading to another spike—point (4) to (5) on the figure.

• Once w reaches its peak value w*, its influence becomes predominant and prevents the neuron from firing; once the effect of the last spike vanishes, it drives a fast hyperpolarization of the neuron down to point (1).

• After V has reached a quasi-equilibrium value along its nullcline, it instantaneously adapts to the slow decay of w and increases progressively until the trajectory reaches the lowest point of the V-nullcline—point (3). This recovery from the strong hyperpolarization is greatly influenced by I_e.

• At this point, the potential starts increasing more rapidly as the first spike is initiated until the bursting starts again with (4), where the first spike predicted by the equivalent model occurs.

Understanding the Initiation and Termination of a Burst

One of the main interests of this equivalent model is that it provides an intuitive understanding of the mathematical conditions describing the initiation and the termination of bursts. As shown on Figures 4B, 5, the whole existence of the short active period can be understood from the position of the neuronal state in phase space compared to the V-nullcline (curve $\dot{V}=0$), which can be seen as an effective threshold. Indeed, the initiation of the burst simply occurs when w becomes low enough so that the trajectory can “pass under” the V-nullcline; this can be understood easily since the excitability of the neuron increases when w decreases. The lowest value w_min represent the situation where the excitability of the neuron has become so high that it spontaneously emits a spike.

FIGURE 5

Figure 5. Trajectory of a “Dirac burst” in dimensionless phase space; the gray numbers indicate the order of the burst initiation. After a reset, the potential first decreases (leftmost parts of the trajectory) until the spike arrives (brown square), at which point the potential is suddenly shifted to the corresponding brown dot on the rightmost part of the trajectory. The decay before the spike arrival becomes more and more significant as w increases since it contributes negatively to $\dot{V}$. Burst continues until w becomes greater than w*, denoted by the green dot, where the V-nullcline (black line, representing the set of points (V_NV(w), w)) is crossed. Once w_max is reached (circled 0), the burst ends and the recovery period starts.

A key result is then the derivation of a condition for burst termination. We show that the end of the spiking sequence that constitutes a burst is ensured by the intrinsic dynamical properties of single neurons—through adaptation mechanisms—and does not require inhibition nor plastic synapses.

To understand the succession of spikes during the burst and why this spiking process comes to an end, we must introduce a description of the dynamic coupling between the neurons. We first explicit this coupling for two limit cases: firstly instantaneous couplings in perfectly regular fixed in-degree networks, using synapses modeled by Dirac delta functions (called Dirac synapses in the following); secondly, mimicking the effect of highly disordered networks, where synapses release a constant current over the entire burst duration. Thirdly, we consider a more biologically relevant coupling using alpha-shaped synapses, detailed in section 6 of the Supplementary Material, which lies between these two previous limits.

In general, the synaptic coupling I_s between the neurons is time-dependent, which makes the resolution of the system's dynamics (Equation 3) highly complex. As a result the V-nullcline ($\dot{V}=0$) is not generally fixed over a whole cycle. This complicates the threshold condition on w* in the case of the “alpha” synapses. Therefore, the Dirac and continuous synaptic models are more convenient since they enable us to get an insight on the bursting mechanisms through a static representation of the phase diagram during a burst.

Regular Networks and Dirac Synapses

The rationale for the condition of burst termination is most easily understandable in the case of regular networks assuming a coupling in the form of Dirac synapses. Indeed, the arrival of a spike then simply results in a step increment of the post-synaptic neuron's membrane potential:

$\begin{array}{ll}V(t_{sp}^{ + })=V(t_{sp}^{ - })+\overline{k}{Q}_s& (9)\end{array}$

where t_sp is the time at which the spike is delivered to the post-synaptic neuron; $t_{sp}^{-}$, $t_{sp}^{+}$ are respectively the instants immediately before and after spike delivery. Q_s is the total charge delivered by the spike and reflects the coupling strength in the network.

The behavior of the neuron can easily be understood by looking at the situation in phase space on Figure 5. Due to the instantaneous coupling through the Dirac function, there is no finite period of time where the equation for V receives a non-zero input. Consequently, in this limit the V-nullcline remains fixed at all times. Therefore, the condition for the occurrence of a new spike during the burst depends only on the position of $V(t_{sp}^{+})$ compared with the value of the V-nullcline at the same w: V_NV(w). During an interspike of duration T_I(w), w can be considered as constant since τ_m, T_I(w) ≪ τ_w (quasi-static approximation). Hence, either $V(t_{sp}^{+})>V_{NV}(w)$ and a new spike occurs, or $V(t_{sp}^{+})\le V_{NV}(w)$ and the burst terminates.

Developing this condition mathematically leads to the following self-consistent equation:

$\begin{array}{ll}{w}^{*}=E_L+I_e-V_r+\left[{{W}}_{-1}\left(-e^{E_L+I_e-w^{*}}\right)+\overline{k}{Q}_s\right]e^d& (10)\end{array}$

where ${{W}}_{-1}$ is the lower branch of the Lambert W function.

Heterogeneous Networks and Continuous Synapses

For very heterogeneous networks, the broad in-degree distribution leads the neurons to fire at seemingly random times during the bursting period. In the limit where the time distribution of the spikes inside a burst becomes completely uniform, we can approximate it through a window-like synaptic current which is zero during the interburst, then jumps to a finite constant value during the burst.

To obtain an effect equivalent to the spikes described in the previous subsection, devoted the Dirac model, the total charge transmitted during the burst should be the same if an equal number of spikes is emitted. This condition reads, for an average in-degree $\overline{k}$, and a mean synaptic current $I_s^{(c)}$ during the burst,

$\begin{array}{ll}{I}_s^{(c)}{T}_B=n_s\overline{k}{Q}_s.& (11)\end{array}$

where n_s is the number of spikes inside the burst. As described previously, the burst termination occurs when the trajectory crosses the V-nullcline. Figure 6 shows this condition in this heterogeneous limit, i.e., as the input received by the neurons becomes continuous during the burst.

FIGURE 6

Figure 6. Trajectory of a burst in dimensionless phase space for neurons coupled via continuous synapses. Once the first spike occurs (marked by 4), the burst is initiated, i.e., a continuous current $I_s^{(c)}$ is injected into the neurons, thus shifting the resting V-nullcline (dashed curve) upwards (solid black). The neuron spikes until the last shift of b brings w above w*, at w_max, where it encounters the nullcline. This marks the end of the burst and the beginning of the recovery period (circled 0).

Because of the quasi-static hypothesis on w during an interspike, burst termination arises when w goes above the lowest point of the V-nullcline, which occurs for $w^{*}=1+E_L+I_s^{(c)}+I_e$. This is obtained by setting V = 0 in w_NV(V). After a few lines of calculation detailed in the Supplementary Material (Equations S8–S10), we obtain the self-consistent equation:

$\begin{array}{ll}{w}^{*}=w_{min}+b\left[\overline{t_s}(w^{*})-d\right]+\overline{k}{Q}_s.& (12)\end{array}$

where $\overline{t_s}(w^{*})$ is the average interspike interval (ISI) in the burst. As in the previous equations (Equations 8, 10), the critical value of the adaptation current wj at which the burst terminates (1) increases when the coupling strength ($\overline{k}{Q}_s$) increases (2) decreases when the transmission delay d increases. Furthermore, this self-consistent equation also shows the effect of the spike-driven adaptation b which increases the maximum value of the adaptation current that can be reached.

Summary of the Theoretical Description

Once w* has been computed using one of the theoretical models, we can derive all the dynamical properties, starting with:

$\begin{array}{ll}{w}_{max}=\{\begin{array}{ll}\lceil \frac{w^{*}-w_{min}}{b}\rceil b\hfill & \text{for fixed in-degree networks}\hfill \\ w^{*}\hfill & \text{with heterogeneity}\hfill \end{array}& (13)\end{array}$

where ⌈·⌉ denotes the ceiling function. Though the self-consistent equations derived above are less easy to interpret compared to the approximated solution (Equation 9), they allow precise quantitative predictions of the network's dynamics without too much computational cost.

Note that the neurons follow a well-defined and unique attractor, with w changing by discrete steps, only in the case of a fixed in-degree network, where they are all equivalent and synchronous, hence the dual form of Equation (14). In the presence of heterogeneity, the attractor has fuzzy boundaries, as shown by numerical simulations on Figure 2. In this case, the average adaptation over all neurons has a smooth dynamics and w_max is closer to the statistical value at which the neurons stop bursting: w*.

The complete dynamics of the model can be completely captured by the relaxation behavior of w, which displays two phases: one resting period where the adaptation variable decreases until it reaches its lowest value, and an active period where w increases rapidly up to its peak value. The duration of the resting period (interburst interval, or IBI) can be approximated as the sum of the following terms:

T_down     characterizes the time necessary for the neuron to undergo its strong hyperpolarization and reach its lowest membrane potential—from (0) to (1) on Figure 4,

T_R           is the duration of the recovery—from (2) to (3),

T_fs           is the time necessary for the initiation of the first spike which is roughly equivalent to the membrane time constant τ_m—from (3) to (4).

This allows us to obtain the characteristic values of the dynamics (see section 8 of the Supplementary Material, “Resting period”, for detailed calculations):

n_s $=\lceil \frac{w^{*}-w_{min}}{b}\rceil$ (average number of spikes in a burst),

T_B $={\displaystyle \sum _{j=1}^{n_s-1}}{t}_s(w_{min}+jb)$, where t_s(w) is the interspike interval (ISI) for a given value of w,

T_down $=ln\left(\frac{\lambda }{\lambda -V_{max}+E_L+I_e-w_{max}}\right)$,

T_R $=\frac{{\tau }_w-a}{1+a}ln\left(\frac{w^{(2)}-\frac{a}{1+a}{I}_e}{w_{min}-\frac{a}{1+a}{I}_e}\right)$, where w⁽²⁾ is the value of w at point (2),

IBI ≈ T_down + T_R + 1.

Because these results are analytic, thus immediate to compute, this has the significant advantage over simulations that it allows us to quickly predict the properties of the collective dynamics for a large number of parameter sets, i.e., of individual neuron's behaviors.

Evolution of the Properties with Neuronal and Synaptic Parameters

In order to assess the separate influence of the different neuronal parameters on the bursting properties, we used the model to test in a systematic way the influence of the separate variables. As can be seen on Figure 7, this allows to compare the relative influence of any desired set of parameters in a fast and systematic way. Thus, it is a valuable tool to make preliminary explorations in order to prepare for subsequent experimental tests.

FIGURE 7

Figure 7. Correlation matrix for the main characteristics of the bursting dynamics vs. neuronal parameters. 〈ISI〉 is the mean value of the interspike over one burst. Correlations were performed over 2 million randomly-drawn neuronal parameter sets using the predictions of the equivalent model. The experimentally observable features are the IBI, T_B, n_s, and 〈ISI〉.

This matrix allows us to confirm obvious trends, such as the negative influence of the driving current I_e on the IBI, as it tends to quicken the depolarization of the neurons. Likewise, τ_w is almost linearly related to the IBI since it dictates the decay time for w. However, this systematic study also revealed less predictable correlations. Indeed, one of the most interesting features is the quasi-absence of influence of the subthreshold adaptation variable a compared to the spike-driven adaptation (characterized by b and V_r) on the most visible features of the activity, namely the IBI and burst duration.

Correlations for 〈ISI〉 should be treated with care as this value is the average of the interspike interval over a burst. An additional spike (increment in n_s) automatically increases 〈ISI〉 since the new interspike interval will be larger than the previous ones. This is due to the monotonic growth of ISI with w as the burst progresses. Thus, if V_r indeed reduces the interspike duration on the whole, the negative correlation between 〈ISI〉 and b mainly comes from the decrease of n_s as b increases.

Due to the sheer amount of calculation this would require, the theoretical values returned by the equivalent model during this large exploration of parameter-space cannot be verified by simulations in a systematic way. However, the distributions of the bursting characteristics (number of spikes, interburst, and burst duration) are in biologically relevant ranges—see Figure 8. This shows that adaptation alone can lead to network bursts with periods varying from a few tenths of milliseconds to several seconds.

FIGURE 8

Figure 8. Distributions of the burst burst duration T_B and of the IBI (both in ms) for 2 million parameter sets.

Predictive Ability of the Burst Model for Heterogeneous Networks

Our description of periodic bursts predicts the main features of the synchronized bursting rhythmic activity such as its period and firing rate, which are significantly influenced by the presence of heterogeneity in the network's structure, as was already visible in Figure 2. Indeed, as the heterogeneity—namely σ_k—increases, the sharpness of the synchronized burst slices decreases until the spikes contained in the burst become more uniformly distributed; this is clearly visible on Figure 9, which shows the comparison of a burst for two Gaussian in-degree networks with different standard deviations.

FIGURE 9

Figure 9. Rasters of the bursting activity for 2 different Gaussian networks with 1,000 neurons and an average in-degree of 100; each inset details the inner structure of a burst with the successive slices. (A) Homogeneous Gaussian in-degree network with $\overline{k}=100$ and σ_k = 5 leads to well-defined synchronized burst slices inside the bursts. (B) Heterogeneous Gaussian in-degree network with $\overline{k}=100$ and σ_k = 20 leads to fuzzy synchronized burst slices.

Our model is able to take this heterogeneity into account through three synaptic descriptions (Dirac, alpha-shaped, or “continuous”): this allows us to predict the interval in which the bursting properties of most networks should be contained. As shown on Figure 10, they fall in between the Dirac and “continuous-synapses” models. This description successfully accounts for dynamics of networks with low heterogeneity. For high levels of heterogeneity and low synaptic strengths, the model tends to overestimate the synchrony, although prediction of the bursting period remains correct.

FIGURE 10

Figure 10. Variation of the IBI depending on maximum value of the PSC (in pA). Values predicted by the equivalent model are shown in dashed, dot-dashed and dotted lines respectively for the Dirac, alpha and continuous models. Simulated values for a Gaussian network with σ_k = 4 (blue) and σ_k = 20 (green) are superimposed: the main curve represents the average value, while the filled area marks the 5th to 95th percentiles.

Discussion

In all the simulations we performed, we observed that oscillating adaptive spiking neurons synchronize, then start emitting bursts of spikes as the coupling increases.

Our model provides a predictive framework which allows us to determine how this bursting behavior is affected by changes in the individual properties of the neurons.

In the following subsections, we first discuss the validity range of the analytic model. Then, through a thorough mapping of the aEIF parameters to ion channels and biological mechanisms, we make experimentally-testable predictions about the possible influence the main adaptation channels on the bursting behavior. Namely, we suggest how adaptation-channel blockers may affect the dynamics when applied on a bursting neuronal culture.

Validity Range of the Equivalent Model

In order to get meaningful results within the framework of the present model, one must take care to use sets of parameters that lead to adaptive spiking neurons.

More importantly, the conceptual boundaries of the model are reached in the limit of either a very weakly or very strongly coupled neuronal network. For strong coupling the discrepancy between the equivalent model and the simulations mostly occurs because PSCs becomes so intense that a single input can generate several spikes. This can occur in silico but has little biological relevance for adaptive spiking neurons. The weak coupling limit, however, is more revealing since a progressive transition from an asynchronous state to a bursting phase occurs. This transition first involves oscillating firing rates, then synchronous slices containing between one and two spikes, before bursts containing multiple spikes appear. Our equivalent model, designed to describe a fully synchronous bursting dynamics, cannot faithfully capture this smooth transition.

Regarding the network structure, more heterogeneous (e.g., scale-free) networks may also be described by the “continuous-synapse” model on some range of the coupling strength as the qualitative bursting behavior is still present on such networks.

The Influence of Adaptation and Its Biological Origin

Despite its simplicity, the aEIF model takes into account most of the adaptation phenomena involved in biological neurons. Thus, voltage-gated subthreshold adaptation currents, like the muscarinic potassium current I_M (Womble and Moises, 1992) are quantified by the constant a in Equation (3). On the other hand, spike-triggered adaptation, which mostly comes from calcium-gated potassium channels leading to after-hyperpolarization (AHP) phenomena (Sah and Louise Faber, 2002), are quantified by the reset conditions. These calcium-activated currents can be separated into three main types (Sah, 1996; Sah and Louise Faber, 2002; Vogalis et al., 2003) according to their timescales. Over a few milliseconds (1–10 ms) the fast hyperpolarization current fAHP contributes to action potential repolarization, and is thus taken into account by the model through the value of V_r in Equation (1). On an intermediate (“medium”) timescale, the current mAHP has a fast rise-time (less than 10 ms), followed by a decay over 50 to several hundred milliseconds (Storm, 1989); it is modeled by the b step of w after a spike, in Equation (2). Finally, the slow hyperpolarization (Shah and Haylett, 2000; Andrade et al., 2012) current sAHP has a slow rise of 100 ms or more, and an even slower decay over several seconds. It is mostly revealed after a train of action potentials and peaks between 400 and 700 ms. Though this current is not explicitly taken into account by the aEIF model, in the case of bursts, its qualitative effect can be obtained approximately by an increase of τ_w, which lengthens the effect of the potassium current after a burst. One of the limits of the model is its unique timescale for all of the adaptation-related features.

From the exploration of parameter-space, we obtain the correlation matrix of Figure 7, which shows a significant influence of spike-triggered adaptation on the dynamics compared to subthreshold adaptation. A previous study (Augustin et al., 2013) also hinted at the importance of a non-zero b value to obtain low-frequency oscillations. Using the equivalent model, this can be explained easily by the quasi-static hypothesis and the shape of Equation (13). Indeed, the second term of the right-hand side involves the average ISI—which is an increasing function of w*—and the spike-driven increment for the adaptation, b. Thus, the higher the effect of the spike-driven adaptation, the higher w*, which leads to longer interbursts. On the other hand, the quasi-static hypothesis states that the evolution of w is slow compared to that of V, meaning that the subthreshold variations given by a are limited by their slow evolution on a timescale of τ_w.

A significant advantage of this simple description is that the mechanisms proposed by our equivalent model, in light of the correlation matrix on Figure 7, allow us to make qualitative predictions that could be tested to validate it experimentally. Thus, we predict that blocking the voltage-gated adaptation (Stiefel et al., 2008) should have only limited influence on the dynamics through a slight increase in the number of spikes during a burst. On the contrary, blocking one of the calcium-gated channels should lead to drastic changes in the collective behavior:

• Blocking the fAHP channels should be equivalent to increasing V_r, hence increasing the number of spikes in a burst, leading to higher w_max, therefore longer IBI.

• Blocking the mAHP channels through apamin (Sah and Louise Faber, 2002) would be equivalent to lowering the value of b, which should strongly impact the number of spikes inside a burst, therefore its duration. Yet, this should not change w_max significantly, so it should not strongly impact the IBI if the sAHP is significant enough. However, in the case of complete blocking, if the sAHP is not strong enough to compensate, this should lead to the complete disappearance of the bursting behavior.

• Specific blocking of sAHP channels via noradrenaline (Sah and Louise Faber, 2002) should lead to a small increase of the number of spikes during a burst, but would mostly be equivalent to lowering τ_w. In situation where adaptation has the strongest influence over the bursting period, this would lead to a significant decrease of the IBI. This is however unlikely to happen in neuronal cultures, as will be explained below.

These experiments would enable to test the adaptation hypothesis and assess the relative strength of the different processes we described. In fact, some previous studies by Empson and Jefferys (2001) and de Sevilla et al. (2006) have shown results that seem to corroborate the previous predictions, at least regarding the effect of apamin on bursting in slices. However, the first study records only from few individual neurons, and the second uses 4-aminopyridine and Mg²⁺-free medium to trigger the epileptiform activity. To assess the general validity of the proposed mechanisms, one would thus need additional measurements using cultures in physiological conditions, and where each ion-channel would be tested independently while recording larger fractions of the network, either through calcium imaging or MEAs.

Moreover, other features, such as slow modulation of extracellular potassium concentration due to neuronal activity (Bazhenov et al., 2008) have been described in the context of rhythmic activities; these experiments would also help determine whether such phenomena are required as driving forces or only contribute to strengthening existing bursting activities. In our simulations, network bursting is very robust against the following modifications of the system: addition of inhibitory neurons in the network or inclusion of short-time depression in the synaptic dynamics, as shown on Figure 11. It is thus likely that the underlying mechanism we detailed for excitatory synapses can be generalized to these cases. Indeed, the mechanism remains unchanged by plasticity, while adding inhibitory neurons in the population essentially translates into an effective decrease of the excitatory coupling; the latter has been pointed out for percolation in networks of integrate-and-fire neurons.

FIGURE 11

Figure 11. Modification of the original dynamics (left), where only excitatory neurons are present, by the introduction of 20% of non-oscillating, fast-spiking inhibitory neurons (middle), or of plastic synapses exhibiting short-term depression (right). The coherence of the qualitative aspect over three very different systems is remarkable.

Eventually, previous studies (Cohen and Segal, 2011) have hinted at the importance of synaptic fatigue in the burst termination: they showed that the duration of an evoked burst was strongly dependent on the elapsed time since the previous burst, due to the time needed to repopulate the pools of neurotransmitters. What we showed here is also compatible with these results, since they can easily be understood in the framework of our model: a smaller recovery time leads to a higher initial value of the adaptation current, thus shortening the burst duration because the maximum value of w is reached sooner. The effects of synaptic plasticity and adaptation should thus be similar; however, given the timescales reported in the literature, termination could be mostly mediated by adaptation, while the IBI might depend more strongly on synaptic recovery time. In such a case, blocking sAHP as proposed above should not dramatically change the IBI of neuronal cultures.

Conclusion

This study explains the dynamical processes determining synchronous network bursting of a population of oscillating neurons coupled through excitatory synapses. In particular we explain why adaptation is a sufficient condition for collective bursting. We reproduce a large range of biological rhythms with burst frequencies spanning almost 3 orders of magnitude, from a few hundred milliseconds to tens of seconds, in agreement with experimental observations.

Thanks to a phase-space analysis, we are able to propose a mechanism for the initiation and termination of the bursting period related to spike-driven adaptation, which we link to the underlying biological phenomena. The derivation of analytic equivalent models describing the complete bursting dynamics allows us to predict the evolution of the characteristics of the global behavior from the properties of the individual units—neurons and synapses. This enables us to propose a set of experiments which should clarify the role of adaptation currents in network bursting, as well as their relative importance compared to other biological processes such as exhaustion of vesicle pools.

In our description, each new spike in the burst is caused by the previous one, which means that the delay between the emission of a spike and its reception by the post-synaptic neuron has a significant influence on the dynamics. Indeed, we understand intuitively that the longer the delay, the lower the excitability of the neurons when the PSC arrives, since the membrane potential can decay to lower values. This fact, added to the effect of heterogeneity—which tends to reduce the interburst interval—hints at the existence of a limit to the spatial extension which can sustain coherent bursting. Exploring the effect of heterogeneity and spatial embedding (through propagation delays) therefore constitutes a natural continuation of this work. This is certainly necessary to address experimental observations in large cultures, such as the tendency of the activity to initiate in specific regions before it propagates to the rest of the network (Orlandi et al., 2013).

Author Contributions

Analytical study was conducted by MB and TF. TF also performed the numerical simulations and redaction of first draft. Additional simulations and analytical work were performed, respectively, by PM and SM. Final restructuring and adjustments were performed together by SB, TF, SM, and PM.

Funding

This work was granted access to HPC resources from the IDRIS (Orsay) under the allocation 2016057701 made by GENCI (Grand équipement national de calcul intensif).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Elisha Moses, Yaron Penn, and Menahem Segal for their help and the invaluable discussions we had, as well as for sharing their data with us; Hans Ekkehard Plesser for his help and thorough reviews during the implementation and improvements of the numerical models in NEST; Sylvie Thérond for helping with the installation of NEST on the Turing supercomputer.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2018.00041/full#supplementary-material

All codes used in this article are available on our GitHub repositories (https://github.com/SENeC-Initiative/PyNeurActiv and https://github.com/SENeC-Initiative/figures_examples). Note that some of the code might require additional libraries such as NNGT or NEST to run.

Footnotes

1. ^Let us insist once again that the term burst always refers here to the concerted activity of a large fraction of the neuronal population (i.e., a network burst) and should not be confused with single neuron bursting behavior—though they share common characteristics—since they have different origins and population bursts occur on much longer timescales.

References

Andrade, R., Foehring, R. C., and Tzingounis, A. V. (2012). The calcium-activated slow AHP: cutting through the Gordian knot. Front. Cell. Neurosci. 6:47. doi: 10.3389/fncel.2012.00047

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, M., Ladenbauer, J., and Obermayer, K. (2013). How adaptation shapes spike rate oscillations in recurrent neuronal networks. Front. Comput. Neurosci. 7:9. doi: 10.3389/fncom.2013.00009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazhenov, M., Timofeev, I., Fröhlich, F., and Sejnowski, T. J. (2008). Cellular and network mechanisms of electrographic seizures. Drug Discov. Today Dis. Models 5, 45–57. doi: 10.1016/j.ddmod.2008.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Blankenship, A. G., and Feller, M. B. (2009). Mechanisms underlying spontaneous patterned activity in developing neural circuits. Nat. Rev. Neurosci. 11, 18–29. doi: 10.1038/nrn2759

PubMed Abstract | CrossRef Full Text | Google Scholar

Borges, F., Protachevicz, P., Lameu, E., Bonetti, R., Iarosz, K., Caldas, I., et al. (2017). Synchronised firing patterns in a random network of adaptive exponential integrate-and-fire neuron model. Neural Netw. 90, 1–7. doi: 10.1016/j.neunet.2017.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bottani, S. (1995). Pulse-coupled relaxation oscillators: from biological synchronization to self-organized criticality. Phys. Rev. Lett. 74, 4189–4192. doi: 10.1103/PhysRevLett.74.4189

PubMed Abstract | CrossRef Full Text | Google Scholar

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, D., and Segal, M. (2011). Network bursts in hippocampal microcultures are terminated by exhaustion of vesicle pools. J. Neurophysiol. 106, 2314–2321. doi: 10.1152/jn.00969.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cohen, O., Keselman, A., Moses, E., Soriano, J., and Tlusty, T. (2010). Quorum percolation in living neural networks. EPL (Europhys. Lett.) 89, 1–6. doi: 10.1209/0295-5075/89/18008

CrossRef Full Text | Google Scholar

Connors, B. W., and Gutnick, M. J. (1990). Intrinsic firing patterns of diverse neocortical neurons. Trends Neurosci. 13, 99–104. doi: 10.1016/0166-2236(90)90185-D

PubMed Abstract | CrossRef Full Text | Google Scholar

de Sevilla, D. F., Garduno, J., Galvan, E., and Buno, W. (2006). Calcium-activated afterhyperpolarizations regulate synchronization and timing of epileptiform bursts in hippocampal CA3 pyramidal neurons. J. Neurophysiol. 96, 3028–3041. doi: 10.1152/jn.00434.2006

CrossRef Full Text | Google Scholar

Derchansky, M., Jahromi, S. S., Mamani, M., Shin, D. S., Sik, A., and Carlen, P. L. (2008). Transition to seizures in the isolated immature mouse hippocampus: a switch from dominant phasic inhibition to dominant phasic excitation. J. Physiol. 586, 477–494. doi: 10.1113/jphysiol.2007.143065

PubMed Abstract | CrossRef Full Text | Google Scholar

Empson, R. M., and Jefferys, J. G. (2001). Ca²⁺ entry through L-type Ca²⁺ channels helps terminate epileptiform activity by activation of a Ca²⁺ dependent afterhyperpolarisation in hippocampal CA3. Neuroscience 102, 297–306. doi: 10.1016/S0306-4522(00)00494-2

CrossRef Full Text | Google Scholar

Ferguson, K. A., Njap, F., Nicola, W., Skinner, F. K., and Campbell, S. A. (2015). Examining the limits of cellular adaptation bursting mechanisms in biologically-based excitatory networks of the hippocampus. J. Comput. Neurosci. 39, 289–309. doi: 10.1007/s10827-015-0577-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D., Donner, K., Shacham, L., Shlosberg, D., Amitai, Y., Hansel, D., et al. (2006). Contribution of persistent Na⁺ current and M-type K⁺ current to somatic bursting in CA1 pyramidal cells: combined experimental and modeling study. J. Neurophysiol. 96, 1912–1926. doi: 10.1152/jn.00205.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Gritsun, T. A., Le Feber, J., Stegenga, J., and Rutten, W. L. C. (2010). Network bursts in cortical cultures are best simulated using pacemaker neurons and adaptive synapses. Biol. Cybern. 102, 293–310. doi: 10.1007/s00422-010-0366-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience. Cambridge, MA; London, UK: The MIT Press.

Google Scholar

Kunkel, S., Morrison, A., Weidel, P., Eppler, J. M., Sinha, A., Schenck, W., et al. (2017). NEST 2.12.0. doi: 10.5281/zenodo.259534

PubMed Abstract | CrossRef Full Text | Google Scholar

Lüthi, A., and McCormick, D. A. (1998). H-current: properties of a neuronal and network pacemaker. Neuron 21, 9–12. doi: 10.1016/S0896-6273(00)80509-7

PubMed Abstract | CrossRef Full Text

Masquelier, T., and Deco, G. (2013). Network bursting dynamics in excitatory cortical neuron cultures results from the combination of different adaptive mechanism. PLoS ONE 8:e75824. doi: 10.1371/journal.pone.0075824

CrossRef Full Text | Google Scholar

Meister, M., Wong, R. O., Baylor, D. A., and Shatz, C. J. (1991). Synchronous bursts of action potentials in ganglion cells of the developing mammalian retina. Science 252, 939–943. doi: 10.1126/science.2035024

PubMed Abstract | CrossRef Full Text | Google Scholar

Naud, R., Marcille, N., Clopath, C., and Gerstner, W. (2008). Firing patterns in the adaptive exponential integrate-and-fire model. Biol. Cybern. 99, 335–347. doi: 10.1007/s00422-008-0264-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Orlandi, J. G., Soriano, J., Alvarez-Lacalle, E., Teller, S., and Casademunt, J. (2013). Noise focusing and the emergence of coherent activity in neuronal cultures. Nat. Phys. 9, 582–590. doi: 10.1038/nphys2686

CrossRef Full Text | Google Scholar

Penn, Y., Segal, M., and Moses, E. (2016). Network synchronization in hippocampal neurons. Proc. Natl. Acad. Sci. U.S.A. 113, 3341–3346. doi: 10.1073/pnas.1515105113

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramirez, J. M., Tryba, A. K., and Peña, F. (2004). Pacemaker neurons and neuronal networks: an integrative view. Curr. Opin. Neurobiol. 14, 665–674. doi: 10.1016/j.conb.2004.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Roth, A., and van Rossum, M. C. W. (2009). Modeling synapses. Comput. Model. Methods Neurosci. 6, 139–160. doi: 10.7551/mitpress/9780262013277.003.0007

CrossRef Full Text | Google Scholar

Rybak, I. A., Molkov, Y. I., Jasinski, P. E., Shevtsova, N. A., and Smith, J. C. (2014). Rhythmic bursting in the pre-Bötzinger complex: mechanisms and models. Progr. Brain Res. 209, 1–23. doi: 10.1016/B978-0-444-63274-6.00001-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sah, P. (1996). Ca²⁺-activated K⁺ currents in neurones: types, physiological roles and modulation. Trends Neurosci. 19, 150–154. doi: 10.1016/S0166-2236(96)80026-9

CrossRef Full Text | Google Scholar

Sah, P., and Louise Faber, E. S. (2002). Channels underlying neuronal calcium-activated potassium currents. Progr. Neurobiol. 66, 345–353. doi: 10.1016/S0301-0082(02)00004-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, M., and Haylett, D. G. (2000). Ca²⁺ channels involved in the generation of the slow afterhyperpolarization in cultured rat hippocampal pyramidal neurons. J. neurophysiol. 83, 2554–2561. doi: 10.1152/jn.2000.83.5.2554

CrossRef Full Text | Google Scholar

Sipilä, S. T., Huttu, K., Voipio, J., and Kaila, K. (2006). Intrinsic bursting of immature CA3 pyramidal neurons and consequent giant depolarizing potentials are driven by a persistent Na⁺ current and terminated by a slow Ca²⁺ -activated K⁺ current. Eur. J. Neurosci. 23, 2330–2338. doi: 10.1111/j.1460-9568.2006.04757.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Somers, D., and Kopell, N. (1993). Rapid synchronization through fast threshold modulation. Biol. Cybern. 68, 393–407. doi: 10.1007/BF00198772

PubMed Abstract | CrossRef Full Text | Google Scholar

Stegenga, J., Le Feber, J., Marani, E., and Rutten, W. L. C. (2008). Analysis of cultured neuronal networks using intraburst firing characteristics. IEEE Trans. Biomed. Eng. 55, 1382–1390. doi: 10.1109/TBME.2007.913987

PubMed Abstract | CrossRef Full Text | Google Scholar

Stiefel, K. M., Gutkin, B. S., and Sejnowski, T. J. (2008). Cholinergic neuromodulation changes phase response curve shape and type in cortical pyramidal neurons. PLoS ONE 3:e3947. doi: 10.1371/journal.pone.0003947

PubMed Abstract | CrossRef Full Text | Google Scholar

Storm, J. F. (1989). An after-hyperpolarization of medium duration in rat hippocampal pyramidal cells. J. Physiol. 409, 171–190. doi: 10.1113/jphysiol.1989.sp017491

PubMed Abstract | CrossRef Full Text | Google Scholar

Suresh, J., Radojicic, M., Pesce, L. L., Bhansali, A., Wang, J., Tryba, A. K., et al. (2016). Network burst activity in hippocampal neuronal cultures: the role of synaptic and intrinsic currents. J. Neurophysiol. 115, 3073–3089. doi: 10.1152/jn.00995.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Touboul, J. (2009). Importance of the cutoff value in the quadratic adaptive integrate-and-fire model. Neural comput. 21, 2114–2122. doi: 10.1162/neco.2009.09-08-853

PubMed Abstract | CrossRef Full Text | Google Scholar

Touboul, J., and Brette, R. (2008). Dynamics and bifurcations of the adaptive exponential integrate-and-fire model. Biol. Cybern. 99, 319–334. doi: 10.1007/s00422-008-0267-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Vreeswijk, C., and Hansel, D. (2001). Patterns of synchrony in neural networks with spike adaptation. Neural Comput. 13, 959–992. doi: 10.1162/08997660151134280

PubMed Abstract | CrossRef Full Text | Google Scholar

Vogalis, F., Storm, J. F., and Lancaster, B. (2003). SK channels and the varieties of slow after-hyperpolarizations in neurons. Eur. J. Neurosci. 18, 3155–3166. doi: 10.1111/j.1460-9568.2003.03040.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagenaar, D. A., Nadasdy, Z., and Potter, S. M. (2006). Persistent dynamic attractors in activity patterns of cultured neuronal networks. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 73:051907. doi: 10.1103/PhysRevE.73.051907

PubMed Abstract | CrossRef Full Text | Google Scholar

Womble, M. D., and Moises, H. C. (1992). Muscarinic inhibition of M-current and a potassium leak conductance in neurones of the rat basolateral amygdala. J. Physiol. 457, 93–114. doi: 10.1113/jphysiol.1992.sp019366

PubMed Abstract | CrossRef Full Text | Google Scholar