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):
Editor: Michele Giugliano,
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: 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 Economia, Industria y Competitividad of Spain (grant SAF2015-69484-R) (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.
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 . It is therefore essential 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 cliques , 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  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 circuitry. Neurons having such a network impact were recently termed operational hubs  or driver cells . 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 connectivity degree . This therefore represents a simple case where the circuit mechanism, promoting 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 developing Entorhinal Cortex (EC) , 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 population bursts we have upgraded and modified a network model based on excitatory leaky integrate-and-fire (LIF) neurons , previously developed to reproduce the functional properties of hub neurons in the developing hippocampal CA3 area . In such a model the population bursts (PBs), corresponding to GDPs in neonatal hippocampus , 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 hippocampus in . The model described in this paper contains two main differences with respect to the hippocampal model .
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 . 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  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 . 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 synchonization 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.
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 ) in the mice EC during developmental stage . Driver cells have been identified using calcium imaging experiments and they were characterized by the capability to impact network synchronization (namely, GDPs′ occurence) 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 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.
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 DO ≃ 7% are reported in (b-d.E), data for a neuron (from a different slice) with DO ≃ 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 νS = 0.33 Hz (ν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 Φ of the GDP versus the GDP index (denoting the successive GDP cycles). Φ 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 el1 which was silent in control condition and once stimulated with a current Istim = 15.135 mV was able to enhance population dynamics. Panels (e-g.S) refer to the excitatory driver LC cell el2 which was active in control condition with DO = 3% and once stimulated with a current Istim = 15.435 mV depressed the PB activity.
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 νS have been injected into the cell; (3) a final recovery period, where the cell is no more stimulated. The frequencies ν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 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 Φ measured with respect to the pre-stimulation phase (for the definition see Methods Eq (1) and in ). 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 D0 ≃ 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 νS, revealed opposite effects on the network behaviour. At lower stimulation frequency (ν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 (ν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 ν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).
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 νS = 0.33 Hz (νS = 1.0 Hz). Namely, the average IGI varied from 6.08 s in control conditions to 5.14 s during stimulation with νS = 0.33 Hz, and from 5.36 s to 6.68 s with νS = 1.0 Hz. Model. All the panels refer to the driver LC cell el3 connected in output to the el1 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 Istim = 15.42 mV and (d-f.S) to current Istim = 15.9 mV.
We also tested the possibility that single neuron stimulation could modify other features of network synchronization besides the GDP frequency, in particular we focused of 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 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 modifications 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) . 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 enthorinal 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  (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 asynchronous firing, thus resembling cortical GDPs’ occurrence in early stage networks. Similarly to the modeling reported in , we considered neuronal intrinsic excitabilities negatively correlated 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  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 Istim 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 DO = 0), that once stimulated could enhance of ≃ 30% the PB emission, thus leading to a decrease of the instantaneous phase Φ with respect to control condition. Panels (e-g.S) refer to a second neuron characterized by a low functional output connectivity, namely DO = 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 dynamics of each single neuron stimulation in relation to neuronal type, current injected and functional 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 stimulation 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 . In particular, the SND experiment consisted in recording the activity of the network in a fixed time interval Δt = 84 s when the considered 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 Istim 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 particular, we tested the response of the network to a broad range of stimulation amplitudes varying from 14.5 mV (slightly below the firing threshold for an isolated neuron Vth = 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 perturbation). As an example, for a stimulation current Istim = 15.90 mV the targeted neuron fired at a frequency ν ≃ 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 at least or more than 50% with respect to control conditions. In the specific, in analogy to what done in , 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 burst decrease at least or larger than 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) reports 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 ih1, eh1, eh2, eh3 was able to arrest completely the bursting dynamics within the considered time window, while in other two cases (for neurons ih2 and eh4) 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 non relevant modification in the number of emitted PBs, within the variability of the bursting activity in control conditions (Fig 3(a)).
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 ν 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 Istim = 15.90 mV for a time interval Δt = 84 s. The central horizontal dashed line shows the average number of PBs emitted in control conditions within an interval Δ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 DO of the neurons versus their total structural degree KT. A double square marks two neurons with overlapping properties. (g) Sketch of the structural connections among driver hub neurons and driver LC1 cells, the gray shaded 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.
The SNS confirmed that the neurons ih1, eh1, eh2, eh3 were capable to arrest the collective dynamics. Neurons eh4 and ih2 poorly impacted PB dynamics for the reported injected current, although for different values of Istim 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 , the SNS revealed also the presence of other 18 driver cells not identified by the SND capable to impact the occurence 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 neglibile 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 (eh1−4, ih1−2) were characterized by high values of functional out-degree, namely with an average functional degree DO = 68% ± 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 DO 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 KT, 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 DO and KT 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 KT = 15 ± 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 eh1, eh2, ih2 were supra-threshold, while neurons eh3, eh4, ih1 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 (eh1, ih1), (ih2, eh2) and (eh3, eh4), as reported in Table 1. The direct structural connections present among these couples (see also 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 (eh3, eh4), 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.
Table 1. Properties of driver hub cells in control condition.
For each driver hub cell (ih1, ih2, eh1, eh2, eh3, eh4) the columns report the intrinsic excitability (Ib), the average spiking frequency in control conditions (ν), the structural out-degree (KO) and in-degree (KI), the functional out-degree (DO) and in-degree (DI).
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 D0 = 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 ). In the following we will refer to them as el… or il1 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 drivers 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 charaterized by a higher structural connectivity degree KT with respect to driver hubs, namely KT = 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 (KT, D0). 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 (Ne = 90, Ni = 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, Ib was in the range [14.55: 15.42] mV (with average 15.0 ± 0.2 mV and ≃ 37% of the hubs supra-threshold), while the total connectivity KT was in the range [6: 31] (with average 16 ± 3 and a single hub with KT = 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 ± 0.1 mV and ≃ 99% of the LC cells below the firing threshold) and by a high KT in the range [14: 32] (with average 22 ± 4 and a single driver LC cell with KT = 14).
It is important to stress that in , 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 eh1 → ih1 → ih2 → eh2 → eh3 → eh4. 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 DI as reported in Table 1 indicating that they were among the first to fire during the PB build-up. In particular, eh1 had a functional in-degree DI 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  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 eh1 fired twice before the PBs (see Fig 3(e)), usually in-between the firing of eh2 and that of the pair (eh3, eh4), and this represents the main route, occurring for ≃ 85% of the PBs. Along the second route (present only for the ≃ 7% of the PBs), eh1 was firing the second time at the end of the sequence. The neuron eh1 fired essentially by following its natural period , 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 eh1 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 occurence 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, eh1 → ih1, ih2 → eh2 and eh2 → (eh3, eh4). 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 eh3 neuron represented the cornerstone of the clique, receiving the inhibitory input coming from the structural pair (eh1, ih1) and the excitatory one from the pair (ih2, eh2), with the activity of the neurons within each pair perfectly frequency locked. More specifically, eh1 entrained the activity of ih1 (below threshold in isolation) so that both neurons before a PB fired with a period quite similar to the natural period of eh1. The other pair (ih2, eh2) was controlled by the inhibitory action of ih2 that slowed down the activity of eh2, whose natural period was 60.6 ms, while before a PB ih2 and eh2 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 eh3 and eh4 and the coordinated activation of eh1 (and ih1) with the pair (ih2, eh2), 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 XOUT as defined in Eq (8). In particular, we considered the available resources only for the hub neurons eh3 and eh4 which were the last neurons of the clique to fire before the PB ignition. We have examined only these two hub neurons, because whenever eh3 and eh4 fired, a burst or an AB was always delivered.
Neurons eh3, eh4 were receiving high frequency excitatory inputs from eh2 (although the natural firing of eh2 was slowed down by the incoming inhibition of ih2) and high frequency inhibitory inputs from ih1 (entrained by the eh1, the neuron with highest firing frequency in the network). This competitive synaptic inputs resulted in a rare activation of eh3 compared to the higher frequency of excitatory inputs arriving from eh2. The period of occurrence of the ABs was comparable to the average interval between PBs (namely, TPB = 1.4±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 eh3 and eh4 was , 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 and should reach a sufficient high value in order to observe a PB, otherwise one had an AB. Furthermore, the value reached by and 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 just before the PB, at least for . At larger values the relationship was no more linear and a saturation was observable, due to the fact that could not overcome one.
We could conclude that the slow firing of the couple (eh3, eh4), moderate by the inhibitory action of ih1 on eh3, 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 and . This could be better understood by reconsidering the SND experiment on ih1, as expected the resection of neuron ih1 from the network led to a much higher activity of neurons eh3 and eh4, as shown in S8 Fig. However, this was not leading to the emission of any PBs, because in this case the value of and 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 . 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 ≤ Istim ≤ 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 silencing, of PBs, apart for two cells (namely, eh1 and eh4) 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 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 el1, el2 and el3, previously considered in subsection Numerical evidences of driver LC cells for comparison with experimentally identified LC cells, belonged to the class LC1 (see Figs 1 and 2), while we have no experimental examples of LC2 cells.
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 eh1, eh2, eh3, eh4, ih1, ih2 (a-f) as well as of the driver LC cells il1,el1,el2,el3, el4, el7 (g-l) for a wide range of the stimulation current Istim (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 Vth = 15 mV.
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 ≤ Istim ≤ 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 el1 (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 Istim of el1. Full circles, resp. open squares, signal the presence, resp. absence, of the corresponding neuron of the functional clique ih1, ih2, eh1, eh2, eh3, eh4, 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 ν of the LC cell (green line and right y-axis) during the current stimulation. The vertical (magenta) line marks the threshold value, Vth, while the vertical (blue) dashed line signals the intrinsic excitability of the LC cell in control condition.
For what concerns the driver hubs’ dynamics, PBs were generated in the network whenever the hubs eh2 and ih1, both structurally connected to eh3, were stimulated with currents smaller than the excitability 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 eh3 should not fire faster than the leader of the clique. If this was not the case, the inhibition (originating from ih1) would not be anymore sufficient to balance the excitation (carried by eh2) or viceversa, thus leading eh3 to operate outside the narrow current window where it should be located to promote collective activity (see Fig 5(c)). In the case of ih2 and eh4 the SNS produced a less pronounced 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 ih2. This is in agreement with what reported in Fig 3(a) for the SND, since the removal of neurons eh4 and ih2 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 il1, which was directly connected to the hub neuron eh3 (as shown in Fig 3(g)). A stimulation of il1 led to a decrease of the activity of eh3 and as a direct consequence of the PB activity. Fig 5(i) is devoted to el2, 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 el2, could be straightforwardly explained by the fact that el2 is directly connected to the inhibitory LC cell il1 and to the inhibitory driver hub ih1, thus performing an effective inhibitory action on the network, even if the stimulated driver el2 was excitatory. For the other six excitatory LC1 drivers acting on il1 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 eh3. In the two cases of driver LC1 cells able to block the population activity, but not acting through il1, these cells were exciting either eh1 or ih1, both belonging to the path with an inhibitory effect on eh3. The four remaining LC1 drivers that were able to reduce, but not to completely silence the population activity, acted either through eh4 (which was unable to block the PBs, even upon SND) or by simultaneously exciting and inhibiting eh3.
It should be remarked that nine of the previously discussed LC1 drivers could either enhance or reduce the PBs for different values of Istim. The double action of these neurons is exemplified by the two examples reported in Fig 5(h) and 5(j), which refer to neuron el1 and el3, 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 el1 was structurally connected to the inhibitory hub cell ih1, el1 was silent in control condition and once stimulated with a current Istim = 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 el1 with Istim ∈ [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 Istim) 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 el1 had no influence at a network level, and this is consistent with the response of ih1 upon SNS reported in Fig 5(e). However for higher stimulation current the clique of functional hubs is modified by the action of el1: 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 eh1 was no more relevant (in most cases) for the PB generation, and this could be explained by the fact that the inhibitory hub ih1 is now controlled directly by el1. 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 ih1. As a matter of fact, some lower activity in the network could be restored for even larger current values above Istim ≃ 15.65 mV. This effect can be explained in terms of the prevalence of facilitation over depression in the synapse connecting el1 with ih1, 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 el3 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: el3 was structurally connected, via the bridge neuron el1, to ih1, whose impact on the network was to arrest the bursting apart a very narrow range of stimulation currents (Fig 5(e)); el3 was also structurally connected to eh4, 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 Ib (below the firing threshold Vth) 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 el4 and el7. 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 el7 we measured <IGI > = 0.9 ± 0.6 s for Istim = 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 disappearence of ABs and as a consequence Inter GDP Interval become more regular, thus leading to an enhancement of the population activity. In particular, the disappearence of ABs was due to the fact that whenever eh3 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 eh3 were below the minimal value required in control conditions, as discussed in the subsection Synaptic resources and population bursts.
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 . 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 out-degrees has been revealed for the first time in . 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 mechanisms 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  with the corresponding disappearence 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 loose their role (see Fig 6(b) and 6(c)). At variance with purely excitatory networks , 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 . 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  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  where it has been shown that the emergence of sharp-wave in adult hippocampal slices was controlled by single perisomatic-targeting interneurons. 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 rhytms [50, 51] and for attentional modulation .
Entorhinal cortex is involved in human mesial temporal lobe epilepsy . A recent in vitro study on the onset of seizure-like events in EC slices  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 . 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 , 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 multiscale computational units whose hierarchical organization can lead to increasingly complex and specialized brain functions and can ground memory formations . 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 . 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 . Control methods based on linear stability analysis for seizure propagation on large-scale brain networks have been recently reported . 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 .
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 ν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.
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.
S3 Fig. Experiment—Robustness of the results with respect to ν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.
S4 Fig. Model—Setup for connectivity and excitability.
(a) Negative correlation between intrinsic excitability Ib and total connectivity KT. The (magenta) line indicates the threshold value Vth = 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.
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 ν under control conditions (shown in panel c)). During SNS experiments each neuron was stimulated with a DC step Istim mV for a time interval Δt = 84 s. The horizontal dashed line shows the average number of PBs emitted in control conditions within an interval Δ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.
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 Ib, (b) total structural connectivity KT. Dots (asteriskes) symbols refer to excitatory (inhibitory) neurons. The critical neurons ih1, ih2, eh1, eh2, eh3, eh4, 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, Vth = 15 mV, dividing supra-threshold from sub-threshold neurons. The neurons are ordered accordingly to their average firing rate in control conditions.
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 eh1 → ih1 → ih2 → eh2 → eh3 → eh4. During the out-of-burst activity, identical time lagged activation are preserved among the structurally connected pairs, namely eh1 → ih1, ih2 → eh2 and eh2 → (eh3, eh4).
S8 Fig. Model—Deletion (SND) of the inhibitory neuron ih1 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 ih1 (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 eh3, eh4. Bottom panel: average synaptic strength of the efferent connections of the two driver hubs neurons eh3, eh4; the output effective synaptic strength is measured by the average value of the fraction (black line), (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 , represented by the dashed lines (see also Fig 4 in the main text).
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 XOUT of the synaptic transmitters in the recovery phase associated to the efferent synapses. (g) Ratio NH/NL 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 el2 and el1, respectively, for the same stimulation protocol employed in Fig 1(b)–1(d.S) and 1(e)–1(g.S) in the main text.
S10 Fig. Model—Results of the SNS of LC drivers on the firing activity of the neurons of the clique.
Firing frequency ν of the driver hub neurons of the clique versus the current stimulation Istim during SNS of LC1 driver el1 (a) and LC2 driver el7 (b): (brown) stars, (red) crosses, (maroon) squares, (black) points, (green) diamonds, (blue) triangles refer respectively to eh1, eh2, eh3, eh4, ih1, ih2. The vertical (magenta) line marks the threshold value Vth, while the (black) arrows signal the firing frequency of the neurons of the clique in control condition .
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 Δ. Full circles, resp. open squares, signal the presence, resp. the absence, of the corresponding neuron of the functional clique (namely, ih1, ih2, eh1, eh2, eh3, eh4) 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 Δ and the number of LC drivers (identified in absence of noise) which survive for finite Δ (red triangles).
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 el1 on the network bursting mediated by short term synaptic plasticity between el1 and the driver hub cell ih1. a) Time average of the synaptic variables (black circles) and (red squares). b) Time average of the synaptic variable . c) Firing frequency of the driver hub neuron ih1. d) Firing frequency of the driver hub neuron eh3. e) Number of PBs emitted in the network. All the measures reported in the panels refer to a time window of Δt = 48 s.
S1 Table. Model—Routes leading to PBs.
Spike time delays ΔT between two successive firing of the neurons forming the functional clique along the main and secondary route leading to bursting. Neurons eh3 and eh4 are assumed to fire essentially at the same time, since eh4 fires almost immediately after eh3 within 0.03 − 0.04 ms. Notice that eh1 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 eh1 and this event is denoted by an asterisk in the table. The arrows indicate the order of firing in the sequence.