# Modeling driver cells in developing neuronal networks

Stefano Luccioli, David Angulo-Garcia, Rosa Cossart, Arnaud Malvache, Laura Mòdol, Vitor Hugo Sousa, Paolo Bonifazi, Alessandro Torcini

<figure><img src="image_1.png" /></figure>

## To cite this version:

Stefano Luccioli, David Angulo-Garcia, Rosa Cossart, Arnaud Malvache, Laura Mòdol, et al., Mod-
eling driver cells in developing neuronal networks. PLoS Computational Biology, 2018, 14 (11),
pp.e1006551. 10.1371/journal.pcbi.1006551 . hal-01963568

HAL Id: hal-01963568

https://amu.hal.science/hal-01963568

Submitted on 21 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# Modeling driver cells in developing neuronal networks

Stefano Luccioli¹,², David Angulo-Garcia³, Rosa Cossart⁴, Arnaud Malvache⁴, Laura Módol⁴, Vitor Hugo Sousa⁴, Paolo Bonifazi⁵,⁶‡, Alessandro Torcini¹,⁴,⁷†*

1. CNR - Consiglio Nazionale delle Ricerche - Istituto dei Sistemi Complessi, Sesto Fiorentino, Italy, 2. INFN Sez. Firenze, Sesto Fiorentino, Italy, 3. Grupo de Modelado Computacional - Dinámica y Complejidad de Sistemas, Instituto de Matemáticas Aplicadas, Universidad de Cartagena, Cartagena de Indias, Colombia, 4. Aix Marseille Univ, INSERM, INMED, Marseille, France, 5. Biocruces Health Research Institute, Bilbao, Vizcaya, Spain, 6. Ikerbasque: The Basque Foundation for Science, Bilbao, Spain, 7. Laboratoire de Physique Théorique et Modélisation, Université de Cergy-Pontoise, CNRS, UMR 8089, Cergy-Pontoise, France

‡ These authors are joint senior authors on this work.

* alessandro.torcini@u-cergy.fr

## Abstract

Spontaneous emergence of synchronized population activity is a characteristic feature of developing brain circuits. Recent experiments in the developing neo-cortex showed the existence of *driver cells* able to impact the synchronization dynamics when single-handedly stimulated. We have developed a spiking network model capable to reproduce the experimental results, thus identifying two classes of driver cells: functional hubs and low functionally connected (LC) neurons. The functional hubs arranged in a *clique* orchestrated the synchronization build-up, while the LC drivers were lately or not at all recruited in the synchronization process. Notwithstanding, they were able to alter the network state when stimulated by modifying the temporal activation of the functional clique or even its composition. LC drivers can lead either to higher population synchrony or even to the arrest of population dynamics, upon stimulation. Noticeably, some LC driver can display both effects depending on the received stimulus. We show that in the model the presence of inhibitory neurons together with the assumption that younger cells are more excitable and less connected is crucial for the emergence of LC drivers. These results provide a further understanding of the structural-functional mechanisms underlying synchronized firings in developing circuits possibly related to the coordinated activity of cell assemblies in the adult brain.

## Author summary

There is timely interest on the impact of peculiar neurons (*driver cells*) and of small neuronal sub-networks (*cliques*) on operational brain dynamics. We first provide experimental data concerning the effect of stimulated driver cells on the bursting activity observable in the developing entorhinal cortex. Secondly, we develop a network model able to fully reproduce the experimental observations. Analogously to the experiments two types of driver cells can be identified: functional hubs and low functionally connected (LC) cells. We explain the role of hub neurons, arranged in a clique, for the orchestration of the

<figure><img src="image_4.png" /></figure>

**OPEN ACCESS**

**Citation:** Luccioli S, Angulo-Garcia D, Cossart R, Malvache A, Módol L, Sousa VH, et al. (2018) Modeling driver cells in developing neuronal networks. PLoS Comput Biol 14(11): e1006551.
https://doi.org/10.1371/journal.pcbi.1006551

**Editor:** Michele Giugliano, BELGIUM

**Received:** February 6, 2018

**Accepted:** October 7, 2018

**Published:** November 2, 2018

**Copyright:** © 2018 Luccioli et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** The authors acknowledge financial support from the European Commission under the program "Marie Curie Network for Initial Training", project N. 289146, "Neural Engineering Transformative Technologies (NETT)" (SL, DAG, AT), by the A* MIDEX grant (No. ANR-11-IDEX-0001-02) and by the I-Site Paris Seine Excellence Initiative (No. ANR-16-IDEX-0008), both funded by the French Government programme "Investissements d'Avenir" (AT and DAG), from
Ikerbasque (The Basque Foundation for Science)
(PB) and from the Ministerio Economía, Industria y
Competitividad of Spain (grant SAF2015-69484-R
(MINECO/FEDER)) (PB). The funders had no role in
study design, data collection and analysis, decision
to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared
that no competing interests exist.

bursting activity in control conditions. Furthermore, we report a new mechanism which
can explain why and how LC drivers emerge in the structural-functional organization of
the entorhinal cortex.

**Introduction**

Coordinated neuronal activity is critical for a proper development and later supports sensory
processing, learning and cognition in the mature brain. Coordinated activity represents also
an important biomarker of pathological brain states such as epilepsy [1]. It is therefore essen-
tial to understand the circuit mechanisms by which neuronal activity becomes coordinated at
a population level. A series of experimental results indicates that non-random features are
clearly expressed in cortical networks [2-4], in particular neuronal sub-networks, termed *cli-
ques* [5], have been shown to play a fundamental role for the network activity and coding both
in experiments [6-10] as well as in models [11-14].

The identification of these small highly active assemblies in the hippocampus [6] and in the
cortex [7-9] poses the question if these small neuronal groups or even single neurons can
indeed control the neural activity at a mesoscopic level. Interestingly, it has been shown that
the stimulation of single neurons can affect population activity in vitro as well as in vivo
[15-25]. The direct impact of single neurons on network and behavioral outputs demonstrates
the importance of the specific structural and functional organization of the underlying cir-
cuitry. Neurons having such a network impact were recently termed *operational hubs* [26] or
*driver cells* [24]. It is thus critical to understand how specific network structures can empower
single driver cells to govern network dynamics. This issue has been addressed experimentally
in some cases. More specifically, in the developing CA3 region of the hippocampus, single
GABAergic hub neurons with an early birthdate were shown to coordinate neuronal activity.
These cells have a high functional connectivity degree, reflecting mainly the fact that they are
activated at the onset of Giant Depolarizing Potentials (GDPs), as well a high effective connec-
tivity degree [17]. This therefore represents a simple case where the circuit mechanism, pro-
moting a cell to the role of hub, is due to their exceptional number of anatomical links. But the
picture can be quite different in other brain regions, as recently demonstrated in the develop-
ing Entorhinal Cortex (EC) [24], where the driver cell population comprises both cells with a
high functional out-degree, as well as low functionally connected (LC) cells.

In order to understand the circuit mechanisms by which even a LC cell can influence popu-
lation bursts we have upgraded and modified a network model based on excitatory leaky inte-
grate-and-fire (LIF) neurons [12], previously developed to reproduce the functional properties
of hub neurons in the developing hippocampal CA3 area [17]. In such a model the *population
bursts* (PBs), corresponding to GDPs in neonatal hippocampus [27], were controlled by the
sequential and coordinated activation of few functional hubs. Notably, the perturbation of one
of these neurons strongly impacted the collective dynamics and brought even to the complete
arrest of the bursting activity, similarly to what experimentally found for the developing hippo-
campus in [17]. The model described in this paper contains two main differences with respect
to the hippocampal model [12].

Firstly, it comprises both inhibitory and excitatory neurons, to account for the fact that,
even though GABA acts as an excitatory neurotransmitter at early postnatal stages, some more
developed neurons have already made the switch to an inhibitory transmission at the end of
the first postnatal week in mice (P8), where most experimental data was obtained [28-30]. In
that respect, the maturation profile of the entorhinal-hippocampal circuit was recently
analyzed and it was shown that stellate cells in layer 2 of the medial EC were the first to mature, followed, in chronological order, by pyramidal cells in MEC layer 2 and CA3 [31]. This indicates a higher probability to find inhibitory GABAergic synapses in EC with respect to CA3 and it justifies our choice to include some inhibitory neuron in the model. Secondly, the developmental profile of the network is regulated only by the correlation between neuronal excitability and connectivity, while in [12] a further correlation was present. The anti-correlation between intrinsic excitability and the synaptic connectivity reproduces to some extent the homeostatic regulation of the intrinsic excitability described during neural development [32]. This model nicely mimicked the experimental observations in the EC similarly displaying the presence of driver cells with both low and high functional connectivity.

The paper is organized as follows. Firstly, we will report experimental evidences of how the stimulation of single LC drivers can impact network synchronization in the developing EC. Secondly, we will show that our simple model can nicely reproduce such results. This to validate its applicability in order to understand the synchronization dynamics of the EC at the early stage of the development. Then, we will present a full characterization of the numerical model leading to a complete understanding of the mechanism underlying the PB generation and the impact of driver LC cells on population dynamics. We will conclude with a discussion of the obtained results and of their possible relevance in the context of neuroscience.

# Results

## Experimental evidences of driver LC cells

The main experimental observation at the rationale of this work is the existence of *driver cells* (or Operational Hubs [26]) in the mice EC during developmental stage [24]. Driver cells have been identified using calcium imaging experiments and they were characterized by the capability to impact network synchronization (namely, GDPs' occurrence) when externally activated/stimulated through intra-cellular current injection. Two classes of driver cells were identified: (i) those with high directed functional connectivity out-degree, early activated and playing a critical role in the network synchronizations (*driver hub cells*) and (ii) those recruited only in the later stages of the synchronization build up, which therefore are low functionally connected (*driver LC cells*).

The experimental setup used to identify, target and probe the single-handedly impact of neurons on spontaneous EC synchronization is schematized in Fig 1(a,E) and S1 Fig. In brief, the functional connectivity of the cells has been measured during the spontaneous activity session, which preceded the single neurons' stimulation session, both lasting two minutes. A directed functional connection from neuron A to B was established whenever the firing activity of A significantly preceded the one of neuron B (more details can be found in Methods). The functional out-degree $D_i^O$ of a neuron i corresponded to the percentage of imaged neurons which were reliably activated after its firing. Neurons in the 90% percentile of the connectivity distribution were classified as hub neurons early activated in the network synchronization.

The protocol used for probing the impact of single neurons on the network dynamics was organized in three phases, each of two minutes duration: (1) a pre-stimulation resting period; (2) a stimulation period, during which a series of supra-threshold current pulses at a specific frequency $v_S$ have been injected into the cell; (3) a final recovery period, where the cell is no more stimulated. The frequencies $v_S$ employed to stimulate the single neurons have been selected to be of the order of magnitude of the average GDP frequency with the aim of revealing cell-network interaction.

A cell was identified as a driver whenever the distributions of the Inter-GDP-Intervals (IGIs) were significantly different in the stimulation period with respect to the pre-stimulation
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Impact of single-handedly stimulation of LC drivers on the collective dynamics of the entorhinal cortex and of the model. The left column (x.E) refers to experimental measurements taken in slices of EC, while the right column (x.S) to the numerical results. **Experiment.** The first panel (a.E) presents the image of a slice loaded with the calcium sensor where the active neurons are encircled and the pipette indicates the neuron targeted for intracellular stimulation. Scale bar: 100 µm. Data for a neuron with functional out-degree $D^O \approx 7\%$ are reported in (b-d.E), data for a neuron (from a different slice) with $D^O \approx 8\%$ are shown in (e-g.E). Panels (b.E) and (e.E) are boxplot of the IGIs for each experimental phase. Panels (c.E) and (f.E) represent the fraction of recruited cells participating in the GDP. During the stimulation period (red curves) a single cell is stimulated with a frequency $v_S = 0.33$ Hz ($v_S = 0.14$ Hz) for the first (second) neuron according to the protocol discussed in the text. Panels (d.E) and (g.E) report the phase $\Phi$ of the GDP versus the GDP index (denoting the successive GDP cycles). $\Phi$ is the difference between the number of observed versus expected GDPs relative to the pre-stimulation frequency (see Methods). The average IGI was 4.2 s (7.8 s) in the pre-stimulation phase, becoming 3.2 s (14 s) during the stimulation period for the first (second) neuron. **Model.** The first row (a.S) displays a cartoon of the performed single neuron stimulation experiment in the network model, where inhibitory (inh.) and excitatory (exc.) neurons and synapses are marked in green and black respectively. Panels (b-d.S) refer to the excitatory driver LC cell $el_1$ which was silent in control condition and once stimulated with a current $I^{stim} = 15.135$ mV was able to enhance population dynamics. Panels (e-g.S) refer to the excitatory driver LC cell $el_2$ which was active in control condition with $D^O = 3\%$ and once stimulated with a current $I^{stim} = 15.435$ mV depressed the PB activity.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1006551.g001

and recovery periods (see Methods for more details). A further indicator that we use to visualize the effect of the stimulation on the GDP deliver is the shift of the IGI phase $\Phi$ measured with respect to the pre-stimulation phase (for the definition see Methods Eq (1) and in [33]). At a population level the stimulation may have an inhibitory (excitatory) effect corresponding
to a slow-down (acceleration) of the GDP frequency associated with an increase (decrease) of the measured IGI corresponding to a positive (negative) phase shift. Two examples of driver LC cells, with $D^0 \approx 7-8\%$, are reported in Fig 1 in the panels (b-d.E) and (e-g.E). In the first case, upon stimulation the network dynamics accelerated, as testified by the decrease of the average IGI (Fig 1(b.E)) and by the negative instantaneous phase shift of GDPs (Fig 1d.E). In the second case, the stimulation led to a pronounced slow down of the average network activity (as shown in Fig 1(e.E)) together with an increase of the instantaneous phase with respect to control conditions (Fig 1(g.E)). In both cases the removal of the stimulation led to a recovery of the dynamics similar to the control ones.

A further extreme case of a silent cell, i.e. not spontaneously active and therefore with a zero (out-degree) functional connectivity, is shown in Fig 2. This cell, when stimulated with different stimulation frequencies $v_S$, revealed opposite effects on the network behaviour. At lower stimulation frequency ($v_S$ = 0.33 Hz) the cell activity induced an acceleration of the population dynamics (see Fig 2(a)-2(c.E)), while at higher stimulation frequency ($v_S$ = 1 Hz) of the same neuron we observed a slowing down of the network dynamics (see Fig 2(d)-2(f.E)). By considering a much larger pool of driver cells we have verified that the value of $v_S$ does not induce a systematic trend towards a slow down or an acceleration of the GDPs (for more details see S2 Text and S3 Fig).

We also tested the possibility that single neuron stimulation could modify other features of network synchronization besides the GDP frequency, in particular we focused on the width of the GDPs, i.e. on their time duration as defined in S1 Text. However, as shown in S2 Fig we could not find any statistical differences between stimulation epoch and pre/post stimulation

<figure><img src="image_6.png" /><figcaption>**Fig 2.** The same driver LC cell upon different stimulations can enhance or depress the population activity. The left column (x.E) refers to experimental measurements taken in slices of entorhinal cortex, while the right column (x.S) to the numerical results. Panels are similar to Fig 1 but reporting the case in which the network acceleration (a-c.E and a-c.S) and slowing down (d-f.E and d-f.S) are observable by stimulating the same neuron at different frequencies. The experiment and the simulation refer to neurons with zero functional out-degree. **Experiment.** Panels (a-c.E) ((d-f.E)) are obtained with $v_S$ = 0.33 Hz ($v_S$ = 1.0 Hz). Namely, the average IGI varied from 6.08 s in control conditions to 5.14 s during stimulation with $v_S$ = 0.33 Hz, and from 5.36 s to 6.68 s with $v_S$ = 1.0 Hz. **Model.** All the panels refer to the driver LC cell $el_3$ connected in output to the $el_1$ neuron discussed in panels (b-d.E) in Fig 1. The results refer to the stimulation of the same neuron with two different currents, namely panels (a-c.S) refer to $I^{stim}$ = 15.42 mV and (d-f.S) to current $I^{stim}$ = 15.9 mV.</figcaption></figure>
periods (see S1 Text for details). This doesn't mean that other more subtle features of the
GDPs could not be impacted by manipulating single neurons, but probably these modifica-
tions could be hardly discernible in view of the limited time resolution associated to calcium
imaging.

Numerical evidences of driver LC cells

In order to mimic the impact of single neurons on the collective dynamics of a neural circuit, we considered a directed random network made of *N* LIF neurons [34, 35] composed of excitatory and inhibitory cells and with synapses regulated by short-term synaptic *depression* and *facilitation*, analogously to the model introduced by Tsodyks-Uziel-Markram (TUM) [36]. In particular, synaptic depression was present in all the connections, while facilitation only in the connections targeting inhibitory neurons. This in agreement with recent experimental investigations of Layer II of the medial entorinal cortex reporting evidences of short-term synaptic depression among excitatory neurons (namely, stellate and pyramidal cells) as well as among fast-spiking interneurons and pyramidal cells and of short-term facilitation in the connections from stellate cells towards low-threshold-spiking interneurons [37] (see *Methods* for more details).

As shown in [36, 38, 39], these networks exhibit a dynamical behavior characterized by an alternance of short periods of quasi-synchronous firing (PBs) and long time intervals of asyn- chronous firing, thus resembling cortical GDPs’ occurrence in early stage networks. Similarly to the modeling reported in [12], we considered neuronal intrinsic excitabilities negatively cor- related with the total connectivity (in-degree plus out-degree) (for more details see *Definition of the Model* in *Methods* and S4 *Fig*). The introduction of these correlations was performed to mimic developing networks, where both mature and young neurons are present at the same time associated to a variability of the structural connectivities and of the intrinsic excitabilities. Experimental evidences point out that younger cells have a more pronounced excitability, most likely due to the fact that their GABAergic inputs are still excitatory [40–42], while mature cells exhibit a higher number of synaptic inputs and they do receive inhibitory or shunting GABAergic inputs [17, 43]. The presence of inhibition and facilitation are the major differences from the model developed in [12] to simulate the dynamics of hippocampal circuits in the early stage of development, justified by the possible presence of mature GABAergic cells in the network.

Using this network model, we studied the effect of single neuron current injection $I^{\text{stim}}$ on network dynamics, thus altering the average firing frequency of the neuron during the stimulation time, similarly to what done in the experiments. In the numerical investigations, at variance with the experiments, the stimulation delivered to the neurons is an unique supra-threshold step of duration of 48 seconds. In Fig 1 two representative driver LC cells are reported for comparison with the experiments. The first cell (panels (b-d.S) of Fig 1) was a silent neuron in control conditions (therefore with $D^O = 0$), that once stimulated could enhance of ≈ 30% the PB emission, thus leading to a decrease of the instantaneous phase $\Phi$ with respect to control condition. Panels (e-g.S) refer to a second neuron characterized by a low functional output connectivity, namely $D^O = 3\%,$ whose stimulation led to a depression in the PB frequency (as shown in panels (e.S) and (f.S)) joined to an increase of the instantaneous phase of the network events with respect to control conditions (as shown in panel (g.S)). These results compare quite well with the experimental findings reported in the same figure.

Furthermore, analogously to what found in the experiment, Fig 2(a)–2(f.S) shows a silent neuron in control condition that once stimulated could lead to both enhancement or depression of the population activity depending on the level of injected current during stimulation.
A full characterization of the network model concerning the impact on the network dynam-
ics of each single neuron stimulation in relation to neuronal type, current injected and func-
tional connectivity is detailed below.

Impact of single neuron stimulation and deletion on network dynamics

In order to explore the full dynamical range associated to the impact of single neuron stimula-
tion on the network dynamics, we examined the response of the model network to two types
of single neuron perturbations, i.e. *single neuron deletion* (SND) and *single neuron stimulation*
(SNS) by employing the protocols introduced in [12]. In particular, the SND experiment con-
sisted in recording the activity of the network in a fixed time interval Δ*t* = 84 s when the con-
sidered neuron was removed from the network itself. While, the stimulation of the single
neuron (SNS) was performed with a step of DC current of amplitude *I*<sup>stim</sup> for a time window
Δ*t* = 84 s. The recording of the activity in control condition was lasting 84 s as well, in order to
compare directly the number of observed PBs during control and perturbation period. In par-
ticular, we tested the response of the network to a broad range of stimulation amplitudes vary-
ing from 14.5 mV (slightly below the firing threshold for an isolated neuron V<sub>th</sub> = 15 mV, see
*Methods*) to 18.0 mV with a step of 0.015 mV, inducing in the stimulated neuron a maximal
firing frequency of ≈70 Hz. Typically the stimulated neuron fired with a frequency much
higher than the frequency of neurons under control conditions (i.e. in absence of any pertur-
bation). As an example, for a stimulation current *I*<sup>stim</sup> = 15.90 mV the targeted neuron fired at
a frequency *v* ≈ 32 – 36 Hz well above the average (≈3 Hz) and the maximal (22 Hz) frequency
of all neurons in control conditions. The SND represented an extreme version of the SNS,
where the neuronal removal corresponded to the injection of an hyperpolarizing current
inhibiting the neurons from firing spontaneously or in response to any synaptic input.

In both SNS and SND experiments the impact of single neuron perturbation on the collective dynamics was measured by the variation of the PB frequency relative to control conditions. In general, we have classified a neuron as a driver cell whenever upon stimulation it is able to modify the PB frequency of at least 50% with respect to control conditions. In the specific, in analogy to what done in [24], for SNS experiments we considered both enhancement and decrease in the PB activity. On the other hand, SND allowed us to directly identify the driver neurons which are fundamental for the PB build-up. Therefore in this case we limited to consider those cells whose SND led to a population decrease of at least 50%. The choice of this threshold is based on the analysis of the distribution of the number of PBs delivered during SNS and SND experiments compared to the PB variability observed in control conditions. This was preliminarily measured by considering the average and the standard deviation of the number of population bursts obtained by considering for each network realization 100 different initial conditions over a time window of 84 s. As we were interested in strong impact on the network dynamics, we decided to consider as significant the variations of the population activity which were well beyond the statistical fluctuations in the population bursting measured by three standard deviations. The choice of 50% fulfilled this condition in all the analyzed networks.

Fig 3(a) and 3(b) report a comparison of the impact of SND and SNS (with representative injected current of 15.90 mV) on the PB activity. The removal of any of the four neurons labeled as ih₁, eh₁, eh₂, eh₃ was able to arrest completely the bursting dynamics within the considered time window, while in other two cases (for neurons ih₂ and eh₄) the activity was reduced of 60% with respect to the one in control conditions. For clarity, the used labels i/e stand for inhibitory/excitatory and h for hub, as we will show later this is related to the functional role played by these cells. For all the other neurons, the SND manipulation induced a
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Model—Network response to single neuron stimulation. (a-b) Response of the network to SND and SNS, respectively. These panels report the number of PBs, recorded during SND (SNS) experiments versus the labels of the removed (stimulated) neuron, ordered accordingly to their average firing rates $v$ under control conditions (shown in panel (c)). Inhibitory neurons are marked with asterisks. In this representative SNS experiment each neuron was stimulated with a DC step $I^{\text{stim}}$ = 15.90 mV for a time interval $\Delta t$ = 84 s. The central horizontal dashed line shows the average number of PBs emitted in control conditions within an interval $\Delta t$ = 84 s, while the lower and upper horizontal dashed lines mark the 50% variation. The vertical dashed line separates firing neurons (on the right side) from silent neurons (on the left side) in control conditions. (d) Raster plot of the network activity. (e) Close ups of population bursts representative of the two principal routes (i.e the firing sequence of hub cells): an instance of the main route is reported in the left panel, while the second most frequent route is displayed in the right panel. Note that neurons are labeled accordingly to their index and are not ordered as in panels (a-c). (f) Scatter plot showing the functional out degree $D^o$ of the neurons versus their total structural degree $K^t$. A double square marks two neurons with overlapping properties. (g) Sketch of the structural connections among driver hub neurons and driver LC1 cells, the large rectangle highlights the clique of the hub cells. In all the panels, circles (squares) denote driver hub (LC) cells, while the green (red) symbols refer to inhibitory (excitatory) driver cells.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1006551.g003

non relevant modification in the number of emitted PBs, within the variability of the bursting activity in control conditions (Fig 3(a)).

The SNS confirmed that the neurons *ih*₁, *eh*₁, *eh*₂, *eh*₃ were capable to arrest the collective dynamics. Neurons *eh*₄ and *ih*₂ poorly impacted PB dynamics for the reported injected current, although for different values of *I*<sup>stim</sup> they were able to strongly influence the network dynamics (as shown in the subsection *Tuning of PBs frequency upon hubs' and driver LC cells' stimulations*). At variance from what found in a purely excitatory network [12], the SNS revealed also the presence of other 18 driver cells not identified by the SND capable to impact the occurrence of PBs in the network (Fig 3(b)).
For an equivalent random network, without any imposed correlation, SNS or SND affected the dynamics in a negligible way producing a maximal variation of the bursting activity of 25-30% with respect to the control conditions (see S5(a) and S5(b) Fig).

To summarize, the presence of correlations among the neuronal intrinsic excitabilities and the corresponding structural connectivities was crucial to render the network sensible to single neuron manipulation. Differently from purely excitatory networks where SNS and SND experiments gave similar results, the inclusion of inhibitory neurons in the network promoted a larger portion of neurons to the role of drivers, and their properties will be investigated in the following.

## Connectivity and excitability of the driver cells

The role played by the neurons in the simulated network was elucidated by performing a directed functional connectivity (FC) analysis. In the case of the spiking network model, in order to focus on the dynamics underlying the PB build-up, the FC analysis was based on the first spike fired by each neuron in correspondence of the PBs. An equivalent information was provided in the analysis of the EC by considering the calcium signal onset to calculate the directed functional connectivity. The six neurons playing a key role in the generation of the PBs (eh₁⁻₄, ih₁⁻₂) were characterized by high values of functional out-degree, namely with an average functional degree $D^O = 68\% \pm 8\%$, ranking them among the 16 neurons with the highest functional degree. Given the high functional out-degree and their fundamental role in the generation of the PBs (as shown by the SND in Fig 3(a)), we identified these neurons as driver hub cells. The high value of $D^O$ reflected their early activation in the PB, thus preceding the activation of the majority of the other neurons.

Next, we examined the structural degree of the neurons, specifically we considered the total structural degree $K^T$, which is the sum of the in-degree and out-degree of the considered cell. As shown in Fig 3(f), we observed an anti-correlation among $D^O$ and $K^T$ where neurons with high functional connectivity are typically less structurally connected than LC neurons. This was particularly true for the six driver hubs, previously examined, since they were characterized by an average $K^T = 15 \pm 3$, well below the average structural connectivity of the neurons in the network (≈ 20).

Concerning the excitability, the six driver hubs despite being in proximity of the firing threshold (slightly above or below) as shown in S6(a) Fig, they were among the 25% fastest spiking neurons in control condition, (as shown in Fig 3(c)). In particular, the three neurons eh₁, eh₂, ih₂ were supra-threshold, while neurons eh₃, eh₄, ih₁ were slightly below the threshold. When embedded in the network their firing activity was modified, in particular three couples of neurons with similar firing rates can be identified, namely (eh₁, ih₁), (eh₂, eh₂) and (eh₃, eh₄), as reported in Table 1. The direct structural connections present among these couples (see also

<table><caption>**Table 1. Properties of driver hub cells in control condition.** For each driver hub cell (ih₁, ih₂, eh₁, eh₂, eh₃, eh₄) the columns report the intrinsic excitability ($I^b$), the average spiking frequency in control conditions ($v$), the structural out-degree ($K^O$) and in-degree ($K^I$), the functional out-degree ($D^O$) and in-degree ($D^I$).</caption><thead><tr><th>Driver hub cell</th><th><i>I</i><sup>b</sup> (mV)</th><th><i>v</i> (Hz)</th><th><i>K</i><sup>O</sup></th><th><i>K</i><sup>I</sup></th><th><i>D</i><sup>O</sup></th><th><i>D</i><sup>I</sup></th></tr></thead><tbody><tr><td>ih<sub>1</sub></td><td>14.91</td><td>19.5</td><td>11</td><td>7</td><td>76%</td><td>1%</td></tr><tr><td>ih<sub>2</sub></td><td>15.02</td><td>15.5</td><td>6</td><td>9</td><td>73%</td><td>3%</td></tr><tr><td>eh<sub>1</sub></td><td>15.32</td><td>20.5</td><td>6</td><td>5</td><td>77%</td><td>0%</td></tr><tr><td>eh<sub>2</sub></td><td>15.23</td><td>16.3</td><td>7</td><td>6</td><td>63%</td><td>11%</td></tr><tr><td>eh<sub>3</sub></td><td>14.93</td><td>1.3</td><td>7</td><td>10</td><td>59%</td><td>12%</td></tr><tr><td>eh<sub>4</sub></td><td>14.99</td><td>1.3</td><td>9</td><td>7</td><td>59%</td><td>12%</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1006551.t001
Fig 3(g)) could explain the observed firing entrainments, as discussed in details in the next subsection. When compared to the other hub neurons, the much lower activity of (*eh*₃, *eh*₄), corresponding to twice the average frequency of the PBs in control condition, was related to the fact that these two neurons fired only in correspondence of the ignition of collective events like PBs and *aborted bursts* (ABs), the latter being associated to an enhancement of the network activity but well below the threshold we fixed to detect PBs. This will become evident from the discussion reported in the subsection *Synaptic resources and population bursts*.

As already mentioned, besides the six driver hubs, the SNS experiments revealed the existence of a different set of 18 drivers, whose activation also impacted the population dynamics, although they had no influence when removed from the network and therefore they were not relevant for the PBs build up. These neurons represented in Fig 3 with squares were characterized by a low FC, namely *D*⁰ = 13% ± 15%. Therefore, we have termed them *driver LC cells* representing the ones which reproduced the behaviour of the driver LC cells identified in the EC (see Figs 1 and 2 and reference [24]). In the following we will refer to them as *el*... or *il*₁ according to the fact that they are excitatory or inhibitory neurons, respectively (note that only one LC driver was inhibitory). As shown in Fig 3(c), LC drivers were not particularly active (with firing frequencies below 1 Hz in control conditions) and in some cases they were even silent. Notably, under current stimulation they could in several cases arrest PBs or strongly reduce/ increase the activity with respect to control conditions as shown in Fig 3(b) for a specific level of current injection and also as discussed in detail in the following sections.

Compared to the driver hubs, driver LC cells had a lower degree of excitability (essentially they were all sub-threshold, see S6(a) Fig), which resulted in a later recruitment in the synchronization build up, and as a consequence in a lower functional out-degree. Therefore, driver LC cells were not necessary for the generation of the PBs, playing the role of followers in the spontaneous network synchronizations. As shown in Fig 3(f), driver LC neurons were characterized by a higher structural connectivity degree *K*<sub>*T*</sub> with respect to driver hubs, namely *K*<sub>*T*</sub> = 23 ± 3, and the most part of them were structurally targeting the driver hubs either directly (i.e. path length one) or via a LC driver (i.e. path length two, centered on a LC driver). In Fig 3(f), the two groups of drivers, hubs and LC cells, can be easily identified as two disjoint groups in the plane (*K*<sup>*T*</sup>, *D*⁰). These results indicated that driver hubs are not structural hubs, while the low functional connectivity neurons are promoted to their role of drivers due to their structural connections. This latter aspect will be exhaustively addressed in subsection *Tuning of PBs frequency upon hubs' and LC cells's stimulation*.

**Generality of the results.** With the purpose of demonstrating that the results reported here are generic, and not due to a fine tuning of the parameters of the model, we have analyzed fifteen different realizations of the network. In particular, we used the same distributions for the intrinsic excitabilities, synaptic parameters and structural connectivities. The parameter values were taken from random distributions with the same averages and standard deviations as defined in *Definition of the model in Methods*. Furthermore, in all the numerical experiments we kept fixed the size of the network (*N* = 100), the number of excitatory/inhibitory neurons (*N*<sub>e</sub> = 90, *N*<sub>i</sub> = 10), the average in-degree, and all the other constraints specified in *Definition of the model in Methods*. In six networks we found no bursting dynamics or number of bursts too small to be significant. While, in the remaining nine network PBs were always present and we could perform significant SND/SNS experiments on all the neurons in each network. This analysis allowed us to identify driver hub cells and driver LC cells in all these networks, with characteristics similar to the ones found in the network analysed in detail in the paper. In particular, we have identified for each network a number of hub cells ranging from two to eight with an average value 5 ± 2, and a number of LC cells ranging from 1 to 27 with an average value 12 ± 9 (apart a peculiar single network where we found just one hub and
one LC cell). By examining the nine networks displaying bursting dynamics we found the presence of inhibitory cells among the hubs in three cases and among the driver LC cells in six cases (with numbers ranging from one to four). As general features, we observed that driver hub cells were characterized by a high intrinsic excitability and a low structural connectivity: namely, $I^b$ was in the range [14.55: 15.42] mV (with average $15.0 \pm 0.2$ mV and $\simeq 37\%$ of the hubs supra-threshold), while the total connectivity $K^T$ was in the range [6: 31] (with average $16 \pm 3$ and a single hub with $K^T = 31$). On the other hand, the LC drivers were characterized by a low intrinsic excitability in the range [14.55: 15.03] mV (with average $14.7 \pm 0.1$ mV and $\simeq 99\%$ of the LC cells below the firing threshold) and by a high $K^T$ in the range [14: 32] (with average $22 \pm 4$ and a single driver LC cell with $K^T = 14$).

It is important to stress that in [12], for a model similar to the one here studied, it has been shown that the population dynamics as well as the sensitivity to SNS/SND are robust to wide variation in several parameters controlling the network design as well as the synaptic dynamics.

## Functional clique of excitatory and inhibitory neurons

In order to deepen the temporal relationship among neural firings leading to a PB, we examined the spikes emitted in a time window of 70 ms preceding the peak of synchronous activation (see Methods for details). The cross correlations between the timing of the first spike emitted by each driver hub neuron during the PB build up are shown in S7 Fig (Upper Sequence of Panels). The cross correlation analysis demonstrated that the sequence of activation of the neurons was $eh_1 \rightarrow ih_1 \rightarrow ih_2 \rightarrow eh_2 \rightarrow eh_3 \rightarrow eh_4$. The labeling previously assigned to these neurons reflected such an order. A common characteristic of these cells was that they had a really low functional in-degree $D^I$ as reported in Table 1 indicating that they were among the first to fire during the PB build-up. In particular, $eh_1$ had a functional in-degree $D^I$ zero, revealing that it was indeed the firing of this neuron to initialize all the bursts and therefore it could be considered as the leader of the clique.

A detailed inspection of the firing times, going beyond the first spike event, revealed the existence of more than one firing sequence leading to the collective neuronal activation: i.e. the existence of different routes to PBs. This is at variance with what found in [12] for a purely excitatory network, where only one route was present and all the PBs were preceded by the same ordered sequential activation of the most critical neurons. In particular, the neuron $eh_1$ fired twice before the PBs (see Fig 3(e)), usually in-between the firing of $eh_2$ and that of the pair ($eh_3, eh_4$), and this represents the main route, occurring for $\simeq 85\%$ of the PBs. Along the second route (present only for the $\simeq 7\%$ of the PBs), $eh_1$ was firing the second time at the end of the sequence. The neuron $eh_1$ fired essentially by following its natural period $T_1 = \tau_m \ln[(I_{eh_1}^b - v_r)/(I_{eh_1}^b - V_{th})] = 52.15$ ms, and its second occurrence in the firing sequence depended on the delay among the firing of the other neurons. As a matter of fact we verified that the elimination of the second spike emitted by $eh_1$ from the network dynamics didn't prevent, and didn't delay, the onset of the PB and had only a marginal effect on the firing of a very limited number of neurons in the PB. Therefore we can conclude that it is not essential to the PB build up. The two routes leading to the PB build-up are shown in Fig 3(e).

To observe a PB the six driver hubs should fire not only in an ordered sequence, as shown in Fig 3(e), but also with defined time delays, their average values with the associated standard deviations are reported in S1 Table for the two principal routes. These results clearly indicate that the six driver hubs are arranged in a *functional clique* whose activation was crucial for the PB build-up. In the period between the occurrence of two PBs, the driver hubs in the clique could be active, but in that case they did not show the precise sequential activation associated
to the main and secondary route, see the out-of-burst results reported in the Lower Sequence of Panels in S7 Fig. A remarkable exception is represented by the case of the ABs, in that case PBs are not triggered despite the presence of the right temporal activation of all the hubs in the clique, due to the lack of synaptic resources (as discussed in details in subsection *Synaptic resources for population bursts*). Out of PBs and ABs, we registered clear time-lagged correlations only for those neuronal pairs sharing direct structural connections (shown in Fig 3(g)): namely, $eh_1 \rightarrow ih_1$, $ih_2 \rightarrow eh_2$ and $eh_2 \rightarrow (eh_3, eh_4)$. The firing delays of these neuronal pairs were not particularly altered also out of burst with respect to those measured during the burst build-up and reported in S1 Table.

As shown in Fig 3(g), the *eh*₃ neuron represented the cornerstone of the clique, receiving the inhibitory input coming from the structural pair (*eh*₁, *ih*₁) and the excitatory one from the pair (*ih*₂, *eh*₂), with the activity of the neurons within each pair perfectly frequency locked. More specifically, *eh*₁ entrained the activity of *ih*₁ (below threshold in isolation) so that both neurons before a PB fired with a period quite similar to the natural period of *eh*₁. The other pair (*ih*₂, *eh*₂) was controlled by the inhibitory action of *ih*₂ that slowed down the activity of *eh*₂, whose natural period was 60.6 ms, while before a PB *ih*₂ and *eh*₂ both fired with a slower period, namely 72 ± 2 ms.

As it will be explained in details in the next two subsections, the two requirements to be fulfilled for the emergence of PBs are the availability of sufficient synaptic resources at neurons *eh*₃ and *eh*₄ and the coordinated activation of *eh*₁ (and *ih*₁) with the pair (*ih*₂, *eh*₂), in the absence of any synaptic connection between the two pairs.

## Synaptic resources for population bursts

Next we analyzed the relation between the evolution of synaptic resources in the hub driver cells and the onset of the PB. The availability of synaptic resources was measured by the effective efferent synaptic strength $X_{eh_3}^{OUT}$ as defined in Eq (8). In particular, we considered the available resources only for the hub neurons *eh*₃ and *eh*₄ which were the last neurons of the clique to fire before the PB ignition. We have examined only these two hub neurons, because whenever *eh*₃ and *eh*₄ fired, a burst or an AB was always delivered.

Neurons *eh*₃, *eh*₄ were receiving high frequency excitatory inputs from *eh*₂ (although the natural firing of *eh*₂ was slowed down by the incoming inhibition of *ih*₂) and high frequency inhibitory inputs from *ih*₁ (entrained by the *eh*₁, the neuron with highest firing frequency in the network). This competitive synaptic inputs resulted in a rare activation of *eh*₃ compared to the higher frequency of excitatory inputs arriving from *eh*₂. The period of occurrence of the ABs was comparable to the average interval between PBs (namely, $T_{PB} = 1.4 \pm 1.0$ s) and ABs were preceded by the sequential activations of the six critical neurons of the clique in the correct order and with the required delays to ignite a PB. The number of observed ABs was 66% of the PBs, thus explaining why the average firing period of *eh*₃ and *eh*₄ was $T_{eh_3} = 0.8 \text{ s} \approx T_{PB}/(1 + 0.66)$, since their firing always triggered a PB or an AB.

To understand why in the case of ABs the sequential activation of the neurons of the clique did not lead to a PB ignition, we examined the value of synaptic resources for regular and aborted bursts, as shown in Fig 4(a). From the figure it is clear that $X_{eh_3}^{OUT}$ and $X_{eh_4}^{OUT}$ should reach a sufficiently high value in order to observe a PB, otherwise one had an AB. Furthermore, the value reached by $X_{eh_3}^{OUT}$ and $X_{eh_4}^{OUT}$ was related to the time passed from the last collective event and thus the requirement of a minimal value of the synaptic resources to observe a PB set a minimal value for the IGI, i.e. the interval between two PBs. As a matter of fact, as shown in Fig 4(b) the IGI values grew almost linearly with the values reached by $X_{eh_3}^{OUT}$ just
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Model—Population bursts and synaptic resources. (a) Top panel: raster plot of the network activity, where population bursts (PBs) and aborted bursts (ABs) are shown. The vertical (red) dashed lines signal the occurrence of aborted burst. Bottom panel: average synaptic strength of the efferent connections of the two hub neurons *eh*3, *eh*4 in control conditions; the output effective synaptic strength is measured by the average value of the fraction *X*<sub>eh3</sub><sup>OUT</sup> (thick black line), *X*<sub>eh4</sub><sup>OUT</sup> (thin blue line) of the synaptic transmitters in the recovered state associated to the efferent synapses. The dashed horizontal lines signals the values of the local maxima of *X*<sub>eh3</sub><sup>OUT</sup> (black line, *X*<sub>eh3</sub><sup>out</sup> = 0.793) and *X*<sub>eh4</sub><sup>OUT</sup> (blue line, *X*<sub>eh4</sub><sup>out</sup> = 0.666) corresponding to the occurrence of the shortest IGI, IGI<sub>min</sub>. (b) Values of the local maxima of *X*<sub>eh3</sub><sup>OUT</sup> (max(*X*<sub>eh3</sub><sup>OUT</sup>)) in correspondence of the latest IGI. In both the figures the (green) arrow marks the occurrence of IGI<sub>min</sub>.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1006551.g004

before the PB, at least for $X_{eh3}^{OUT} < 0.9$. At larger values the relationship was no more linear and a saturation was observable, due to the fact that $X_{eh3}^{OUT}$ could not overcome one.

We could conclude that the slow firing of the couple (*eh*3, *eh*4), moderate by the inhibitory action of *ih*₁ on *eh*₃, was essential to ignite a PB, since a faster activity would not leave to the synapses the time to reach the minimal value required for a PB ignition, namely *X*<sub>eh3</sub><sup>out</sup> = 0.793 and *X*<sub>eh4</sub><sup>out</sup> = 0.666. This could be better understood by reconsidering the SND experiment on *ih*₁, as expected the resection of neuron *ih*₁ from the network led to a much higher activity of neurons *eh*₃ and *eh*₄, as shown in S8 Fig. However, this was not leading to the emission of any PBs, because in this case the value of *X*<sub>eh3</sub><sup>OUT</sup> and *X*<sub>eh4</sub><sup>OUT</sup> remained always well below the value required for a PB ignition.

Moreover, also among the network synchronization events classified as PBs we observed a variability in the number of neurons taking part to a PB. A detailed analysis reported in S3 Text reveals that at least two kind of events can be identified: one which sees the contribution of almost all the active neurons and another one where the participation is more limited (see S9 Fig). These results resemble the ones reported experimentally for the EC in [24]. Furthermore, the analysis shows a clear correlation between the available synaptic resources of the drivers controlling the population activity and the entity of the observed PBs.

## Tuning of PBs frequency upon hubs' and LC cells' stimulation

In order to better understand the role played by the hub and LC drivers for the collective dynamics of the network, we performed SNS experiments for a wide range of stimulation currents. The results of this analysis for currents in the range 14.5 mV ≤ *I*<sup>stim</sup> ≤ 18 mV are shown in Fig 5 (where all the driver hubs and six representative cases of driver LC cells are reported) and in Fig 6(a). The driver hub neurons could, upon SNS, usually lead to a reduction, or
<figure><img src="image_9.png" /><figcaption>**Fig 5.** Model—PBs frequency is tuned by current stimulation of driver cells. The plots report the number of PBs emitted during SNS of the hub neurons *eh*₁, *eh*₂, *eh*₃, *eh*₄, *ih*₁, *ih*₂ (a-f) as well as of the driver LC cells *il*₁, *el*₁, *el*₂, *el*₃, *el*₄, *el*₇ (g-l) for a wide range of the stimulation current *I*<sup>stim</sup> (over a time interval Δ*t* = 84 s). The blue vertical dashed lines, resp. the horizontal magenta solid line, refer to the value of the intrinsic excitability, resp. to the bursting activity, when the network is in control condition. The threshold value of the current is set to *V*<sub>th</sub> = 15 mV.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1006551.g005

silencing, of PBs, apart for two cells (namely, *eh*₁ and *eh*₄) which, for specific stimulation currents, could even enhance the population activity. On the other hand, the 18 driver LC cells can be divided in two classes LC1 and LC2 according to their influence on the network dynamics upon SNS: a first group of 14 driver LC1 cells able mainly to reduce/stop the collective
<figure><img src="image_10.png" /><figcaption>**Fig 6. Model—Response to SNS of the driver cells.** (a) Quantification of the response to SNS for each of the driver cells, sorted in three groups, respectively hubs, LC1 and LC2. The heatmap displays how many times the SNS over a wide range of currents (namely, 14.5 mV ≤ I<sup>stim</sup> ≤ 18 mV) induced a given number of PBs (x-axis) in the network. To facilitate the visualization, each row of the heatmap has been smoothed with a gaussian function of 1.58 standard deviation and unitary area. The red vertical lines denote the limits of activity in control condition: one standard deviation around the average. **Model—Current stimulation of driver LC cells can modify the functional clique of the network.** The panels (b),(c) refer to driver LC1 cell *el*₁ (see Fig 1(b,S)–1(d,S)). In the top panel (b) the configuration of the functional clique is reported for some sample stimulation current I<sup>stim</sup> of *el*₁. Full circles, resp. open squares, signal the presence, resp. absence, of the corresponding neuron of the functional clique *ih*₁, *ih*₂, *eh*₁, *eh*₂, *eh*₃, *eh*₄, while the open (red) circles indicate the presence of new neurons *NH* in the functional clique (the number of new neurons is reported inside the red circles). In the bottom panel (c) it is shown the number of PBs emitted by the network (black line with dots and left y-axis) and the firing frequency *v* of the LC cell (green line and right y-axis) during the current stimulation. The vertical (magenta) line marks the threshold value, *V*<sub>th</sub>, while the vertical (blue) dashed line signals the intrinsic excitability of the LC cell in control condition.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1006551.g006

activity, and in few cases to increase the PB frequency, and a group of 4 LC2 neurons capable
only to enhance the PB frequency. The three neurons *el*₁, *el*₂ and *el*₃, previously considered in
subsection *Numerical evidences of driver LC cells* for comparison with experimentally identi-
fied LC cells, belonged to the class LC1 (see Figs 1 and 2), while we have no experimental
examples of LC2 cells.

For what concerns the driver hubs’ dynamics, PBs were generated in the network whenever
the hubs *eh*₂ and *ih*₁, both structurally connected to *eh*₃, were stimulated with currents smaller
than the excitability *I*<sub>*eh*₁</sub><sup>b</sup> of the leader of the clique and within a specific interval (see Fig 5(b)
and 5(e)). This means that in order to have a PB both neurons controlling *eh*₃ should not fire
faster than the leader of the clique. If this was not the case, the inhibition (originating from *ih*₁)
would not be anymore sufficient to balance the excitation (carried by *eh*₂) or viceversa, thus
leading *eh*₃ to operate outside the narrow current window where it should be located to pro-
mote collective activity (see Fig 5(c)). In the case of *ih*₂ and *eh*₄ the SNS produced a less pro-
nounced impact on the PB activity, their stimulation could never silence the network (as
shown in Fig 5(d) and 5(f)), apart in two narrow stimulation windows for *ih*₂. This is in agree-
ment with what reported in Fig 3(a) for the SND, since the removal of neurons *eh*₄ and *ih*₂
only reduced the occurrence of PBs of ≈ 60%.
**LC drivers impact hub neurons.** The SNS of LC1 drivers could induce, in 10 cases over 14 identified LC1 cells, a complete silencing in the network. A peculiar feature of eight out of these ten cases was that the PBs were completely suppressed as soon as these LC1 driver were brought supra-threshold: two examples are reported in Fig 5(g) and 5(i). The first example in Fig 5(g) refers to the unique inhibitory driver LC cell we have identified, namely *il*₁, which was directly connected to the hub neuron *eh*₃ (as shown in Fig 3(g)). A stimulation of *il*₁ led to a decrease of the activity of *eh*₃ and as a direct consequence of the PB activity. Fig 5(i) is devoted to *el*₂, previously examined in Section *Numerical evidences of driver LC cells* and reported in Fig 1(e.S)–1(g.S). The depressive effect on the network activity due to the stimulation of *el*₂, could be straightforwardly explained by the fact that *el*₂ is directly connected to the inhibitory LC cell *il*₁ and to the inhibitory driver hub *ih*₁, thus performing an effective inhibitory action on the network, even if the stimulated driver *el*₂ was excitatory. For the other six excitatory LC1 drivers acting on *il*₁ only in two cases the PBs could not be completely blocked, and this happened when the cells were also directly connected to the driver hub *eh*₃. In the two cases of driver LC1 cells able to block the population activity, but not acting through *il*₁, these cells were exciting either *eh*₁ or *ih*₁, both belonging to the path with an inhibitory effect on *eh*₃. The four remaining LC1 drivers that were able to reduce, but not to completely silence the population activity, acted either through *eh*₄ (which was unable to block the PBs, even upon SND) or by simultaneously exciting and inhibiting *eh*₃.

It should be remarked that nine of the previously discussed LC1 drivers could either enhance or reduce the PBs for different values of $I^{\text{stim}}$. The double action of these neurons is exemplified by the two examples reported in Fig 5(h) and 5(j), which refer to neuron *el*₁ and *el*₃, already examined in connection with experimental data in Fig 1(b.S)–1(d.S) and in Fig 2(a. S)–2(f.S), respectively.

The neuron *el*₁ was structurally connected to the inhibitory hub cell *ih*₁, *el*₁ was silent in control condition and once stimulated with a current $I^{\text{stim}} = 15.135$ mV was able to enhance population dynamics as shown in Fig 1(b.S)–1(d.S). However, depending on the stimulation current it could even completely silence the PB activity, as shown in Fig 5(h). In order to understand in deeper details the mechanisms underlying both the enhancement and suppression of PBs, we stimulated *el*₁ with $I^{\text{stim}} \in [14.5: 16]$ mV and for each value of the stimulation current we performed SND of the all neurons in the network. These analyses were aimed at identifying (for each value of $I^{\text{stim}}$) the driver hub cells involved in the PB generation, i.e the neurons that upon SND reduced the population activity at least or more than 50%. The results of these experiments are shown in Fig 6(b) and 6(c), for sufficiently low stimulation currents (even above threshold) the activity of *el*₁ had no influence at a network level, and this is consistent with the response of *ih*₁ upon SNS reported in Fig 5(e). However for higher stimulation current the clique of functional hubs is modified by the action of *el*₁: not all the hub cells previously identified remained relevant for the network activity and in some cases some new driver hubs was identified, as reported in Fig 6(b). The most significant modification is that the neuron *eh*₁ was no more relevant (in most cases) for the PB generation, and this could be explained by the fact that the inhibitory hub *ih*₁ is now controlled directly by *el*₁. This is further confirmed by the fact that when the stimulation became sufficiently large the collective dynamics was completely silenced due to the high activity of *ih*₁. As a matter of fact, some lower activity in the network could be restored for even larger current values above $I^{\text{stim}} \approx 15.65$ mV. This effect can be explained in terms of the prevalence of facilitation over depression in the synapse connecting *el*₁ with *ih*₁, as discussed in S5 Text and in S12 Fig. Furthermore, at these large stimulation currents one can observe also modified functional cliques.

The LC1 driver *el*₃ had also a double action leading to enhancement or depression of the collective activity as shown in Fig 5(j). This double action grounded in the following network
architecture: *el*₃ was structurally connected, via the bridge neuron *el*₁, to *ih*₁, whose impact on the network was to arrest the bursting apart a very narrow range of stimulation currents (Fig 5(e)); *el*₃ was also structurally connected to *eh*₄, which in some ranges of stimulation currents could enhance network dynamics (Fig 5(d)).

To conclude the analysis of the driver cells, we consider LC2 cells. These were excitatory neurons characterized by a low *I*<sup>b</sup> (below the firing threshold *V*<sub>th</sub>) and, given the imposed correlation in the network model, by a high global structural connectivity (see S4 Fig). In control conditions these neurons were not active and did not participate to PBs. Two examples of the SNS of these neurons are reported in Fig 5(k) and 5(l) for LC2 drivers *el*₄ and *el*₇. Whenever they were stimulated above threshold they induced a sharp increase in the PB activity in the order of 50%. Furthermore, in the case of LC2 drivers the SNS led in general to a more regular bursting dynamics, characterized by a smaller average Inter GDP Interval (of the order of the recovery time for the synapses, see Methods for details) and a smaller standard deviation (i.e. for *el*₇ we measured <IGI> = 0.9 ± 0.6 s for *I*<sup>stim</sup> = 15.9 mV) with respect to control conditions (where <IGI> = 1.4 ± 1.0 s). When LC2 cells were current stimulated, the action of enhancement of PBs activity was not mediated by the impact on other driver cells. As shown in S10 Fig, while the stimulation of LC1 drivers led to a noticeable modification of the firing rates of the hub cells, the SNS of driver LC2 cells had essentially no influence on the hubs. Therefore due to their high structural out-degree we can safely affirm that their influence on the network dynamics should be related to a cooperative excitatory effect. As a matter of fact, during the SNS of driver LC2 cells we observed the disappearance of ABs and as a consequence Inter GDP Interval become more regular, thus leading to an enhancement of the population activity. In particular, the disappearance of ABs was due to the fact that whenever *eh*₃ was firing, a burst was emitted due to the presence of a higher level of excitation in the network, even when the synaptic resources of *eh*₃ were below the minimal value required in control conditions, as discussed in the subsection *Synaptic resources and population bursts*.

Discussion

We have developed a simple brain circuit model to mimic recent experimental results obtained in cortical slices of the mice Entorinhal Cortex during developmental stages [24]. These analysis revealed the existence of high and low functionally connected driver cells able to control the network dynamics. The fact that functional hubs can orchestrate the network dynamics is somehow expected [17, 44], while the existence of driver neurons with low functional outdegrees has been revealed for the first time in [24]. In this paper, we focused mainly on the analysis of these latter class of drivers, which in control conditions were essentially irrelevant for the build-up of the GDPs. On the contrary, if single-handedly stimulated they could nevertheless strongly modify the frequency of occurrence of GDPs, as evident from the experimental findings reported in Figs 1(b)–1(g.E) and 2(a)–2(f.E). In particular, their stimulation could lead both to an enhancement as well as to a reduction of the population activity (GDPs’ frequency). Quite remarkably, some of the driver LC cell were able to perform both these tasks as an effect of different stimulation frequencies as revealed by the experiment shown in Fig 2(a)–2(f.E).

Furthermore we have demonstrated that the experimental findings could be replicated in a
simple spiking neural network model made of excitatory and inhibitory neurons with short-
term synaptic plasticity and developmentally inspired correlations (see Figs 1(b)–1(g.S) and
2(a)–2(f.S)). The analysis of the model has allowed to understand the fundamental mecha-
nisms able to promote a single neuron to the role of network driver without being a functional
hub, as usually expected.
In the model, all the driver neurons able to influence the network dynamics could be identified and they could be distinguished in neuronal hubs characterized by high out-degree or low functionally connected drivers. Functional hubs are highly excitable excitatory and inhibitory neurons arranged in a clique, whose sequential activation triggered the Population Bursts (analogous to GDPs). This in agreement with recent experimental evidences that small neuronal groups of highly active neurons can impact and control cortical dynamics [7-10]. On the other hand, driver LC cells are characterized by a lower level of excitability, but a higher structural connectivity with respect to driver hubs. Due to their low activity and functional connectivity in control conditions, these neurons were not fundamental for the PBs development, but were passively recruited during the burst, or even completely silent. The LC drivers can be divided in two classes LC1 and LC2 according to their influence on the population activity whenever stimulated with different values of DC current: the majority of them were able both to enhance and reduce (or even set to zero) the frequency of occurrence of the PBs (LC1), while a small group was able only to enhance the PBs' frequency with respect to control conditions (LC2). Noticeably, driver LC1 cells were structurally connected to the hubs (directly or via a bridge LC cell). Therefore, whenever stimulated they can influence the network activity by acting on the clique dynamics. In most cases, even if these cells were excitatory, their action on the network was mainly depressive, since either they stimulated directly inhibitory hubs or the inhibitory LC1 driver, which acted as a bridge over the clique. In more than the 50% of the cases (8 over 14) whenever brought over threshold driver LC1 cells led to a complete arrest of the PB activity.

Driver LC2 cells instead were silent in control conditions and highly structurally connected, therefore they were putative structural hubs. As a matter of fact, whenever brought supra threshold they favoured a more regular collective dynamics. The activation of the many efferent connections of LC2 drivers led to the creation of many alternative pathways for the PB igniton, in a sort of homeostatic regulation of the network which led to an optimal employ of the synaptic resources [45] with the corresponding disappearance of the aborted bursts, largely present in control conditions.

Furthermore, we have shown that the stimulation of single driver LC cells was not only able to alter the collective activity but also to deeply modify the role of neurons in the network, such that some neurons can be promoted to the role of driver hubs or driver hubs can even lose their role (see Fig 6(b) and 6(c)). At variance with purely excitatory networks [12], the synchronized dynamics of the present network, composed of excitatory and inhibitory neurons, is less vulnerable to targeted attacks to the hubs [46, 47]. As demonstrated by the fact that different firing sequences of hub neurons can lead to population burst ignitions (see Fig 3(e)) and that hubs can be easily substituted in their role by driver LC cells when properly stimulated. The robustness of the synchronized dynamics is confirmed by the fact that the presence of channel noise, up to quite large noise strength, does not substantially modify the composition neither of the functional clique nor of the LC drivers (for more details on the analysis see S4 Text and S11 Fig).

The experimental analysis of the population activity in the EC has revealed the existence of two kind of synchronization events: global, where most of the neurons were involved, and local, where only a sub-group took part to the burst [24]. In our model we have shown a similar effect that can be explained in terms of the available synaptic resources of the driver cells (for more details see S3 Text and S9 Fig). It has been reported in different contexts [48, 49] that the strength of the connectivity between neurons is a key variable defining neuronal ensemble dynamics. However, whether the mechanism that we have found to be at the basis of this variability could explain also the experimental findings in [24] is an interesting question
that cannot be answered with our experimental set-up, and that we leave open for future investigation.

Another relevant aspect is that the inclusion of inhibitory neurons in the network did not cause a trivial depressing action on the bursting activity, as it could be naively expected, but instead they can play an active role in the PB build-up. Our analysis clearly demonstrate that their presence among the driver cells is crucial in determining and controlling the PB activity, somehow similarly to what found in [20] where it has been shown that the emergence of sharp-wave in adult hippocampal slices was controlled by single perisomatic-targeting inter-neurons. We expect that the model described in this paper could also find application for other developing circuits, beyond EC, where inhibitory GABAergic synapses are also present.

Our results suggest that inhibitory neurons can have a major role in information encoding by rendering on one side the population dynamics more robust to perturbations of input stimuli and on another side much richer in terms of possible repertoire of neuronal firings. These indications confirm the key role of inhibitory neurons in neural dynamics, already demonstrated for the generation of brain rhythms [50, 51] and for attentional modulation [52].

Entorhinal cortex is involved in human mesial temporal lobe epilepsy [53]. A recent in vitro study on the onset of seizure-like events in EC slices [54] revealed fast-rising and sustained extracellular potassium increases associated to interneuronal network activity consistently preceding the initiation of seizures, supporting a key role of interneuron activity in the EC in focal seizure generation. There results are consistent with our observation that inhibitory neurons play a key role in the generation and orchestration of network synchronizations.

The main ingredient introduced in our model to mimic the developmental stage was an anti-correlation between intrinsic excitability and the synaptic connectivity inspired by homeostatic regulation mechanisms observed during neural development [32]. Interestingly a similar phenomenon co-regulating neuronal connectivity and excitability could possibly occur in adult circuits due to functional switches or functional losses (such as those following an injury). In such cases the prolonged absence of functional input in neurons could synergistically lead to synaptic depletion and increased neural excitability. Therefore, massive pre-synaptic neural loss in adults circuits could as a consequence enhance the neural excitability of the post-synaptic targets [9], possibly upgrading the role of these cells to functional hubs with eventual pathological consequences as we could imagine in the case of epilepsy.

Recently there has been a renovated interest on the existence and role of neuronal cliques within the brain circuitries [13, 55]. Cliques have been proposed as structural functional multi-scale computational units whose hierarchical organization can lead to increasingly complex and specialized brain functions and can ground memory formations [55]. In addition, activation of neuronal cliques as in response to external stimuli or feedforward excitation can lead to a cascade of neuronal network synchronizations with distinct spatio-temporal profiles [13]. Our results provide a further understanding on how cliques can emerge (spontaneously during development) and modify (in response to stimuli similarly to the SNS here discussed) with a consequent reshaping of the spatio-temporal profile of the dynamics of the network in which the clique is embedded. Notably, it is the presence of inhibitory neurons within the network to favour the emergence of different cliques by empowering drivers with different functional connectivity degree. While driver functional hubs guarantee the functioning of the network synchronization in absence of stimuli (such as during development and in non-stimulated conditions), LC drivers widen the ability of the network to play distinct synchronization profiles (i.e. spatio-temporal activations) possibly underlying emergent functions within the brain networks.

Finally, our results could be of some relevance also for the control of collective dynamics in complex networks [56, 57]. Usually the controllability of complex networks is addressed with
linear dynamics [58, 59]. However, at present there is not a general framework to address controllability in nonlinear systems, in this context the SND and SNS protocols we developed for pulse-coupled networks could be extended to general complex networks as a tool to classify driver nodes and as a measure of controllability [60]. Control methods based on linear stability analysis for seizure propagation on large-scale brain networks have been recently reported [61]. Our methodology based on the SNS and SND protocols could be useful to go beyond this approach in order to single out the epileptogenic zone and the brain areas involved in the propagation of the seizure in patient specific brain models [62].

# Methods

## Experiment

**Ethics statement.** All animal use protocols were performed under the guidelines of the French National Ethic Committee for Sciences and Health report on “Ethical Principles for Animal Experimentation” in agreement with the European Community Directive 86/609/EEC.

**Animal treatment.** Double-homozygous Mash1BACCreER/CreER/RCE:LoxP+/+ and Dlx1/2CreER/CreER/ RCE:LoxP+/+ [63, 64] male mice were crossed with 7 to 8-week-old wild-type Swiss females (C.E Janvier, France) for offspring production. To induce CreER activity, we administered a tamoxifen solution (Sigma, St. Louis, MO) by gavaging (force-feeding) pregnant mice with a silicon-protected needle (Fine Science Tools, Foster City, CA).

**Slice preparation and calcium imaging.** Horizontal cortical slices (400 mm thick) were prepared from 8 day old (P8) GAD67<sup>GFP</sup> (n = 29), Lhx6<sup>iCre::</sup>RCE:LoxP (n = 23) or 5-HT3aR-BAC<sup>EGFP</sup> (n = 15) mouse pups with a Leica VT1200 S vibratome using the Vibrocheck module in ice-cold oxygenated modified artificial cerebrospinal fluid (0.5 mM CaCl<sub>2</sub> and 7 mM MgSO<sub>4</sub>; NaCl replaced by an equimolar concentration of choline). Slices were then transferred for rest (1 hr) in oxygenated normal ACSF containing (in mM): 126 NaCl, 3.5 KCl, 1.2 NaH<sub>2</sub>PO<sub>4</sub>, 26 NaHCO<sub>3</sub>, 1.3 MgCl<sub>2</sub>, 2.0 CaCl<sub>2</sub>, and 10 D-glucose, pH 7.4. For AM-loading, slices were incubated in a small vial containing 2.5 ml of oxygenated ACSF with 25 ml of a 1 mM Fura2-AM solution (in 100% DMSO) for 20–30 min. Slices were incubated in the dark, and the incubation solution was maintained at 35–37°C. Slices were perfused with continuously aerated (3 ml/min; O<sub>2</sub>/CO<sub>2</sub>-95/5%) normal ACSF at 35–37 °C. Imaging was performed with a multibeam multiphoton pulsed laser scanning system (LaVision Biotech) coupled to a microscope as previously described (see [65]). Cell activity is expressed as relative intensity changes in the fluorescence transients of the FURA indicator. Images were acquired through a CCD camera, which typically resulted in a time resolution of 50–150 ms per frame. Slices were imaged using a 20X, NA 0.95 objective (Olympus). Imaging depth was on average 80 µm below the surface (range: 50–100 µm).

**Experimental design.** We considered the same data set analyzed in [24] composed of n = 64 neurons, plus the three cells examined in Figs 1 and 2. Therefore, a total of n = 67 neurons were electrophysiologically stimulated and recorded following the criteria: (1) stable electrophysiological recordings at resting membrane potential (i.e., the holding current did not change by more than 15 pA); (2) stable network dynamics measured with calcium imaging (i.e., the coefficient of variation of the inter-GDP interval did not exceed 1); (3) complete labeling of the recorded cell; and (4) good quality calcium imaging while recording. Neurons were held in current-clamp using a patch-clamp amplifier (HEKA, EPC10) in the whole-cell configuration. Intracellular solution composition was (in mM): 130 K-methylSO<sub>4</sub>, 5 KCl, 5 NaCl, 10 HEPES, 2.5 Mg-ATP, 0.3 GTP, and 0.5% neurobiotin. No correction for liquid junction potential was applied. The osmolarity was 265–275 mOsm, pH 7.3. Microelectrodes resistance
was 6–8 MOhms. Uncompensated access resistance was monitored throughout the recordings. Recordings were digitized online (10 kHz) with an interface card to a personal computer and acquired using Axoscope 7.0 software. Spontaneous EPSPs were detected and analyzed using the MiniAnalysis software. A measurement of the average period of occurrence of GDPs was usually performed even before patching the neuron. For most stimulation experiments, the movie acquisition time was separated evenly in three epochs: (1) a 2 min resting period during which the cell was held close to $V_{rest}$ (i.e., zero current injection); (2) a 2 min stimulation period during which phasic stimulation protocols were applied; and (3) a 2 min recovery period where the cell was brought back to resting membrane potential. Stimulation protocol: supra-threshold current pulses (amplitude: 100–200 pA, duration: 200 ms) repeated at one of the following frequencies: $v_S = \{1/IGI, 0.5Hz, 1Hz\}$, where $IGI = 4.9 \pm 1.5$ s is the average GDP period in the phase before patching the considered neuron. The chosen values are equal or larger than the average GDP frequency to maximize the impact of the stimulation on the collective dynamics, slower frequency of stimulation would not allow to interact with all the GDP events.

**Analysis of the data.** *Signal detection.* We used custom designed MATLAB software [17] that allowed: (1) automatic identification of loaded cells; (2) measuring and normalizing between 0 and 1 the average fluorescence transients from each cell as a function of time (FURA signal); (3) detecting the onsets and offsets of calcium signals, and (4) reconstructing the functional connectivity of the imaged network.

**Statistical analysis.** Network synchronizations (GDPs) were detected as synchronous onsets peaks including more neurons than expected by chance, and their time stamp denoted by $t_G$. The threshold is defined as the 95th percentile of peak heights after reshuffling 1000 times the cell fluorescence traces. The Inter-GDP-interval (IGI) is defined as the interval between two consecutive GDPs. To establish whether the stimulation of a single neuron was able to influence the frequency of GDPs occurrence, we first calculate the average IGI in the three epochs: pre-stimulaion (control), during the stimulation period, and post-stimulation (recovery). Due to the variability distribution of IGI in each interval, we calculate, the average IGI in a window of $t_s = 60s$ calculated starting from each $t_G$, eliminating the data corresponding to overlaps between epochs. To test for the significance of the change in the period of GDP due to single neuron stimulation a Kolmogorov-Smirnov test is applied between all the 3 resulting distributions of average IGI and a significance level of $p < 0.05$ is chosen to be realiable.

Moreover, a phase shift indicator $\Phi_i$ that tracks the change in frequency of the GDP with respect to the control conditions is defined as follows:

$$ \Phi_i = (t_G^i - \langle t_G^i \rangle)/\Delta t \qquad (1) $$

where $\Delta t$ is the average IGI interval measured in control conditions, and $\langle t_G^i \rangle = i \times \Delta t$ is the expected time of occurrence of the i-th GDP, according to the control condition. A meaningful definition of $\Phi_i$ requires that the process of GDP's emission is stationary.

## Model

**Definition of the model.** To study the response of bursting neural networks to single neuron stimulation and removal, we employed the Tsodyks-Uziel-Markram (TUM) model [36]. 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 [38, 39, 66]. We have considered such a model to mimic the dynamics of developing brain circuits, which is characterized by coherent bursting activities, such as giant depolarizing potentials [17, 67]. These coherent oscillations emerge, instead of abnormal synchronization,
despite the fact that the GABA transmitter has essentially an excitatory effect on immature neurons [30].

In this paper we consider a network of N leaky-integrate-and-fire (LIF) neurons interacting via synaptic currents regulated by short-term synaptic plasticity (depression and facilitation) according to the model introduced in [36]. In particular the facilitation mechanism is present only for synapses targeting inhibitory neurons.

The time evolution of the membrane potential $V_i$ of each neuron reads as

$$ \tau_m \dot{V}_i = -V_i + I_i^{\text{syn}} + I_i^{\text{b}} \quad (2) $$

where $\tau_m$ is the membrane time constant, $I_i^{\text{syn}}$ is the synaptic current received by neuron i from all its presynaptic inputs and $I_i^{\text{b}}$ represents its level of intrinsic excitability. By following the original definitions of the TUM model [36], 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 $\delta$-like function of time. Accordingly, the spike-train $S_j(t)$ produced by neuron j, is defined as,

$$ S_j(t) = \sum_m \delta(t - t_j(m)), \quad (3) $$

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 [68] the state of the synaptic connection between the jth presynaptic neuron and the ith 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} + Y_{ij} + Z_{ij} = 1$. The evolution equations for these variables read as

$$ \dot{Y}_{ij} = -\frac{Y_{ij}}{T_{ij}^{\text{f}}} + u_{ij}X_{ij}S_j \quad (4) $$

$$ \dot{Z}_{ij} = \frac{Y_{ij}}{T_{ij}^{\text{f}}} - \frac{Z_{ij}}{T_{ij}^{\text{R}}}. \quad (5) $$

Only the active transmitters react to the incoming spikes $S_j$: the adimensional parameters $u_{ij}$ tune their effectiveness. For the synapses targeting excitatory neurons, therefore displaying only synaptic depression, $u_{ij} \equiv U_{ij}$ stay constant, while in presence of facilitation, i.e. for the synapses targeting inhibitory neurons, $u_{ij}$ have the following dynamical evolution

$$ \dot{u}_{ij} = \frac{-(u_{ij} - U_{ij})}{T_{ij}^{\text{f}}} + U_{ij}(1 - u_{ij})S_j; \quad (6) $$

where $\{T_{ij}^{\text{f}}\}$ control the decay of the facilitation variables. Moreover, $\{T_{ij}^{\text{f}}\}$ represent the characteristic decay times of the postsynaptic current, while $\{T_{ij}^{\text{R}}\}$ 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

$$ I_i^{\text{syn}} = \frac{1}{K_i^l} \sum_j G_{ij} Y_{ij}, \quad (7) $$
where $G_{ij}$ are the coupling strengths, whose values can be finite (zero) if the presynaptic neuron j is connected to (disconnected from) the postsynaptic neuron i. Furthermore, if the presynaptic neuron is excitatory (inhibitory) the sign of $G_{ij}$ will be positive (negative).

In this paper, we consider a *diluted* network made of $N = N_e + N_i = 100$ neurons, where $N_e = 90$ ($N_i = 10$) is the number of excitatory (inhibitory) cells. The i-th neuron has $K_i^I$ ($K_i^O$) afferent (efferent) synaptic connections distributed as in a directed Erdös-Rényi graph with average in-degree $\bar{K}^I = 10$, as a matter of fact also the average out-degree was $\bar{K}^O = 10$. The sum appearing in (7) is normalized by the input degree $K_i^I$ to ensure homeostatic synaptic inputs [69, 70].

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_i^{OUT}$ ($X_i^{IN}$) in the recovered state associated to its efferent (afferent) synapses, namely

$$X_i^{OUT} = \frac{1}{K_i^O} \sum_k \epsilon_{ki} X_{ki} \quad , \quad X_i^{IN} = \frac{1}{K_i^I} \sum_j \epsilon_{ij} X_{ij} \qquad (8)$$

where $\epsilon_{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.

The intrinsic excitabilities of the single neurons {$I_i^b$} are randomly chosen from a flat distribution of width 0.45 mV centered around the value $V_{th} = 15$ mV, with the constraint that 10% of neurons are above threshold. This requirement was needed to obtain bursting behavior in the network. With this choice the distribution of the single neuron firing rates under control conditions is in the range [0.05;22] Hz.

Furthermore, we have considered networks where a negative correlation between the intrinsic neuronal excitability $I_i^b$ and the total connectivity (in-degree plus out-degree) $K_i^T = K_i^I + K_i^O$ is embedded. 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_i^T$, thus to ensure an inverse correlation between $I_i^b$ and $K_i^T$. The correlation is visualized in S4 Fig.

For the other parameters, we use the following set of values: $\tau_m = 30$ ms, $V_r = 13.5$ mV, $V_{th} = 15$ mV. The synaptic parameters {$T_{ij}^I$}, {$T_{ij}^R$}, {$T_{ij}^F$}, {$U_{ij}$} and {$G_{ij}$} are Gaussian distributed with averages $\bar{T}^I = 3$ ms, $\bar{T}_{ee}^R = \bar{T}_{ei}^R = 800$ ms, $\bar{T}_{ii}^R = \bar{T}_{ie}^R = 100$ ms, $\bar{T}_{ii}^F = \bar{T}_{ie}^F = 1000$ ms, $\bar{U}_{ee} = \bar{U}_{ei} = 0.5$, $\bar{U}_{ii} = \bar{U}_{ie} = 0.04$, and $\bar{G}_{ee} = 45$ mV, $\bar{G}_{ei} = 135$ mV, $\bar{G}_{ii} = \bar{G}_{ie} = 180$ mV respectively, and with standard deviation equal to the half of the average. These parameter values are analogous to the ones employed in [36] and have a phenomenological origin.

In order to have an accurate and fast integration scheme, we transformed the set of ordinary differential Eqs (2), (4), (5) and (6) into an event-driven map [71] ruling the evolution of the network from a spike emission to the next one. 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.

A simple way to include channel noise in the single neuron LIF model is to choose randomly the reset value $V_r$ whenever the neuron fires [72]. This approach guarantees that we can still perform exact event-driven integration of the model even in presence of noise. In particular, we select the random reset value for each neuron from a uniform distribution defined on the interval $[V_r - \Delta:V_r + \Delta]$, where the parameter $\Delta$ tunes the strength of the noise.

**Population bursts and aborted 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 70 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 only 17–20% of neurons were typically fir- ing, while ≈80% of the neuronal population contributed to the bursting event.

The average interburst interval for the network with (without) correlations under control conditions was 1.4 ± 1 s (0.3 ± 0.1 s), while the burst duration was 17 ± 3 ms (30 ± 6 ms) for a network made of N = 100 neurons.

Aborted bursts were collective events associated to an observable enhancement of the net- work activity, but well below the threshold we fixed to detect PBs. The period of occurrence of the ABs was comparable to the average interburst interval, while their number corresponded to the 66% of the PBs in the network with correlations and to the 47% in random networks.

**Functional connectivity.** In order to highlight statistically significant time-lagged activa- tions 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 [17, 73] 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:

$$C_{ab}(\tau) = \frac{\sum_{t=t}^{T-\tau} a_{t+\tau} b_t}{min(\sum_{i=1}^{T} a_i, \sum_{k=1}^{T} b_k)} \quad (9)$$

where {$a_t$},{$b_t$} represented the considered time series and T was their total duration. When- ever $C_{ab}(\tau)$ presented a maximum at some finite time value $\tau_{max}$ a functional connection was assigned between the two neurons: for $\tau_{max} < 0$ ($\tau_{max} > 0$) directed from a to b (from b to a). A directed functional connection cannot be defined for an uniform cross-correlation corre- sponding to uncorrelated neurons or for synchronous firing of the two neurons associated to a Gaussian $C_{ab}(\tau)$ 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 func- tional out-degree $D_i^o$ (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 occur- ring 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_i^o$, 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 emit- ted within a window of 70 ms preceding the PB peak was taken into account;
3. for the analysis of the network activity during inter-burst periods, all action potentials emit-
ted out of the PBs were considered.

**Supporting information**

S1 Text. Experiment—Effect of cell stimulation on the GDP width.

(PDF)

S2 Text. Experiment—Dependence of the results on the stimulation frequency.

(PDF)

S3 Text. Model—Population Burst variability.

(PDF)

S4 Text. Model—Robustness of the results to noise.

(PDF)

S5 Text. Model—Interplay between depression and facilitation for the control of network synchronization.

(PDF)

**S1 Fig. Experimental setup.** (a) Slice of Entorhinal cortex with calcium indicator. Contoured cells are the active cells. (b) Single neuron activity as a function of time during the three phases of the stimulation protocol: 1) pre-stimulation period where only spontaneous activity is recorded (2 min.); 2) single neuron is injected with pulses of fixed amplitude at a certain frequency $v_S$ (2 min.); 3) post-stimulation period without stimulation (2 min.). (c) Calcium trace for a selected neuron during the whole protocol. A time point is plotted in the upper part of the calcium trace whenever an onset of activity is present. Red (blue) traces denotes stimulation (control) epochs.

(EPS)

**S2 Fig. Experiment—GDP width distribution.** a) Average time trace of the negative FURA Intensity used for peak detection as a proxy of population activity indicating the three epochs of the experimental protocol as in S1 Fig. Peaks of activity (GDPs) are denoted with triangles and the widths of the GDP are indicated with black horizontal segments. b) Boxplot for the distribution of the GDP widths for the time trace in a) showing no significant differences between the three periods (KS-test). c) Distribution of the pooled data for the GDP widths during the three periods considered in each experiment.

(EPS)

**S3 Fig. Experiment—Robustness of the results with respect to $v_S$.** Paired t-test comparing IGI in control conditions and during stimulation of single neuron for each considered stimulation frequency. No significant trends are observable in any of the three cases, this observation is supported by the reported *p*-value.

(EPS)

**S4 Fig. Model—Setup for connectivity and excitability.** (a) Negative correlation between intrinsic excitability $I^b$ and total connectivity $K^T$. The (magenta) line indicates the threshold value $V_{th}$ = 15 mV, dividing supra-threshold from sub-threshold neurons. (b) Scatter plot of the in-degrees and out-degrees for each neuron in the network (no correlation). In both the figures dots (asterisks) refer to excitatory (inhibitory) neurons. The data refer to $N$ = 100 and all the parameter values are defined as in *Methods*.

(EPS)
**S5 Fig. Model—Response of the network without correlations to single neuron deletion (SND) and stimulation (SNS).** (a), (b) Number of PBs recorded during SND (SNS) experiments versus the label of the removed (stimulated) neurons, ordered accordingly to their average firing rates $v$ under control conditions (shown in panel c)). During SNS experiments each neuron was stimulated with a DC step $I^{\text{stim}}$ mV for a time interval $\Delta t = 84$ s. The horizontal dashed line shows the average number of PBs emitted in control conditions within an interval $\Delta t = 84$ s, while the horizontal dotted lines mark the 50% variation. The vertical dashed red line separates firing neurons (on the right side) from silent neurons (on the left side) in control conditions. In all the panels, dots (asterisks) symbols refer to excitatory (inhibitory) neurons. (EPS)

**S6 Fig. Model—Structural properties of the neurons.** Scatter plots showing the structural properties of the neurons of the network in control conditions, (a) intrinsic excitability $I^b$, (b) total structural connectivity $K^T$. Dots (asterisk) symbols refer to excitatory (inhibitory) neurons. The critical neurons $ih_1, ih_2, eh_1, eh_2, eh_3, eh_4$, belonging to the functional clique responsible for the PB-build up, are signaled by open circles, while the driver LC cells are denoted by open squares, in both cases red (green) contour codes for excitatory (inhibitory) neurons. The vertical dashed line separates firing neurons (on the right side) from silent neurons (on the left side) in control conditions, while the horizontal (magenta) line marks the threshold value, $V_{th} = 15$ mV, dividing supra-threshold from sub-threshold neurons. The neurons are ordered accordingly to their average firing rate in control conditions.
(EPS)

**S7 Fig. Model—The activity of driver hub cells.** Cross correlation functions between the driver hubs. The blue histograms are calculated using the first spike fired by each neuron during the PBs build-up. The red histograms are calculated using the spikes fired out of the PBs and the ABs. Note that during the PB build-up, neurons activate reliably in the following order $eh_1 \rightarrow ih_1 \rightarrow ih_2 \rightarrow eh_2 \rightarrow eh_3 \rightarrow eh_4$. During the out-of-burst activity, identical time lagged activation are preserved among the structurally connected pairs, namely $eh_1 \rightarrow ih_1$, $ih_2 \rightarrow eh_2$ and $eh_2 \rightarrow (eh_3, eh_4)$.
(EPS)

**S8 Fig. Model—Deletion (SND) of the inhibitory neuron $ih_1$ of the clique leads to the arrest of the bursting activity.** In the top panel raster plot of the network activity during SND experiment on $ih_1$ (neurons are labeled accordingly to the natural index); dots (asterisks) symbols refer to excitatory (inhibitory) neurons, while large dots and dashed (red) lines refer to the driver hubs $eh_3, eh_4$. Bottom panel: average synaptic strength of the efferent connections of the two driver hubs neurons $eh_3, eh_4$; the output effective synaptic strength is measured by the average value of the fraction $X_{eh_3}^{\text{out}}$ (black line), $X_{eh_4}^{\text{out}}$ (blue line) of the synaptic transmitters in the recovered state associated to the efferent synapses. The output effective synaptic strengths are always under the minimal values for PB ignition $X_{eh_3}^{\text{out}*}, X_{eh_1}^{\text{out}*}$ represented by the dashed lines (see also Fig 4 in the main text).
(EPS)

**S9 Fig. Model—Population Burst variability.** (a-d) Population activity in sample experiments (for the employed protocol see the main text), the time trace associated to the stimulation period is denoted in red. (b-e) Similarity matrices for the PBs showing the emergence of two clusters of events: those with high participation (denoted by circles in (a-d)) and the ones with low participation (denoted by asterisks in (a-d)). (c-f) Average value of the fraction $X^{\text{OUT}}$ of the synaptic transmitters in the recovery phase associated to the efferent synapses. (g) Ratio
$N_H/N_L$ as a function of the average PB frequency showing a high negative rank-correlation (Spearman rank). In this panel, results for drivers LC and hub cells are reported as blue and red symbols, respectively. Panels (a-c) and (d-f) correspond to the driver cells $el_2$ and $el_1$, respectively, for the same stimulation protocol employed in Fig 1(b)-1(d.S) and 1(e)-1(g.S) in the main text.

(EPS)

**S10 Fig. Model—Results of the SNS of LC drivers on the firing activity of the neurons of the clique.** Firing frequency $v$ of the driver hub neurons of the clique versus the current stimulation $I^{\text{stim}}$ during SNS of LC1 driver $el_1$ (a) and LC2 driver $el_7$ (b): (brown) stars, (red) crosses, (maroon) squares, (black) points, (green) diamonds, (blue) triangles refer respectively to $eh_1$, $eh_2$, $eh_3$, $eh_4$, $ih_1$, $ih_2$. The vertical (magenta) line marks the threshold value $V_{\text{th}}$, while the (black) arrows signal the firing frequency of the neurons of the clique in control condition $\bar{v}_{eh_1}, \bar{v}_{eh_2}, \bar{v}_{eh_3}, \bar{v}_{eh_4}, \bar{v}_{ih_1}, \bar{v}_{ih_2}$.

(EPS)

**S11 Fig. Model—Drivers in presence of noise.** In the top panel (a) the configuration of the functional clique is reported for increasing strength of the noise $\Delta$. Full circles, resp. open squares, signal the presence, resp. the absence, of the corresponding neuron of the functional clique (namely, $ih_1$, $ih_2$, $eh_1$, $eh_2$, $eh_3$, $eh_4$) found in absence of noise, while the red open circle indicates the presence of new neurons NH in the functional clique (whose number is reported inside the circle). In the bottom panel (b) it is shown the total number of LC drivers (black diamonds) identified as a function of $\Delta$ and the number of LC drivers (identified in absence of noise) which survive for finite $\Delta$ (red triangles).

(EPS)

**S12 Fig. Model—Interplay between depression and facilitation for the control of network synchronization.** In the figure it is shown the impact of the current stimulation of the LC driver $el_1$ on the network bursting mediated by short term synaptic plasticity between $el_1$ and the driver hub cell $ih_1$. a) Time average of the synaptic variables $u_{ih_1,el_1}$ (black circles) and $X_{ih_1,el_1}$ (red squares). b) Time average of the synaptic variable $Y_{ih_1,el_1}$. c) Firing frequency of the driver hub neuron $ih_1$. d) Firing frequency of the driver hub neuron $eh_3$. e) Number of PBs emitted in the network. All the measures reported in the panels refer to a time window of $\Delta t = 48$ s.

(EPS)

**S1 Table. Model—Routes leading to PBs.** Spike time delays $\Delta T$ between two successive firing of the neurons forming the functional clique along the main and secondary route leading to bursting. Neurons $eh_3$ and $eh_4$ are assumed to fire essentially at the same time, since $eh_4$ fires almost immediately after $eh_3$ within 0.03 – 0.04 ms. Notice that $eh_1$ is the only driver hub cell firing twice before a PB: the two routes could be distinguished by the time occurrence of the second spike of $eh_1$ and this event is denoted by an asterisk in the table. The arrows indicate the order of firing in the sequence.

(PDF)

## Acknowledgments

The authors had useful interactions with M. Di Bernardo, S. Olmi, M. Timme. This work has been realized within the activity of the Joint Italian-Israeli Laboratory on Integrative Network Neuroscience financed by the Italian Ministry of Foreign Affairs (S.L., P.B., A.T.).
## Author Contributions

**Conceptualization:** Arnaud Malvache, Paolo Bonifazi, Alessandro Torcini.

**Data curation:** David Angulo-Garcia, Rosa Cossart, Arnaud Malvache, Laura Módol, Vitor Hugo Sousa, Paolo Bonifazi.

**Funding acquisition:** Paolo Bonifazi, Alessandro Torcini.

**Investigation:** Stefano Luccioli, David Angulo-Garcia, Rosa Cossart, Arnaud Malvache, Alessandro Torcini.

**Methodology:** Stefano Luccioli, David Angulo-Garcia, Arnaud Malvache, Laura Módol, Paolo Bonifazi, Alessandro Torcini.

**Software:** Stefano Luccioli.

**Supervision:** Rosa Cossart, Paolo Bonifazi, Alessandro Torcini.

**Validation:** Stefano Luccioli, Paolo Bonifazi, Alessandro Torcini.

**Visualization:** David Angulo-Garcia.

**Writing – original draft:** Stefano Luccioli, Paolo Bonifazi, Alessandro Torcini.

**Writing – review & editing:** Paolo Bonifazi, Alessandro Torcini.

## References

1. Bragin A, Engel J Jr, Staba RJ. High-frequency oscillations in epileptic brain. Current opinion in neurol- ogy. 2010; 23(2):151. https://doi.org/10.1097/WCO.0b013e3283373ac8 PMID: 20160649

2. Song S, Sjöström PJ, Reigl M, Nelson S, Chklovskii DB. Highly nonrandom features of synaptic connec- tivity in local cortical circuits. PLoS biology. 2005; 3(3):e68. https://doi.org/10.1371/journal.pbio. 0030068 PMID: 15737062

3. Yoshimura Y, Dantzker JL, Callaway EM. Excitatory cortical neurons form fine-scale functional net- works. Nature. 2005; 433(7028):868–873. https://doi.org/10.1038/nature03252 PMID: 15729343

4. Kampa BM, Letzkus JJ, Stuart GJ. Cortical feed-forward networks for binding different streams of sen- sory information. Nature neuroscience. 2006; 9(12):1472–1473. https://doi.org/10.1038/nn1798 PMID: 17099707

5. Following [12], we use the term clique in an informal way to indicate a group of neurons whose activation is essential for the network dynamics. It should be stressed that the term clique here is not indicating a fully connected undirected sub-network as in graph theory.

6. Lin L, Osan R, Shoham S, Jin W, Zuo W, Tsien JZ. Identification of network-level coding units for real- time representation of episodic experiences in the hippocampus. Proceedings of the National Academy of Sciences. 2005; 102(17):6125–6130. https://doi.org/10.1073/pnas.0408233102

7. Yassin L, Benedetti BL, Jouhanneau JS, Wen JA, Poulet JF, Barth AL. An embedded subnetwork of highly active neurons in the neocortex. Neuron. 2010; 68(6):1043–1050. https://doi.org/10.1016/j. neuron.2010.11.029 PMID: 21172607

8. Miller JeK, Ayzenshtat I, Carrillo-Reid L, Yuste R. Visual stimuli recruit intrinsically generated cortical ensembles. Proceedings of the National Academy of Sciences. 2014; 111(38):E4053–E4061. https:// doi.org/10.1073/pnas.1406077111

9. Seabrook TA, Huberman AD. Cortical Cliques: A Few Plastic Neurons Get All the Action. Neuron. 2015; 86(5):1113–1116. https://doi.org/10.1016/j.neuron.2015.05.039 PMID: 26050030

10. Barnes SJ, Sammons RP, Jacobsen RI, Mackie J, Keller GB, Keck T. Subnetwork-Specific Homeo- static Plasticity in Mouse Visual Cortex In Vivo. Neuron. 2015; 86(5):1290–1303. https://doi.org/10. 1016/j.neuron.2015.05.010 PMID: 26050045

11. Fuchs E, Ayali A, Ben-Jacob E, Boccaletti S. The formation of synchronization cliques during the devel- opment of modular neural networks. Physical biology. 2009; 6(3):036018. https://doi.org/10.1088/1478- 3975/6/3/036018 PMID: 19648664

12. Luccioli S, Ben-Jacob E, Barzilai A, Bonifazi P, Torcini A. Clique of functional hubs orchestrates popula- tion bursts in developmentally regulated neural networks. PLoS Comput Biol. 2014; 10(9):e1003823. https://doi.org/10.1371/journal.pcbi.1003823 PMID: 25255443
13. Reimann MW, Nolte M, Scolamiero M, Turner K, Perin R, Chindemi G, et al. Cliques of Neurons Bound into Cavities Provide a Missing Link between Structure and Function. Frontiers in Computational Neuroscience. 2017; 11:48. https://doi.org/10.3389/fncom.2017.00048 PMID: 28659782

14. Setareh H, Deger M, Petersen CC, Gerstner W. Cortical dynamics in presence of assemblies of densely connected weight-hub neurons. Frontiers in computational neuroscience. 2017;11. https://doi.org/10.3389/fncom.2017.00052 PMID: 28690508

15. Miles R, Wong RK. Single neurones can initiate synchronized population discharge in the hippocampus. Nature. 1983; https://doi.org/10.1038/306371a0 PMID: 6316152

16. Cy Li, Mm Poo, Dan Y. Burst spiking of a single cortical neuron modifies global brain state. Science. 2009; 324(5927):643-646. https://doi.org/10.1126/science.1169957

17. Bonifazi P, Goldin M, Picardo MA, Jorquera I, Cattani A, Bianconi G, et al. GABAergic hub neurons orchestrate synchrony in developing hippocampal networks. Science. 2009; 326(5958):1419-1424. https://doi.org/10.1126/science.1175509 PMID: 19965761

18. Brecht M, Schneider M, Sakmann B, Margrie TW. Whisker movements evoked by stimulation of single pyramidal cells in rat motor cortex. Nature. 2004; 427(6976):704-710. https://doi.org/10.1038/nature02266 PMID: 14973477

19. Cardin JA, Carlén M, Meletis K, Knoblich U, Zhang F, Deisseroth K, et al. Targeted optogenetic stimulation and recording of neurons in vivo using cell-type-specific expression of Channelrhodopsin-2. Nature protocols. 2010; 5(2):247. https://doi.org/10.1038/nprot.2009.228 PMID: 20134425

20. Ellender TJ, Nissen W, Colgin LL, Mann EO, Paulsen O. Priming of hippocampal population bursts by individual perisomatic-targeting interneurons. The Journal of neuroscience. 2010; 30(17):5979-5991. https://doi.org/10.1523/JNEUROSCI.3962-09.2010 PMID: 20427657

21. Houweling AR, Brecht M. Behavioural report of single neuron stimulation in somatosensory cortex. Nature. 2007; 451(7174):65-68. https://doi.org/10.1038/nature06447 PMID: 18094684

22. Picardo MA, Guigue P, Bonifazi P, Batista-Brito R, Allene C, Ribas A, et al. Pioneer GABA cells comprise a subpopulation of hub neurons in the developing hippocampus. Neuron. 2011; 71(4):695-709. https://doi.org/10.1016/j.neuron.2011.06.018 PMID: 21867885

23. Wester JC, McBain CJ. Interneurons differentially contribute to spontaneous network activity in the developing hippocampus dependent on their embryonic lineage. Journal of Neuroscience. 2016; 36(9):2646-2662. https://doi.org/10.1523/JNEUROSCI.4000-15.2016 PMID: 26937006

24. Mòdol L, Sousa VH, Malvache A, Tressard T, Baude A, Cossart R. Spatial Embryonic Origin Delineates GABAergic Hub Neurons Driving Network Dynamics in the Developing Entorhinal Cortex. Cerebral Cortex. 2017; 27(9):4649-4661. https://doi.org/10.1093/cercor/bhx198 PMID: 28922859

25. Carrillo-Reid L, Yang W, Bando Y, Peterka DS, Yuste R. Imprinting and recalling cortical ensembles. Science. 2016; 353(6300):691-694. https://doi.org/10.1126/science.aaf7560 PMID: 27516599

26. Cossart R. Operational hub cells: a morpho-physiologically diverse class of GABAergic neurons united by a common function. Current Opinion in Neurobiology. 2014; 26:51-56. https://doi.org/10.1016/j.conb.2013.12.002 PMID: 24650504

27. Ben-Ari Y. Developing networks play a similar melody. Trends in neurosciences. 2001; 24(6):353-360. https://doi.org/10.1016/S0166-2236(00)01813-0 PMID: 11356508

28. Yamada J, Okabe A, Toyoda H, Kilb W, Luhmann HJ, Fukuda A. Cl- uptake promoting depolarizing GABA actions in immature rat neocortical neurones is mediated by NKCC1. The Journal of physiology. 2004; 557(3):829-841. https://doi.org/10.1113/jphysiol.2004.062471 PMID: 15090604

29. Rheims S, Minlebaev M, Ivanov A, Represa A, Khazipov R, Holmes GL, et al. Excitatory GABA in rodent developing neocortex in vitro. Journal of neurophysiology. 2008; 100(2):609-619. https://doi.org/10.1152/jn.90402.2008 PMID: 18497364

30. Ben-Ari Y, Gaiarsa JL, Tyzio R, Khazipov R. GABA: a pioneer transmitter that excites immature neurons and generates primitive oscillations. Physiological Reviews. 2007; 87(4):1215-1284. https://doi.org/10.1152/physrev.00017.2006 PMID: 17928584

31. Donato F, Jacobsen RI, Moser MB, Moser EI. Stellate cells drive maturation of the entorhinal-hippocampal circuit. Science. 2017; 355(6330):eaai8178. https://doi.org/10.1126/science.aai8178 PMID: 28154241

32. Tien NW, Kerschensteiner D. Homeostatic plasticity in neural development. Neural development. 2018; 13(1):9. https://doi.org/10.1186/s13064-018-0105-x PMID: 29855353

33. Soloperto A, Bisio M, Palazzolo G, Chiappalone M, Bonifazi P, Difato F. Modulation of Neural Network Activity through Single Cell Ablation: An in Vitro Model of Minimally Invasive Neurosurgery. Molecules. 2016; 21(8). Available from: http://www.mdpi.com/1420-3049/21/8/1018. PMID: 27527143

34. Burkitt AN. A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input. Biological cybernetics. 2006; 95(1):1-19. https://doi.org/10.1007/s00422-006-0068-6 PMID: 16622699
35. Burkitt AN. A review of the integrate-and-fire neuron model: II. Inhomogeneous synaptic input and network properties. Biological cybernetics. 2006; 95(2):97-112. https://doi.org/10.1007/s00422-006-0082-8 PMID: 16821035

36. Tsodyks M, Uziel A, Markram H. Synchrony generation in recurrent networks with frequency-dependent synapses. The Journal of neuroscience: the official journal of the Society for Neuroscience. 2000 Jan; 20(1):50RC+. Available from: http://www.jneurosci.org/content/20/1/RC50 abstract. https://doi.org/10.1523/JNEUROSCI.20-01-j0003.2000

37. Fuchs EC, Neitz A, Pinna R, Melzer S, Caputi A, Monyer H. Local and distant input controlling excitation in layer II of the medial entorhinal cortex. Neuron. 2016; 89(1):194-208. https://doi.org/10.1016/j.neuron.2015.11.029 PMID: 26711115

38. Loebel A, Tsodyks M. Computation by ensemble synchronization in recurrent networks with synaptic depression. Journal of computational neuroscience. 2002; 13(2):111-124. https://doi.org/10.1023/A:1020110223441 PMID: 12215725

39. Stetter O, Battaglia D, Soriano J, Geisel T. Model-free reconstruction of excitatory neuronal connectivity from calcium imaging signals. PLoS computational biology. 2012; 8(8):e1002653. https://doi.org/10.1371/journal.pcbi.1002653 PMID: 22927808

40. Ge S, Goh EL, Sailor KA, Kitabatake Y, Ming GI, Song H. GABA regulates synaptic integration of newly generated neurons in the adult brain. Nature. 2005; 439(7076):589-593. https://doi.org/10.1038/nature04404 PMID: 16341203

41. Doetsch F, Hen R. Young and excitable: the function of new neurons in the adult mammalian brain. Current opinion in neurobiology. 2005; 15(1):121-128. https://doi.org/10.1016/j.conb.2005.01.018 PMID: 15721754

42. Karayannis T, García NVDM, Fishell GJ. Functional adaptation of cortical interneurons to attenuated activity is subtype-specific. Frontiers in neural circuits. 2012; 6. https://doi.org/10.3389/fncir.2012.00066 PMID: 23015781

43. Marissal T, Bonifazi P, Picardo MA, Nardou R, Petit LF, Baude A, et al. Pioneer glutamatergic cells develop into a morpho-functionally distinct population in the juvenile CA3 hippocampus. Nature communications. 2012; 3:1316. https://doi.org/10.1038/ncomms2318 PMID: 23271650

44. van den Heuvel MP, Sporns O. Network hubs in the human brain. Trends in cognitive sciences. 2013; 17(12):683-696. https://doi.org/10.1016/j.tics.2013.09.012 PMID: 24231140

45. Turrigiano G. Too many cooks? Intrinsic and synaptic homeostatic mechanisms in cortical circuit refinement. Annual review of neuroscience. 2011; 34:89-103. https://doi.org/10.1146/annurev-neuro-060909-153238 PMID: 21438687

46. Tanizawa T, Paul G, Cohen R, Havlin S, Stanley HE. Optimization of network robustness to waves of targeted and random attacks. Phys Rev E. 2005 Apr; 71:047101. Available from: https://link.aps.org/doi/10.1103/PhysRevE.71.047101.

47. Joyce KE, Hayasaka S, Laurienti PJ. The human functional brain network demonstrates structural and dynamical resilience to targeted attack. PLoS computational biology. 2013; 9(1):e1002885. https://doi.org/10.1371/journal.pcbi.1002885 PMID: 23358557

48. Ko H, Hofer SB, Pichler B, Buchanan KA, Sjöström PJ, Mrsic-Flogel TD. Functional specificity of local synaptic connections in neocortical networks. Nature. 2011; 473(7345):87. https://doi.org/10.1038/nature09880 PMID: 21478872

49. Carrillo-Reid L, Lopez-Huerta VG, Garcia-Munoz M, Theiss S, Arbuthnott GW. Cell assembly signatures defined by short-term synaptic plasticity in cortical networks. International journal of neural systems. 2015; 25(07):1550026. https://doi.org/10.1142/S0129065715500264 PMID: 26173906

50. Whittington MA, Traub R, Kopell N, Ermentrout B, Buhl E. Inhibition-based rhythms: experimental and mathematical observations on network dynamics. International journal of psychophysiology. 2000; 38(3):315-336. https://doi.org/10.1016/S0167-8760(00)00173-2 PMID: 11102670

51. Buzsáki G, Draguhn A. Neuronal oscillations in cortical networks. science. 2004; 304(5679): 1926-1929. https://doi.org/10.1126/science.1099745 PMID: 15218136

52. Snyder AC, Morais MJ, Smith MA. Dynamics of excitatory and inhibitory networks are differentially altered by selective attention. Journal of neurophysiology. 2016; 116(4):1807-1820. https://doi.org/10.1152/jn.00343.2016 PMID: 27466133

53. Bartolomei F, Khalil M, Wendling F, Sontheimer A, Régis J, Ranjeva JP, et al. Entorhinal cortex involvement in human mesial temporal lobe epilepsy: an electrophysiologic and volumetric study. Epilepsia. 2005; 46(5):677-687. https://doi.org/10.1111/j.1528-1167.2005.43804.x PMID: 15857433

54. Librizzi L, Losi G, Marcon I, Sessolo M, Scalmani P, Carmignoto G, et al. Interneuronal Network Activity at the Onset of Seizure-Like Events in Entorhinal Cortex Slices. Journal of Neuroscience. 2017; 37(43): 10398-10407. https://doi.org/10.1523/JNEUROSCI.3906-16.2017 PMID: 28947576
55. Tsien JZ. A Postulate on the Brain's Basic Wiring Logic. Trends in neurosciences. 2015; 38(11): 669-671. https://doi.org/10.1016/j.tins.2015.09.002 PMID: 26482260

56. Sorrentino F, di Bernardo M, Garofalo F, Chen G. Controllability of complex networks via pinning. Physical Review E. 2007; 75(4):046103. https://doi.org/10.1103/PhysRevE.75.046103

57. Fradkov AL, Miroshnik IV, Nikiforov VO. Nonlinear and adaptive control of complex systems. vol. 491. Springer Science & Business Media; 2013.

58. Liu YY, Slotine JJ, Barabási AL. Controllability of complex networks. Nature. 2011; 473(7346):167-173. https://doi.org/10.1038/nature10011 PMID: 21562557

59. Menichetti G, Dall'Asta L, Bianconi G. Network controllability is determined by the density of low in-degree and out-degree nodes. Physical review letters. 2014; 113(7):078701. https://doi.org/10.1103/PhysRevLett.113.078701 PMID: 25170736

60. Yuan Z, Zhao C, Di Z, Wang WX, Lai YC. Exact controllability of complex networks. Nature communications. 2013; 4. https://doi.org/10.1038/ncomms3447

61. Olmi S, Petkoski S, Guye M, Bartolomei F, Jirsa V. Controlling seizure propagation in large-scale brain networks. arXiv preprint arXiv:180403588. 2018.

62. Proix T, Jirsa VK, Bartolomei F, Guye M, Truccolo W. Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy. Nature communications. 2018; 9(1):1088. https://doi.org/10.1038/s41467-018-02973-y PMID: 29540685

63. Batista-Brito R, Close J, Machold R, Fishell G. The distinct temporal origins of olfactory bulb interneuron subtypes. The Journal of Neuroscience. 2008; 28(15):3966-3975. https://doi.org/10.1523/JNEUROSCI.5625-07.2008 PMID: 18400896

64. Miyoshi K, Miyoshi T, Siomi H. Many ways to generate microRNA-like small RNAs: non-canonical pathways for microRNA production. Molecular Genetics and Genomics. 2010; 284(2):95-103. https://doi.org/10.1007/s00438-010-0556-1 PMID: 20596726

65. Crépel V, Aronov D, Jorquera I, Represa A, Ben-Ari Y, Cossart R. A parturition-associated nonsynaptic coherent activity pattern in the developing hippocampus. Neuron. 2007; 54(1):105-120. https://doi.org/10.1016/j.neuron.2007.03.007 PMID: 17408581

66. Mongillo G, Barak O, Tsodyks M. Synaptic theory of working memory. Science. 2008; 319(5869): 1543-1546. https://doi.org/10.1126/science.1150769 PMID: 18339943

67. Allène C, Cattani A, Ackman JB, Bonifazi P, Aniksztejn L, Ben-Ari Y, et al. Sequential generation of two distinct synapse-driven network patterns in developing neocortex. The Journal of Neuroscience. 2008; 28(48):12851-12863. https://doi.org/10.1523/JNEUROSCI.3733-08.2008 PMID: 19036979

68. Tsodyks MV, Markram H. The neural code between neocortical pyramidal neurons depends on neurotransmitter release probability. Proceedings of the National Academy of Sciences. 1997; 94(2): 719-723. https://doi.org/10.1073/pnas.94.2.719

69. Turrigiano GG, Leslie KR, Desai NS, Rutherford LC, Nelson SB. Activity-dependent scaling of quantal amplitude in neocortical neurons. Nature. 1998; 391(6670):892-896. https://doi.org/10.1038/36103 PMID: 9495341

70. Turrigiano GG. The self-tuning neuron: synaptic scaling of excitatory synapses. Cell. 2008; 135(3): 422-435. https://doi.org/10.1016/j.cell.2008.10.008 PMID: 18984155

71. Zillmer R, Livi R, Politi A, Torcini A. Stability of the splay state in pulse-coupled networks. Phys Rev E. 2007 Oct; 76:046102. Available from: http://link.aps.org/doi/10.1103/PhysRevE.76.046102.

72. Gerstner W. Population dynamics of spiking neurons: fast transients, asynchronous states, and locking. Neural computation. 2000; 12(1):43-89. https://doi.org/10.1162/089976600300015899 PMID: 10636933

73. Bonifazi P, Difato F, Massobrio P, Breschi GL, Pasquale V, Levi T, et al. In vitro large-scale experimental and theoretical studies for the realization of bi-directional brain-prostheses. Frontiers in neural circuits. 2013; 7. https://doi.org/10.3389/fncir.2013.00040 PMID: 23503997