Clique of Functional Hubs Orchestrates Population Bursts in Developmentally Regulated Neural Networks

It has recently been discovered that single neuron stimulation can impact network dynamics in immature and adult neuronal circuits. Here we report a novel mechanism which can explain in neuronal circuits, at an early stage of development, the peculiar role played by a few specific neurons in promoting/arresting the population activity. For this purpose, we consider a standard neuronal network model, with short-term synaptic plasticity, whose population activity is characterized by bursting behavior. The addition of developmentally inspired constraints and correlations in the distribution of the neuronal connectivities and excitabilities leads to the emergence of functional hub neurons, whose stimulation/deletion is critical for the network activity. Functional hubs form a clique, where a precise sequential activation of the neurons is essential to ignite collective events without any need for a specific topological architecture. Unsupervised time-lagged firings of supra-threshold cells, in connection with coordinated entrainments of near-threshold neurons, are the key ingredients to orchestrate population activity.


Introduction
There is increasing experimental evidence that single neuron firing can impact brain circuits dynamics [1]. It has been shown that a single pyramidal cell can trigger whisker deflection [2], drive sensory perception [3] and modify brain states [4]. Similarly, a single GABAergic hub cell can affect collective activity within the developing hippocampal circuitries [5]. In vivo cortical studies have shown that a single extra action potential (AP) can generate a few dozens extra spikes in its postsynaptic targets [6]. Furthermore, a burst of APs, evoked in a pyramidal cell, can propagate through the network activating locally a high fraction of Somatostatine GABAergic cells (a subset of inhibitory neurons) and a few excitatory cells [7]. The capability of single neurons to evoke sparse [6] and network-wide neuronal events [1,4,5] in brain circuits can be interpreted within the framework of the single neuron doctrine, firstly postulated on sensorial perception by Barlow in 1972 [8]. According to this doctrine, the spiking of a single neuron in a network has a high functional relevance being able to code very specifically for high level features of abstraction such as concepts. Face selective cells [9] are a typical example of sparse object representation in the brain and of putative grandmother cells [10]. The sensitivity of neuronal networks to small perturbations, such as those introduced by the firing of a single cell, can also find an explanation within the self-organized criticality (SOC) framework [11]. In the last decade, SOC has widely been proposed as the mechanism underlying power law distributions, with characteristic exponents, featuring the size and duration of population events. These distributions have been measured in-vivo and in-vitro experiments on neuronal networks from invertebrates, rodents, monkeys and humans [12][13][14][15]. The hypothesis underlying the SOC interpretation is that neuronal networks self-organize into a critical state where responses, over temporal and spatial scales of any size (the so-called ''avalanches''), can be triggered by small perturbations. Despite the theoretical frameworks above introduced, one of the main open question is how and why only specific neurons can affect the global network dynamics as observed in [2][3][4][5]. Two main approaches can be foreseen: a ''structural-functional'' approach [16][17][18][19][20], where the specific topology of the network and the connectivity pattern of the cells are responsible for the relevance of the single neuron or a ''dynamical'' approach, where the single neuron becomes relevant due to the nonlinear evolution of neuronal excitability and synaptic connectivity in the network [21,22].
A recent computational study on the synchronization properties of a specific neural circuit [23], has pointed out that the level of burst synchrony is a function of both the network topology and the intrinsic dynamics of peculiar neurons, which have a central location in the network graph. This led the authors to conclude that in realistic neuronal systems the choice of a specific topology is not sufficient to induce an unequivocal dynamical behavior in network activity. To further deepen the comprehension of the interplay among cell excitability and synaptic connectivity in promoting network burst synchrony, in this paper we study the effect of single neurons perturbations on the collective dynamics of a network of leaky-integrate-and-fire neurons with short-term synaptic-plasticity [24]. The relevance of these network models for neuroscience have been demonstrated in many contexts ranging from the modelization of working memory [25] to the possibility to perform computation by ensemble synchronization [26]. Although these models have extensively been studied for their capability to generate spontaneous population bursting, little is known about the influence of single cell perturbations on their global dynamics [24].
In order to analyze the population dynamics in a neural circuit at the initial stage of its development, when both mature and young cells are simultaneously present, we consider a random diluted network presenting developmentally inspired correlations between neuronal excitability and connectivity. The presence of these correlations can render the network sensitive to single neuron perturbation of a few peculiar neurons. The coherent activity of the network can be even arrested by removing or stimulating any of these neurons, which are functional hubs arranged in a clique regulating the neuronal bursting. We show that the level of available synaptic resources influences the reciprocal firing times of the synaptically connected neurons of the clique. However, the fundamental mechanism responsible for the burst triggering relies on an unsupervised process leading to a precise firing sequence between the neurons which are not structurally connected. Furthermore, frequency locking of the same neurons led, counter-intuitively, to anti-resonances [27,28], inducing reduced bursting activity or even complete silence in the circuit. Notably, although obtained in a developmentally regulated framework, these results can also be extended to a more general context where the effective connectivity and excitability of the neurons are dynamically regulated by the different states of brain processing.

Results
In this paper we intend to mimic an immature neuronal network frozen at a certain stage of its initial development, similar to the one examined in the experimental work on developmental hippocampal circuits [5] which inspired this work. At early postnatal stages, the main features characterizing such networks are the excitatory action of GABAergic transmission (which is instead the most common inhibitory source in mature circuits) and the presence of synchronized network events, as largely documented in central and peripheral nervous circuits [29]. According to that, we consider a network model composed of only excitatory neurons and displaying bursting activity.
In particular, we considered a directed random network made of N leaky-integrate-and-fire (LIF) neurons [30,31] interacting via excitatory synapses and regulated by short-term-synaptic-plasticity (see Methods for more details), similarly to the model introduced by Tsodyks-Uziel-Markram (TUM) [24]. As previously shown in [24,26,32], these networks exhibit a dynamical behavior characterized by an alternance of short periods of quasi-synchronous firing (population bursts, PBs) and long time intervals of asynchronous firing. Notably, the presence of short-term-synaptic-plasticity is the crucial ingredient to observe PBs, even without an inhibitory population [24,26,32]. Therefore, the TUM model with excitatory synapses can be considered as a minimal model to mimic the experimentally described stereotypical/characteristic condition of developing neuronal networks [33].
Furthermore, in developing networks, both mature and young neurons are present at the same time, and this feature is reflected in the variability of the structural connectivities and of the intrinsic excitabilities. Experimental observations indicate that younger cells have a more pronounced excitability [34,35], while mature cells exhibit a higher number of synaptic inputs [5,36]. Thus suggesting that the number of afferent and efferent synaptic connections [5,36,37] as well as their level of hyperpolarization [38] are positively correlated with the maturation stage of the cells. The gradient of excitability -with younger neurons more excitable than older ones -could be explained by a gradient in the excitatory action of GABAergic transmission, i.e. older neurons receive a less depolarizing action by GABAergic input [33].
The presence at the same time of younger and older neurons can be modeled by considering correlations among the in-degree and out-degree of each cell as well as among their intrinsic excitability and connectivity. In particular, in the attempt to find the network organization which is more sensitive to single neuron perturbations, we compare the dynamics of networks where none, one or more of the following correlations have been embedded (for more details see Methods and Supplementary Information): N setup T1: positive correlation between the in-degree and outdegree of each neuron; N setup T2: negative correlation between the intrinsic neuronal excitability and the total connectivity (in-degree plus outdegree); N setup T3: positive correlation between the intrinsic neuronal excitability and the total connectivity (in-degree plus outdegree).
Correlated networks with all possible combinations of the setups T1-T3 have been examined. However, the paper is mainly devoted to the comparison of the properties of the network with

Author Summary
To which extent a single neuron can influence brain circuits/networks dynamics? Why only a few neurons display such a strong power? These open questions are inspired by recent experimental observations in developing and adult neuronal circuits, as well as by classical debates within the framework of the single neuron doctrine. In this work we identify and present a mechanism which can explain in neuronal circuits, at some early stage of their development, how and why only a few specific neurons can exhibit such power. For this purpose, we consider a standard neuronal network model whose population activity is characterized by bursting behavior. The introduction of a distribution of correlated neuronal excitabilities and degrees, inspired by the simultaneous presence of younger and older neurons in the network, leads to the emergence of functional hub neurons. These critical cells, whenever perturbed, are capable of suppressing network synchronization. Notably, we show that their strong influence on the population dynamics is not related to their structural properties, but to their operational and structural integration into a clique. These results highlight how network-wide effects can be induced by single neurons without any need for a specific topological architecture.
correlations of type T1 and T2 (as displayed in Fig. S1) with the completely uncorrelated one, which is a directed Erdös-Rényi graph (see Figs. S3A, S3B). In order to test the possible influence of hub neurons on the network dynamics, also few structural hubs have been added to the network whenever correlations of type T1 were embedded (see Methods and Fig. S1 for more details).
It is important to stress that correlations of type T1 and T2 have a justification on the fact that we consider networks at their developmental stage, as explained above. Furthermore, the correlation of type T2 can also represent a homeostatic regulation of the neuronal firing to cope with different levels of synaptic inputs [39].
For clarity reasons, the paper will mainly deal with a specific realization of a network, made of N~100 neurons and embedding correlations of type T1 and T2. However, we have verified the validity of our findings in other five realizations of the network with correlations T1 and T2: four for N~100 (examined in Text S2) and one corresponding to N~200 (discussed in details in Text S3).

Single neuron stimulation/deletion impacts population bursting in developmentally correlated networks
In the developing hippocampus it has been shown how the stimulation of specific single neurons can drastically reduce the frequency of the PBs [5,16]. These neurons have been identified as hub cells for their high degree of functional, effective, and structural connectivity [40]. Stimulation consisted of phasic or tonic current injection capable of inducing sustained high firing regime of the stimulated neuron over a period of a few minutes. Based on this experimental observations, we tested the impact of prolonged single neuron stimulation (SNS) on the occurrence of PBs on our network model. SNS was obtained by adding abruptly a DC current term to the considered neuron. For illustrative purpose, we report in Fig. 1 A-B the stimulation protocol for a specific neuron capable of suppressing the occurrence of PBs for all the duration of the SNS (in this case limited to 12 s). During the current stimulation the neuron fired with a frequency of^36 Hz well above the average (6+5 Hz) and the maximal (24 Hz) firing rate of the neurons in the network under control conditions (see the bottom panel in Fig. 1 C).
The stimulation process was completely reversible and after the end of the SNS both the firing rate of the cell and the PBs frequency returned to the pre-stimulation control level. In order to evaluate the impact of single neuron perturbation on the collective dynamics, we considered the variation of PB frequency relative to control conditions (i.e. in absence of any stimulation). In Figs. 1 C,D the impact of a single neuron stimulation on the PBs frequency is reported for a classical Erdös-Rényi network (no correlations) and a network with embedded correlations T1 plus T2. Please notice that the neurons in Figs. 1C-D are ordered according to their average firing rate n under control conditions, which covered the interval ½0:03; 24:80 Hz. The comparison of panels C and D in Fig. 1 clearly shows that the correlated network is much more sensitive to single neuron stimulation. In fact, the SNS was able, for three neurons, to suppress the occurrence of PBs during their stimulation, while for approximately another half dozen of neurons the PBs were halved with respect to control conditions. The three most critical neurons c 1 , c 2 and c 3 were characterized, before stimulation, by firing frequencies larger than 3.3 Hz, and they lay among the top 33% fastest spiking neurons. On the contrary, in a network where no correlations were present, the SNS had only extremely marginal influence on the population activity (see Fig. 1D), although the firing rate distributions in the correlated and uncorrelated network were extremely similar (under control conditions) as shown in the bottom panels of Fig. 1 C and D.
In [24] it was shown that the elimination of a pool of neurons from an uncorrelated network encompassing short-term synaptic plasticity caused a strong reduction of the population bursts. In this work we repeated such numerical experiment with single cell resolution, i.e. we considered the influence of single neuron deletion, SND, on the neuronal response of the network (results reported in Fig. 2). For the network with correlations T1 plus T2, in four cases the SND led to the complete disappearance of PBs within the examined time interval, while for five other neurons their individual removal led to a decrease of the order of 30{40% in the frequency of occurrence of the PBs (Fig. 2 A). Three among these four critical neurons (namely, c 1 , c 2 and c 3 ) were also responsible for silencing the network during the SNS experiment performed with a stimulation current I stim~1 5:90 mV, as shown in Fig. 1C and Fig. 2 A. The same SNS experiment on the fourth critical neuron, labeled c 0 , reduced the PB frequency of about 40% with respect to control conditions (Fig. 1C). However, as we will show in the following, a SNS experiment performed on c 0 with a different injected current can also lead to a complete silence in the network. Notably, in analogy with the SNS, also this additional critical neuron c 0 impacting the PB occurrence under SND lays among the top 33% fastest spiking neurons. Differently from SNS, the removal of neurons with lower frequencies had almost no impact on the network dynamics. For uncorrelated networks the effect of SND was almost negligible, inducing a maximal variation in the bursting activity of the order of 10-15% with respect to the activity under control conditions (see Fig. S3C).
Other correlation setups. We have also analyzed the response to SNS and SND experiments in networks embedding all the possible combinations of the correlations setups T1-T3. In particular, we considered networks with positive correlation between structural in-degree (K I ) and out-degree As we are looking for strong impacts on the network dynamics, we identified the network as ''sensitive'' to SND (or SNS) whenever the PBs frequency was altered more than 90% with respect to the corresponding PB activity in control conditions. Therefore, we considered as significative only the modifications of the activity which were well beyond the statistical fluctuations in the population bursting, shown by the shaded gray area in panels C and D in Figs. S4, S5, S6, and S7.
In all the examined cases, despite the fact that the firing frequencies distributions were quite similar to the ones measured in the correlated network embedding setups T1 and T2, we did not observe significant modifications of the bursting activity by performing SNS and SND experiments on any neuron of the network (see panels C and D in Figs. S4,S5,S6,S7). The situation where SNS and SND had a larger effect on the network activity was for the correlations of type T2. In that specific case we observed that SND on 2 neurons (lying among the top 33% fastest spiking neurons, namely nw5 Hz) halved the bursting frequency of the network, and SNS on one of these 2 neurons had a similar effect. In all the other cases the PB activity was never perturbed more than 20-30%. Only the simultaneous presence of type T1 and type T2 correlations noticeably enhanced the sensitivity to SNS and SND, leading to the possibility to silence the network.

Structural and functional properties of the network
In order to gain some insight into the mechanisms underlying the reported response of the network, with correlations of type T1 plus T2, to SNS and SND experiments, we analyzed the structural and functional connectivity of the network in relation to the intrinsic excitability of the neurons. Functional connectivity (FC) analysis [41] was aimed at revealing time-lagged firing correlations between neuronal pairs, similarly to what described in [5] for the developing hippocampus. In particular, for every possible pairs of neurons (i, j) we cross-correlated their spike time series, with the exclusion of the spikes occurring within bursts, for which only the timestamp of the first spike was kept (see Methods). A functional connection directed from i to j was established whenever the activation of i reliably preceded the activation of j and viceversa  (see Methods). For each cell i, we calculated the functional outdegree (in-degree) D O i (D I i ), i.e. the number of cells which were reliably activated after (before) its firing.
As shown in the top panel of Fig. 2B, the four critical neurons, c 0 -c 3 , identified during the SNS and SND experiments, have very high functional out-degree, namely 83ƒD O i ƒ90. In particular, three of them (c 0 , c 1 and c 2 ) are ranked among the first four neurons with the highest functional out-degree. Therefore the critical neurons are reliably preceding the activation of most of the other neurons in the network. In addition, neurons c 0 and c 2 were supra-threshold (I b wV th , see Methods) and therefore firing tonically even if isolated from the network, while neuron c 1 was at threshold and c 3 below it (as shown in the central panel of Fig. 2B).
In contrast to their high functional out-degree, critical neurons were characterized by a low structural degree K T (total number of afferent and efferent connections), namely K T v16 with respect to an average value 23+13, as shown in the bottom panel of Fig. 2B. This result was a direct consequence of the anti-correlation imposed between total degree and excitability and this represented a crucial aspect for the emergence of the critical neurons.
In Fig. 2 C we report the results of SNS (SND) experiments as a function of D 0 , I b and K T of the stimulated (removed) neurons. The experiments on the neurons with high K T (the structural hubs, shown in Fig. S1 A and S1 B) influenced marginally the network bursting, apart for the single neuron stimulation of the two principal hubs which led to a moderate increase of the activity (see the bottom panels in Fig. 2 C). However SND on the same neurons had no significant effect. On the other hand, neurons with high functional out-degree D 0 (functional hubs) were quite relevant to sustain the collective dynamics. The removal of neurons with low D 0 (including the structural hubs) seemed almost not affecting the bursting properties of the network. Altogether, apart the stressed differences, the SNS and SND experiments appeared to give quite similar results.
The generality of these findings have been tested by performing SNS/SND experiments on other five different realizations of the network with embedded correlations of type T1 and T2, in all cases a small subset of neurons resulted to be critical in the sense discussed above (for more details see Text S2 and S3).

Network response during SNS: Dependence on the injected current
In order to further clarify the impact of varying the intrinsic excitability of single neurons on the network bursting activity, we have performed extensive analysis of the network response under SNS experiments for a wide range of stimulation currents, namely I stim [½14:7; 18:0 mV. In panels A and B of Fig. 3 it is summarized the impact on the bursting activity of the SNS for networks with type T1 and T2 correlations and without any correlations. SNS had really a minimal effect on the uncorrelated network: in this case the number of emitted PBs varied only up to a 20% with respect to control conditions. On the contrary, for the correlated network, SNS was able to silence the network over a wide range of currents when c 1 , c 2 and c 3 were stimulated. For the other neurons, SNS with high stimulation currents could also have the effect of promoting an increase of PBs up to 130-140% with respect to control conditions. In particular, an increase in the PB activity has been observed consistently for two structural hubs, whenever they are brought above the firing threshold, as shown in Fig. S1 C, and for other two neurons directly connected to these hubs. This behaviour is expected for an excitatory network without correlations, where the neurons with higher out-degree have usually the highest impact on the network [19]. However, the removal of each of these four neurons from the network did not influence the PB activity, furthermore they were passively recruited during bursting events.These results had an explanation in the fact that in control conditions the structural hubs were well below threshold, due to the anti-correlation between total degree and excitability, while by increasing the stimulation on these hubs we violated such constraint.
As shown in panels D, E, F of Fig. 3, for neurons c 1 , c 2 , and c 3 the bursting activity survived only in narrow stimulation windows located around, or just above the firing threshold value. A current variation DI b^0 :1{0:3 mV was, for these three neurons, sufficient to silence the network. The stimulation of the neuron c 0 , the one critical for SND but not for SNS (when I stim~1 5:90 mV), revealed the existence of very narrow anti-resonance windows (i.e. minima in the number of emitted PBs), as shown in Fig. 3 C. For very specific intrinsic excitability this neuron could effectively silence the network, but for generic excitation its influence on PBs activity was limited. The anti-resonances occurred (for I b w15:30 mV) at almost regular intervals: initially of width^0:2 mV and, at larger intrinsic excitability, of widtĥ 0:4 mV. This point will be further discussed and clarified in Sect. Time Orchestration.

The functional clique
The results reported above suggest that the four critical neurons c 0 ,c 1 ,c 2 ,c 3 , identified in the network with correlations T1 plus T2 should have a key role in the onset of the collective bursting. Therefore, we focused our analysis on the PB build up, i.e. we examined the events occurring in a time window of 25 ms preceding the peak of synchronous activation (for more details see Methods). In particular, we quantified how many times each single neuron participated in the build up of a PB. As we have verified, for the correlated network all the bursts were preceded by the firing of the four critical neurons, while in absence of correlations there was no neuron capable of reliably preceding every burst activation. The cross correlations between the timing of the first spike emitted by each critical neuron during the PB build up (see Methods) are shown in Fig. 4 A (blue histograms). This analysis revealed a precise temporal sequence in the neuronal activation, respectively c 0 ?c 1 ?c 2 ?c 3 , as shown also for a few representative bursts in Fig. 5 A,B (therefore the labeling assigned to these neurons). Interestingly, the same neurons did not show this precise temporal activation out of the PBs, as revealed by the red histograms in Fig. 4 A (see also Methods). Furthermore, the time sequence of the firing events of the critical neurons during the build up of the PB was quite well determined: c 0 anticipated the firing of c 1 of 3:94+0:5 ms, c 1 anticipated c 2 of 9:6+3:3 ms and c 2 anticipated c 3 of 3:3+1:0 ms. During the inter-burst periods we observed clear time lagged correlations only for the pair c 0 ?c 1 , presenting a direct synaptic connection, and in a weaker manner also for the pair c 2 ?c 3 . On the basis of the reported data, we can safely affirm that the critical neurons form a functional clique responsible for the onset of the PBs.

The role of plasticity
As clarified in [24], the bursting activity was due to the shortterm-synaptic depression. In particular PB emission could be related to the evolution of the fraction of synaptic resources in the recovered state, characterized by the variable X IN i (X OUT i ), averaged over the afferent (efferent) synapses of each neuron i (see Methods). The authors in [24] have shown that the fraction of synaptic resources, averaged over all the excitatory synapses, had a  In this Section, we want to address the question whether the variation of the effective strength of the synapses could be also responsible for the silencing of the network (with correlations of type T1 plus T2) during SNS experiments. So far we have clarified that the removal of any of the four neurons in the functional clique blocked the bursting activity, however it is not clear why a small stimulation of c 1 ,c 2 ,c 3 was capable also of blocking the PBs. As reported in Fig. 6 A, B the stimulation of neuron c 3 with a large current I stim~1 5:9 mV (as in the experiment reported in Fig. 1 A,B) reduced noticeably X OUT c3 , due to the high firing activity of the stimulated neuron. Analogous results have been found for all the other three critical neurons. For neurons c 1 , c 2 and c 3 this stimulation blocked the bursting activity of the network, thus inducing an almost complete recovery of the available resources of the afferent synapses, measured by X IN (as shown in Fig. 6 A for c 3 ). These results could suggest that SNS and SND experiments are indeed equivalent, since if the efferent synapses are extremely depressed, this could correspond somehow to remove the neuron from the network. However, SNS of neuron c 0 did not lead generically to the suppression of the bursting activity even if its efferent synapses were similarly depressed (as shown in Fig. 3 C). Furthermore, the synaptic depression could not explain the antiresonances in the bursting activity observed for SNS of c 0 and c 2 with different I stim . Since the time averaged synaptic strength, SX OUT T, exhibited only a smooth decrease as a function of I stim for all the four critical neurons (as well as for any generic neuron in the network), as shown in Fig. 6 C.

Time orchestration
As already mentioned, the roles of the four neurons in the functional clique of the network with type T1 and T2 correlations were quite well established, and just a precise firing time sequence could induce the population avalanche. To better understand the role of each critical neuron, it is necessary to point out that, under control conditions, the neurons c 0 and c 2 could fire even if isolated (since I b c0~1 5:07 mV and I b c2~1 5:30 mV were larger than V th ), c 1 was at threshold (I b c1~1 4:99 mV) and c 3 was the only neuron below threshold (I b c3~1 4:89 mV). This clearly explains, given the existing synaptic connection from c 0 to c 1 (see Fig. 4 B), the reason why c 0 entrained c 1 , both during the burst build up as well as during the inter-burst periods (see Fig. 4). Furthermore, from the results of the SNS experiments performed on c 1 and c 3 (Panels D and F in Fig. 3) one can observe that the network activity arrested whenever I stim wI b c0~1 5:07 mV for both these neurons (for comparison, note that the range of I stim reported in panel C is different from panel D,E and F). Therefore, whenever these two neurons fired faster than the clique leader c 0 , the burst activity, which should be triggered by a well determined sequence of events, would be terminated. Thus we can conclude that c 1 and c 3 could only be the followers of the dynamics dictated by the two supra-threshold neurons, and in particular by the leader c 0 .
As clearly shown in Fig. 7 A, exactly before a burst event (i.e. in the PB build up phase) neuron c 1 fired with a precise time lag after neuron c 0 (blue dashed line in the figure). However, the time lag DT c1,c0 between the firing of c 0 and c 1 needed some time after each bursting event to adjust to its pre-burst value. This could be interpreted also as an effective refractory period needed to the pair c 0 -c 1 to recover the proper entrainment favorable to the burst discharge. As shown in Fig. 7 A, the time evolution of the variable X c1,c0 , which measured the effective strength of the synapse connecting c 0 to c 1 , is directly connected to the duration of the time interval DT c1,c0 (or analogously to the effective refractory time of the entrainment c 0 -c 1 ). After a burst, X c1,c0 was noticeably depressed (reaching almost zero) and it slowly recovered its asymptotic value over a time scale dictated by T R c1,c0 . Indeed X c1,c0 was strongly oscillating due to the firings of c 0 , however the recovery of the pre-burst condition can be assessed by considering its extreme values (minima and maxima) both slowly increasing after the burst. The recover of the effective synaptic strength was associated to the adjustment of DT c1,c0 to the value taken during the build up of a PB. From Fig. 7 A, it is also evident that the fulfillment of this condition was not sufficient to induce another PB, since the PB could occur even a long time after the favorable pre-burst value was reached by DT c1,c0 .
Similar behaviors had been observed also for the synapse connecting c 2 to c 3 , although the firing of neuron c 2 alone was not sufficient to bring c 3 above threshold and therefore to initiate the PB. Indeed, the activation of c 3 , whose firing was fundamental to trigger the avalanche, was more complex. From a structural point of view, the neuron c 3 received inputs directly from c 1 and c 2 , while there were no synaptic connections between c 1 and c 2 (see Fig. 4 B). The entrained firing of the pair c 0 {c 1 followed by the firing of c 2 , within a precise time window, was required to induce c 3 to emit a spike and therefore a PB. This can be clearly appreciated from Fig. 7 B and C. In particular, in Fig. 7 B is reported a situation where c 2 fired at the right time after c 0 , but c 1 has fired too late to start an avalanche in the network (as previously explained the firing of c 1 was not yet entrained to that of c 0 ). Much more common is the situation reported in Fig. 7 C, where c 1 fired essentially always at the same time after c 0 , but instead the time delayDT c2,c0 in the firing of c 2 was extremely variable ranging from an almost coincidence with c 0 to a delay of 100 ms. The PB could occur only when c 2 fired in a precise time window following the activation of c 0 . Once noticed that the most part of the PB failures are due to c 2 and in a first attempt to understand the emergence of bursts in the network, we can focus only on the firing times of neuron c 0 and c 2 .
To get a deeper insight on this issue, let us consider the antiresonances (corresponding to minima in the PB activity) observed during the SNS experiments performed on c 0 (see Fig. 3 C). To interpret such minima we examined the firing periods T 0 and T 2 of the neuron c 0 and c 2 once isolated from the network. For the LIF model [30] these are simply given by T 0~tm ln½(I stim {V r )= (I stim {V th ) and T 2~tm ln½(I b c2 {V r )=(I b c2 {V th ), where I stim is the stimulation current acting on c 0 and I b c2 the intrinsic excitability of c 2 . As shown in Table 1 the PB minima were associated to rational ratios of these periods. This amounts to exact frequency locking of the firing of the two neurons [42], whenever this occurs the bursting activity is depressed or even suppressed. This because the build up of a burst relies on a precise temporal mismatch between the firing of neuron c 0 and c 2 , which, in the case of exact locking, can be achieved quite rarely or even never. Therefore, given the absence of any structural connection among these two neurons, the clique functionality relied on unsupervised coordinated firing of c 0 and c 2 .
In order to confirm this hypothesis, we developed a simple model to reproduce the results of the SNS experiment on c 0 . In particular, we assumed that c 0 and c 2 could be considered as two independently spiking neurons with their own firing periods determined by the stimulation current I stim for c 0 and by the intrinsic excitability for c 2 . Furthermore, we assumed that a PB is emitted with a certain probability (related to the synaptic depression induced by the stimulation) whenever c 0 and c 2 fired in the correct order and with a prescribed time delay (for more details see Methods). The results are reported in Fig. 8 and in the Table 1, the agreement is quite surprising due to the limited ingredients employed in the model. Furthermore, the fact that more than the 60% of the ''anti-resonances'' as well as the level of the PB activity were reproduced was a clear indication that the simple ingredients at the basis of the model represented the main mechanisms behind the PB build up process in this network. These mechanisms could be summarized as follows: the functional clique can be assumed to be composed of two structurally connected pairs c 0 {c 1 and c 2 {c 3 , where c 0 and c 2 fired tonically and independently one from the other. Any spike emission of c 0 induced a firing of c 1 , however to recruit c 3 and therefore to initiate the PB, also c 2 should deliver a spike, with the right time delay after c 1 . Therefore, if c 0 and c 2 fired with periods which were rational multiples one of the other it was unlikely to build up the PB. Since the synchronism among the two neurons did not allow c 1 to participate to the build up of the PB. The spike delivered by c 1 is fundamental to lead c 3 above threshold and to trigger the avalanche, but it should be emitted at the right moment, as clearly shown in Fig. 7 B and C.

Discussion
The aim of the present work was to identify neuronal network arrangements sensitive to single neuron perturbations, such as those induced by single neuron stimulation or deletion (or forced silencing). We choose as a benchmark model a random network of excitatory LIF neurons, connected via depressive synapses regulated by the TUM mechanism [24]. Such networks displayed spontaneous bursting activity also in absence of inhibition, as extensively described in the literature [24,26,32,43]. The choice of random topology was aimed at revealing the role of developmentally regulated neuronal excitability and connectivity gradients [5,[35][36][37][38], rather than specific topological configurations, in rendering network organization sensitive to single neuron perturbations.
The introduction of a positive correlation between in-and outdegree (T1) and a negative correlation between intrinsic neuronal excitability and total degree (T2), besides being justified from a developmental point of view, favors also the stabilization of the network activity. This because, as pointed out in [19], in an excitatory network the sensitivity to fluctuations is mainly due to cells with a high out-degree. Therefore, to avoid that their activation during spontaneous activity can cause network destabilization, a possible strategy is to impose an anti-correlation between their level of excitability and their degree, as done in the present work, or between in-and out-degree as shown in [19]. Furthermore, when correlations T1 and T2 were embedded in the network, single neuron deletion/stimulation of a few peculiar neurons strongly impacted the frequency of occurrence of population bursts. Most critical neurons, i.e. those capable of silencing the network when deleted or stimulated, shared common features: they constantly/reliably participated in the PB build up (i.e. they were functional hubs) and they had a quite low structural degree. These functional hubs formed a clique, where their precise ordered temporal activation was necessary for the burst generation. In the specific case here described, the clique was composed by two synaptically connected pairs, each composed of one neuron above and one below threshold. The burst could be triggered only when the first three neurons operated at precise time lags and the last neuron of the clique (which is just below threshold) was led to fire.
Each population burst caused the depletion of the synaptic resources, therefore another PB could occur only when the synaptic resources would be recovered, thus inducing an effective refractory time between two successive PBs. However, this is only In (B) a clear failure is indicated by red circles, in this case c 2 fired at the right time, but c 1 was too slow; in (C) neuron c 1 fires at the right moment several times (black dots are within the gray shaded area in the top panel), but the avalanche is not initiated until c 2 does not emit a spike within a precise time interval after the firing of c 0 . In all the figures, the data refer to control conditions. The blue horizontal dashed lines refer to the average value of DT c2,c0 or DT c1,c0 at the PB onset, while the shaded gray areas indicate the corresponding standard deviations. doi:10.1371/journal.pcbi.1003823.g007 a necessary, but not sufficient condition for PB triggering. The key element responsible for generating PBs was the unsupervised occurrence of a precise sequence of firing times of the two suprathreshold critical neurons, i.e. not mediated by any structural synaptic connection. On the other hand, the mode locking of the firing frequencies of these two neurons was instead responsible for anti-resonances associated to a drastic reduction of the PBs. For random networks, i.e. with no correlations, or embedding just one of the correlations of type T1, T2, T3 or the combination of type T1 and T3, we did not find any evidence of functional cliques and the mechanisms of network synchronization were much more robust and immune from single neuron perturbations (see Supplementary Information).
The activity of random uncorrelated networks has been previously examined in [24], in particular the authors have shown that the elimination of a pool of neurons (namely, 30 neurons, corresponding to the^8% of the excitatory population) led to the interruption of the bursting activity. The PBs were suppressed whenever the removed pool was composed by neurons with an intermediate firing rate (^1:3{2:5 Hz). These neurons were responsible for triggering the avalanches in the network, due to their effectively strong excitatory synapses and to their proximity to the firing threshold. From these findings, it is clear that in an uncorrelated network the PBs emerge due to a cooperative effect involving a large portion of neurons. On the contrary, the introduction of correlations of type T1 and T2 induces single neuron sensitivity as discussed in this paper.
Furthermore, our results show that the integration into a clique is the key element that can enable single neurons to impact the population dynamics, without any further topological requirements for the network architecture. The functional hubs forming and operating within the clique, are actively involved in generating network synchronizations and, as a consequence, capable to impact the network dynamics when stimulated. Therefore, without necessarily being structural or effective hubs, i.e. capable to cause a direct influence on the activity of many other nodes [40,44], they operate as operational hubs accordingly to the definition recently introduced in [45]. Similarly to the hub cells in the developing hippocampus whose stimulation was capable to drastically reduce the frequency of spontaneous network synchronization [5], the critical neurons presented in this paper have a very high functional connectivity and several of them are close to the firing threshold.
At variance with hippocampal hubs, critical neurons do not have a high structural degree. This is a consequence of the correlation imposed on the network where the excitability of the neurons is anti-correlated to the total structural degree of the cells. Indeed, in the correlated network studied in this work, the orchestration of the neuronal activity relies on the coordinated firing of a few critical ''young'' neurons (i.e. with a low structural The first column reports the stimulation currents I stim for which pronounced minima (anti-resonances) are observed in the stimulated PB activity during SNS experiment on c 0 (same data as in Fig. 3 C and red curve in Fig. 8  degree) mediated by their inter-connections. However, in real biological developing networks, it is possible that a further developmental connectivity regulation is fulfilled, with the chance of finding synaptic connections in a pair of young cells much lower compared to a pair composed of a young and a mature cells. This would be also in line to the rich gets richer rule which can generate scale-free networks [46]. In such case, the orchestration between unconnected young neurons would require the presence of a structural connector or hub, i.e. a more developed neuron, capable to receive and promptly activate in the presence of a few synchronized inputs. Therefore, our study supports the hypothesis that, in developmentally constrained networks, PBs are triggered by a precise time activation of a few around threshold oscillators. This is indeed the case of neurons c 1 and c 3 , which are fundamental for the ignition of the neuronal avalanche, but they need to be activated by a precise firing sequence involving c 0 and c 2 . This evidence is even more striking in the example discussed in Text S3 for N~200 neurons, where the functional clique is composed of a small group of neurons all just below threshold, apart the leader who activates the neurons in cascade leading to the burst.
We have verified that the main ingredients required to observe strong sensitivity to single neuron stimulation and deletion are, besides the presence of the correlations of type T1 and T2, a small number of neurons supra-threshold as well as a strongly diluted network. This can find an explanation in the fact that by increasing the degree of the neurons as well as the number of neurons suprathreshold the network dynamics becomes more cooperative. Furthermore, the synaptic time scales seem not to be crucial for the emergence of single neuron sensitivity (for more details see the subsection Dependence on the Model Parameters in Methods).
Although presented within a developmental framework, our results can also have elements of interest in the wider context of brain processing. In fact, in this work we show that the introduction of an excitability gradient in the network can lead to the emergence of functional cliques capable to shape the neuronal population activity. Indeed, different brain states could dynamically modulate the level of excitability or the gain function of the neurons within a circuit (as clearly discussed in [47]) and in this way instantaneously induce the emergence of functional cliques. Furthermore, functional chains of neural activation have been reported also in the different framework of feed-forward networks [48,49]. In this context, in Ref. [20] the authors found that structural hubs (i.e. highly connected neurons) have a peculiar role in promoting the signal transmission across sequences of nonhub sub-networks.
The previously discussed anti-resonance effect leading to the silence of the bursting activity resembles recent results reported in literature [28], where the authors have shown that abnormally synchronous neuronal populations can be desynchronized by administrating stimulations at resonant frequencies to an ensemble of spiking neurons. In that context, the desynchronization of the neuronal activity can be achieved by delivering a periodic stimulation at few sites, with a period which was an integer multiple of the fundamental period of the synchronized system. This is in striking contrast with what usually observed for a resonant forcing of a population of coupled oscillators [50,51]. Our results, revealing population desynchronization associated to anti-resonances at the level of single neuron frequencies, are even more intriguing. On one side these findings suggest the possibility of extremely non-invasive procedures to treat pathological neuronal synchronization, which is associated to several neurological disorders [52,53]. On the other side they reveal the potentiality of a brain circuit able to adapt to external stimuli on the basis of unsupervised mechanisms, which can switch the network activity from coherent to incoherent.
The numerical results here presented predict a primary role for supra-threshold and near-threshold cells capable to impact network synchronizations when organized into functional clique. Probing the existence of such cliques, whose emergence could be also dynamically regulated by varying the gradient of excitability in the circuits [47], is experimentally challenging, but surely feasible in in-vitro biological preparations. Cultured networks allow for an easier access and probing of the circuits [36,54], while representing a general model of unsupervised (or self-organized) spontaneous network synchronization in circuits under development, analogously to what observed in central and peripheral brain circuits [55,56]. In these circuits high functionally connected neurons (mostly activated at the build up of bursting) and highly active (i.e. supra-threshold) neurons could be identified by using both multi-electrode recordings and/or calcium imaging [54,57]. Furthermore, by manipulating the frequency of firing of such cells through multi-site optical or electrical stimulation [21,36] it is possible both to disrupt the sequential activation necessary for triggering network synchronizations (as displayed in Figs. 5 and 4) and to test the anti-resonance effects, as described in Fig. 8 and Table 1.
In this work, we considered a deterministic model of short term synaptic depression based on a trial-averaged representation. In recent papers, the stochastic processes involved in vesicle release and synaptic recovery time have been also taken into account to model short-term synaptic depression [58][59][60]. In particular, in Ref. [59] the authors compare deterministic and stochastic model for short-term plasticity. They found that for supra-threshold neurons the two setups give essentially the same behavior, while for sub-threshold neurons, whose spiking activity is fluctuation driven, the results of the deterministic and stochastic models essentially coincide for low frequencies (up to^10 Hz). In our study the functional hubs are found to be or supra-threshold or to have a relatively low firing rate (see bottom panel in Fig. 1 C). Therefore we expect that the implementation of stochastic shortterm plasticity would not affect qualitatively our main findings, but further investigations are required to fully clarify this issue.

Methods
To study the response of bursting neural networks to single neuron stimulation and removal, we employed the Tsodyks-Uziel-Markram (TUM) model [24]. Despite being sufficiently simple to allow for extensive numerical simulations and theoretical analysis, this model has been fruitfully utilized in neuroscience to interpret several phenomena [25,26,32]. We have considered such a model, restricted to excitatory synapses, to somehow mimic the dynamics of developing brain circuitries, which is characterized by coherent bursting activities, such as giant depolarizing potentials [5,29]. These coherent oscillations emerge, instead of abnormal synchronization, despite the fact that the GABA transmitter has essentially an excitatory effect on immature neurons [61]. The model uses leaky-integrate-and-fire (LIF) neurons with excitatory synapses displaying short-term synaptic depression [24] arranged in a directed random network. It should be stressed that we do not consider a network under topological development, which is typically characterized by a dynamical evolution (addition/ deletion) of the links among the neurons.

The model
In this paper we consider a network of N excitatory LIF neurons, interacting via synaptic currents regulated by short-term-plasticity according to the model introduced in [24]. The time evolution of the membrane potential V i of each neuron reads as where t m is the membrane time constant, I syn i is the synaptic current received by neuron i from all its presynaptic inputs and I b i represents its level of intrinsic excitability. The membrane input resistance is incorporated into the currents, which therefore are measured in voltage units (mV).
Whenever the membrane potential V i (t) reaches the threshold value V th , it is reset to V r , and a spike is sent towards the postsynaptic neurons. For the sake of simplicity the spike is assumed to be a d-like function of time. Accordingly, the spiketrain S j (t) produced by neuron j, is defined as, where t j (m) represent the m-th spike time emission of neuron j.
The transmission of the spike train S j to the efferent neurons is mediated by the synaptic evolution. In particular, by following [62] the state of the synaptic connection between the j-th presynaptic neuron and the i-th postsynaptic neuron is described by three adimensional variables, X ij , Y ij , and Z ij , which represent the fractions of synaptic transmitters in the recovered, active, and inactive state, respectively and which are linked by the constraint X ij zY ij zZ ij~1 . The evolution equations for these variables read as Only the active transmitters react to the incoming spikes S j : the adimensional parameters u ij tune their effectiveness. Moreover, fT I ij g represent the characteristic decay times of the postsynaptic current, while fT R ij g are the recovery times from synaptic depression.
Finally, the synaptic current is expressed as the sum of all the active transmitters (post-synaptic currents) delivered to neuron i where G i is the coupling strength, while E ij is the connectivity matrix whose entries are set equal to 1 (0) if the presynaptic neuron j is connected to (disconnected from) the postsynaptic neuron i. At variance with [24], we assume that the coupling strengths are the same for all the synapses afferent to a certain neuron i. We have verified that this simplification does not alter the main dynamical features of the TUM model under control conditions. In this paper we study the case of excitatory coupling between neurons, i.e. G i w0. Moreover, we consider a diluted network made of N~100 neurons where the i-th neuron has K I i (K O i ) afferent (efferent) synaptic connections distributed as in a directed Erdös-Rényi graph with average in-degree { K I~1 0, as a matter of fact also the average out-degree was { K 0~1 0. The sum appearing in (5) is normalized by the input degree K I i to ensure homeostatic synaptic inputs [39,63].
The propensity of neuron i to transmit (receive) a spike can be measured in terms of the average value of the fraction of the synaptic transmitters X OUT i (X IN i ) in the recovered state associated to its efferent (afferent) synapses, namely The intrinsic excitabilities of the single neurons fI b i g are randomly chosen from a flat distribution of width 0.45 mV centered around the value V th~1 5 mV, with the constraint that 10% of neurons are above threshold. This requirement was needed to obtain bursting behavior in the network. This choice ensures under control conditions that the distribution of the single neuron firing rates is in the range ½0:03; 25 Hz.
For the other parameters, we use the following set of values: t m~3 0 ms, V r~1 3:5 mV, V th~1 5 mV. The synaptic parameters fT I ij g, fT R ij g, fu ij g and fG i g are Gaussian distributed with averages T I~3 ms, T R~8 00 ms, u u~0:5 and G~45 mV, respectively, and with standard deviation equal to the half of the average. These parameter values are analogous to the ones employed in [24] and have a phenomenological origin.

Correlations
Furthermore, we have considered networks where correlations of type T1, T2 or T3 are embedded. Correlation T1 is obtained by generating randomly two pools of N{4 input and output degrees from an Erdös-Rényi distribution with average degree equal to 10. The degrees are ordered within each pool and then assigned to N{4 neurons in order to obtain a positive correlation between K O i and K I i . Finally, four hubs with total degree K T w50 are added to this N{4 neurons. The final total degree distribution is shown in Fig. S1B.
Correlation of type T2 (T3) imposes a negative (positive) correlation between excitability I b i and the total degree of the single neuron K T i~K I i zK O i . To generate this kind of correlation the intrinsic excitabilities are randomly generated, as explained above, and then assigned to the various neurons accordingly to their total connectivities K T i , thus to ensure an inverse (direct) correlation between I b i and K T i . Correlations of type T2 (T3) are visualized in Fig. S1 A and Fig. S6 A.

Numerical integration of the model
In order to have an accurate and fast integration scheme, we transformed the set of ordinary differential equations (1), (3) and (4) into an event-driven map [64] ruling the evolution of the network from a spike emission to the next one (see Text S1 for more details on the implementation of the event-driven map). It is worth to stress that the event-driven formulation is an exact rewriting of the dynamical evolution and that it does not involve any approximation.

Population bursts
In order to identify a population burst we have binned the spiking activity of the network in time windows of 10 ms. A population burst is identified whenever the spike count involves more than 25% of the neural population. In order to study the PB build up, a higher temporal resolution was needed and the spiking activity was binned in time windows of 1 ms. The peak of the activation was used as time origin (or center of the PB) and it was characterized by more than 5% of the neurons firing within a 1 ms bin. The time window of 25 ms preceding the peak of the PB was considered as the build up period for the burst. In particular, the threshold crossing times have been defined via a simple linear interpolation based on the spike counts measured in successive time bins.
These PB definitions gave consistent results for all the studied properties of the network. The employed burst detection procedure did not depend significantly on the precise choice of the threshold, since during the inter-burst periods (lasting hundreds of milliseconds) only 10-15% of neurons were typically firing, while more than 80% of the neuronal population contributed to the bursting event (lasting^30 ms).
The average interburst interval for the network with (without) correlations under control conditions was 586+183 ms (208+74 ms) for a network made of N~100 neurons, while the burst duration was 27+3 ms. Doubling the number of neurons in the correlated network did not affect particularly neither the average interburst, which became^500 ms, nor the burst duration (^24 ms). For a more detailed discussion of the dynamics of this larger network see the Text S3.

Dependence on the model parameters
In this subsection, we summarize the crucial ingredients needed to observe strong sensitivity to SNS/SND. In particular, we have checked, for different model parameters, when SNS/SND experiments were still able to noticeably modify the PB activity (i.e. more than 90% with respect to control conditions). Firstly, we considered the sparseness of the network, the reported results refer to a diluted network with a probability of 10% to have a link among two neurons. We have observed that the results of the SNS/SND experiments strongly depend on the level of dilution, however they can still be effective up to a connection probability of 50%. Another crucial aspect was the small number of neurons supra-threshold, in the studied case this number corresponded to the 10% of the neurons. By varying this percentage up to 20% we still observed that the network can be silenced by single neuron stimulation/removal. Furthermore, the dependence on the system size seemed not be crucial, since as described in Text S3 by doubling the system size a functional clique can still be identified.
We also tested the influence of the synaptic time constants on the population sensitivity to SND/SNS. As a matter of fact, the time scale ruling the depletion of the neurotransmitter T I affects the duration of the PB, while the recovery time from the synaptic depression T R influences the intervals between consecutive PBs [24]. By varying T I within the interval ½1:5; 4:5 ms, while keeping fixed the ratio with T R , we do not observe substantial modifications on the network dynamics. Furthermore we found strong response to SNS and SND, leading to the population silence for several stimulated neurons, both for faster and slower synaptic time scale than the one actually employed in the studied example. However, it should be noticed that by increasing T I to extremely large values (namely, T I w20{30 ms) this destroys the bursting behaviour in the network, which is then substituted by an asynchronous activity [43].

Functional connectivity
In order to highlight statistically significant time-lagged activations of neurons, for every possible neuronal pair we measured the cross-correlation between their spike time series. On the basis of this cross-correlation we eventually assign a directed functional connection among the two considered neurons, similarly to what reported in [5,54] for calcium imaging studies.
Let us explain how we proceeded in more details. For every neuron, the action potentials timestamps were first converted into a binary time series with one millisecond time resolution, where ones (zeros) marked the occurrence (absence) of the action potentials. Given the binary time series of two neurons a and b, the cross correlation was then calculated as follows: where fa t g,fb t g represented the considered time series and T was their total duration. Whenever C ab (t) presented a maximum at some finite time value t max a functional connection was assigned between the two neurons: for t max v0 (t max w0) directed from a to b (from b to a). A directed functional connection cannot be defined for an uniform cross-correlation corresponding to uncorrelated neurons or for synchronous firing of the two neurons associated to a Gaussian C ab (t) centered at zero. To exclude the possibility that the cross correlation could be described by a Gaussian with zero mean or by a uniform distribution we employed both the Student's t-test and the Kolmogorov-Smirnov test with a level of confidence of 5%. The functional out-degree D O i (in-degree D I i ) of a neuron i corresponded to the number of neurons which were reliably activated after (before) its firing.
Time series surrogates. In order to treat as an unique event multiple spike emissions occurring within a PB, different time series surrogates were defined for different kind of analysis according to the following procedures: 1. for the definition of the functional in-degree D I i and out-degree D O i , all the spiking events associated to an inter-spike interval longer than 35 ms were considered. Since we observed that this was the minimal duration of an inter-spike outside a PB and it was larger than the average duration of the PBs. This implies that for each neuron only the timestamp of the first spike within a PB was kept; 2. for the description of the PBs build up only the timestamps of the first action potential emitted within a window of 25 ms preceding the PB peak was taken into account; 3. for the analysis of the network activity during inter-burst periods, all action potentials emitted out of the PBs were considered.

A simple model for SNS
The model here reported has been developed to reproduce the network response during the SNS experiments on c 0 for a range of stimulation current I stim , which is displayed in Fig. 3 C. To mimic this activity, we only considered the dynamics of the two neurons of the clique c 0 and c 2 which were supra-threshold. In particular, we assumed that these two neurons fired tonically and independently as they would be two isolated LIF neurons (oscillators) [30]. Therefore, as a first step we generated two regular spike trains, one for c 0 and one for c 2 with constant inter-spike time intervals T 0~tm ln½(I stim {V r )=(I stim {V th ) and T 2~tm ln½(I b c2 {V r )= (I b c2 {V th ), respectively. Successively, by examining the two spike trains, we assumed that a PB in the network could occur whenever the neuron c 2 fired after c 0 with a certain time delay t D .
Furthermore, we also assumed that the PB emission was a probabilistic event with a finite probability P~P(I stim ). P is simply given by the average efferent synaptic strength SX OUT c0 T measured during the SNS experiment on c 0 suitably rescaled in order to get probability one for I stim~V th . This probability has been introduced to mimic the decrease of the effective synaptic strength induced by the increasing stimulation, as shown in Fig. 6 C.
In summary, for each stimulation current I stim we considered the sequence of the firing times of c 0 and c 2 and we registered the occurrence of a PB whenever the two following conditions were both fulfilled N c 2 fired after c 0 within a time window ½t D {s D ; t D zs D , where t D~1 4 ms is the average time delay DT c2,c0 measured immediately before a population burst and s D~4 ms is its standard deviation, both measured in control conditions; N a random number r extracted by a flat random distribution with support ½0 : 1 was smaller than P(I stim ).
Furthermore, each time a PB was registered, for the subsequent 27+3 ms (corresponding to the duration of a PB under control conditions) no further PB could be counted. The results of this simple model are reported in Fig. 8 and Table 1 together with the numerical results obtained by the simulation of the network activity. with various stimulation currents in the interval I stim [½14:5 : 16:0 mV. The blue vertical dashed lines and the magenta horizontal solid lines mark, resp., the value of the intrinsic excitability and the bursting activity of the network during control conditions. The number of PBs are measured over a time interval Dt~84 s. (E) Color coded rates of emission of PBs during SNS experiment performed for a range of injected DC currents I stim (y-axis). The PB rates during SNS are normalized to the PB rate in resting conditions. Neurons are ordered according to their functional outdegree rank (x-axis). The number of PBs are measured over a time interval