Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Transition between Functional Regimes in an Integrate-And-Fire Network Model of the Thalamus

Abstract

The thalamus is a key brain element in the processing of sensory information. During the sleep and awake states, this brain area is characterized by the presence of two distinct dynamical regimes: in the sleep state activity is dominated by spindle oscillations (7 − 15 Hz) weakly affected by external stimuli, while in the awake state the activity is primarily driven by external stimuli. Here we develop a simple and computationally efficient model of the thalamus that exhibits two dynamical regimes with different information-processing capabilities, and study the transition between them. The network model includes glutamatergic thalamocortical (TC) relay neurons and GABAergic reticular (RE) neurons described by adaptive integrate-and-fire models in which spikes are induced by either depolarization or hyperpolarization rebound. We found a range of connectivity conditions under which the thalamic network composed by these neurons displays the two aforementioned dynamical regimes. Our results show that TC-RE loops generate spindle-like oscillations and that a minimum level of clustering (i.e. local connectivity density) in the RE-RE connections is necessary for the coexistence of the two regimes. We also observe that the transition between the two regimes occurs when the external excitatory input on TC neurons (mimicking sensory stimulation) is large enough to cause a significant fraction of them to switch from hyperpolarization-rebound-driven firing to depolarization-driven firing. Overall, our model gives a novel and clear description of the role that the two types of neurons and their connectivity play in the dynamical regimes observed in the thalamus, and in the transition between them. These results pave the way for the development of efficient models of the transmission of sensory information from periphery to cortex.

Introduction

The thalamus is often identified as a relay station between subcortical and cortical areas, since all sensory pathways of the nervous system pass through it before reaching the cortex. Indeed, sensory inputs from visual, auditory and somato-sensory receptors reach the cortex through synapses on thalamocortical relay neurons in a specific region of the thalamus, which in turn projects into the corresponding area in the primary visual cortex. Along with these forward projections, there are local inhibitory neurons receiving inputs from feedback fibers from layer 6 to the corresponding thalamic nuclei [1]. It is thus reasonable to think that thalamus does not limit its activity to faithfully transmit information to the cortex, but it might play a role in gating and modulating the flow of information towards the cortex [24], i.e. in selecting which external information is supposed to reach the cortex and when. In particular, this view is coherent with the important role found to be played by the thalamus in sleep/arousal/wake process [57], and attention [810].

The main kind of excitatory neurons in the thalamus are the above-mentioned thalamocortical relay (TC) neurons. In vitro studies [11, 12] have revealed that these neurons can operate in different firing modes depending on their voltage level. Near the resting membrane potential, TC neurons can produce trains of spikes with frequency proportional to the amplitude of the injected current, due to voltage-dependent currents that generate action potentials [1]. This is usually called tonic mode. Alternatively, when TC neurons are hyperpolarized they can operate in a bursting mode, characterized by high-frequency bursts of action potentials (300 Hz) in response to hyperpolarization.

During slow-wave sleep, TC neurons display strong spindle oscillations (7 − 15 Hz) independently from external stimuli [1, 13]. In contrast, in the awake state TC neurons are known to vary their activity according to inputs coming from the associated receptor layers, and to affect in turn the activity of the associated primary sensory cortex. For instance, TC neurons belonging to the lateral geniculate nucleus (LGN) and the ventral posterior nucleus (VPN) are modulated by the retina [14] and by the tactile afferents [15], respectively, and modulate in turn the activity of primary visual and somatosensory cortical areas [4, 16, 17]. TC neurons are also key components of the above-mentioned gating role of the thalamus, contributing to the selection of salient information during selective attention [9].

As suggested by Crick in his seminal paper [2], the role of modulating the efficacy of sensory transmission of TC neurons is mainly played by the neurons of the reticular nucleus of the thalamus (RE neurons). In particular, the activation of RE neurons can strongly hyperpolarize TC neurons, which consequently undergo inhibitory rebound that gives rise to an endogenous oscillatory activity [18]. Specifically, spindles can be originated by TC bursts eliciting firing activity in RE cells. In turn RE bursts hyperpolarize TC cells, which consequently stop firing. When RE cells, lacking excitatory drive, stop firing too, the rebound of TC cells from hyperpolarization causes them to emit a burst of spikes and the cycle starts again. The overall process takes about 100 ms and generates rhythmic spindle oscillations. Therefore spindle generation is due to an interplay between TC and RE cells [19, 20]. Coherently with this fact, manipulating the activity of RE neurons was found to have behavioral consequences in attention tasks [9, 21].

During the awake state, TC cells undergo a transition and alternate this bursting mode with a tonic mode. As mentioned above, both modes are typical of TC neurons, and they could provide different frameworks for information processing, since during the bursting mode action potentials in the TC cell are not linked directly to EPSPs in that cell, whereas the opposite is true in the tonic mode. Therefore we expect that the bursting mode transmits information less efficiently than the tonic mode, in which an increase in the extra-thalamic inputs on TC neurons leads to a direct increase in the response of TC neurons [3]. How the thalamus exhibits the functional transition between the two regimes is not clear. In fact, a coherent view accounting for both TC and RE interactions and the resulting functional behavior of the thalamic network is still missing, due in particular to the relative paucity of simultaneous neurophysiological recordings of the two neuron types in vivo. In this context, the role of modeling becomes very relevant, due to its capacity to suggest candidate mechanisms for the generation of the observed behavior.

Modeling of thalamic networks has been investigated for more than 20 years [22, 23], during which network models have been developed that capture a wealth of thalamic phenomena [24]. However, almost all studies to date have adopted neuron models at least as complex as the Hodgkin and Huxley model [22, 25], probably due to the aforementioned role of rebound currents. We are aware of only one attempt to model realistically thalamic interactions with integrate-and-fire (IF) neurons [26]. Building upon that model, here we decided to study in a systematic manner the mechanism for spindle generation and the dynamics of the circuit in general. In particular, we focus on the ability of the thalamic network as a whole to switch between two dynamical regimes that display different external input sensitivity. We also study the role played in this phenomenon by the network architecture (connectivity and synaptic strength), ranging from loops of two neurons to the effect of sensory and cortical input on the whole thalamic network. To that end we have developed a thalamus TC-RE network model based on a particularly simple spiking neuron model, namely a different suited version of the adaptive exponential integrate-and-fire (aeIF) neuron model [27] for each neuron type. The choice of this neuron model allowed us to focus on a restricted number of parameters, specifically those related to physiological quantities influencing the rebound-driven oscillations and the tonic state [28]. Moreover, our aeIF thalamic model will be particularly suited to be interfaced with cortical column LIF models (e.g. [29]) to model accurately the whole thalamocortical loop.

In the Results section we build our network progressively. First, we show how our aeIF neurons reproduce the two activity modes of TC and RE neurons: the standard depolarizing regime [12] and the rebound from hyperpolarization [18]. Then we investigate how spindle oscillations are generated through TC-RE interaction as a function of their coupling and of the presence of external inputs, and how heterogeneity can be tamed by the interaction of different TC-RE loops. Analysis of the complete network leads to our two main results: (i) a critical value of RE-RE clustering (defined as the probability that two neurons connected to a third one are connected to each other) favors the presence of large scale spindle oscillations, and (ii) in the presence of clustering the network displays two dynamical regimes as the sensory input increases: it is insensitive to stimuli below a given intensity threshold, while above this threshold TC neurons (but not RE neurons) modulate their activity as a function of the input. Finally, we test that these conclusions hold also in the presence of cortical inputs impinging on the reticular neurons.

Results

The presence of two different dynamical regimes in the thalamus has been known for decades [3032]. This behavior can be linked to a specific property of the two main kinds of neurons in the thalamus, described above, glutamatergic thalamocortical relay (TC) neurons and GABAergic thalamic reticular (RE) neurons. Both types of neurons can fire either as a result of depolarizing driving or as a rebound due to hyperpolarizing driving. In the following we will show how we modified an existing (aeIF) model of thalamic neurons [26] to reproduce the two types of responses for both kinds of neurons (Fig 1). We investigated how and when the connectivity between the neurons displaying these properties induces a regime dominated by spindle oscillations, or responding to stimuli in a tonic-like mode. The analysis starts from two-neuron loops and extends up to full networks receiving input from the periphery and the cortical areas.

thumbnail
Fig 1. Dynamical properties of single RE and TC neurons as a function of input current.

(A) Depolarization induced activity of a RE neuron. Membrane voltage (top) and adaptation variable (middle) of a RE neuron in response to a depolarizing current (bottom). (B) Corresponding post-stimulus time histograms for increasing depolarizing currents. (C) Hyperpolarization-rebound activity of a RE neuron and (D) corresponding post-stimulus time histograms for increasing hyperpolarizing currents. Parameters a and b, representing respectively the dynamics and the strength of adaptation (see Eq (2)) of RE neurons are defined in this way: a = 0.4 μS and b = 0.02 nA. (E) Depolarization induced activity of a TC neuron and (F) corresponding post-stimulus time histograms for increasing depolarizing currents. (G) Hyperpolarization-rebound activity of a TC neuron and (H) corresponding post-stimulus time histograms for increasing hyperpolarizing currents. The values a and b are 0.2 μS and 0 nA. The current intensity in (A,C,E,G) is 1000 nA, while it varies between 1000 nA and 5000 nA in panels (B,D,F,H). VT = −50 mV is the threshold potential for both types of neurons. Other parameters are defined in the Materials and Methods section.

https://doi.org/10.1371/journal.pone.0161934.g001

Dynamics of single neurons

The first step towards reproducing the two dynamical regimes of the thalamus described above, and the transition between them, is to choose a single-neuron model able to capture the peculiar properties of thalamic neurons, and in particular the firing induced by hyperpolarization-driven rebound. To that end we selected a properly tuned adaptive exponential integrate-and-fire (aeIF) spiking neuron model [27, 33, 34] (see Materials and Methods Section) for each of the two thalamic neuron types considered. By tuning the key parameters of the aeIF model it is possible to adjust the dynamics and the strength of adaptation (parameters a and b in Eq (2), respectively, in the Materials and Methods Section) to reproduce the intrinsic dynamical modes typical of thalamic neurons.

For a = 0.4 μS and b = 0.02 nA, the RE aeIF neuron models (RE neuron from now on) exhibits regular firing activity in response to depolarizing stimuli (Fig 1A and 1B), while they display bursting activity in response to hyperpolarizing stimuli (Fig 1C and 1D), consistently with experimental findings [3537]. In particular, in response to a sustained depolarizing stimulus (Fig 1A), RE neurons display firing activity with a certain degree of spike-frequency adaptation that saturates before the end of the stimulus and stops neuronal firing. For sufficiently large applied currents, the response extends for the whole duration of the stimulus (Fig 1B). In response to a hyperpolarizing stimulus (Fig 1C and 1D), and due to the relatively large value of a, the neuron exhibits rebound bursting activity, also with spike-frequency adaptation, for the same spike threshold used in the depolarizing case.

TC neurons generally show a more robust bursting activity and a negligible level of spike-frequency adaptation [11] (see [1] for a review). This is achieved in the model by imposing a smaller value of a = 0.2 μS and b = 0 nA, thus making the adaptation strength negligible. In particular, in response to a depolarizing stimulus our TC neurons model produce patterns of firing activity (Fig 1E) with negligible spike-frequency adaptation (Fig 1F) (leading thereby to high firing activity for all the duration of the stimulus). In contrast, a hyperpolarizing stimulus leads to rebound bursting (Fig 1G) and moderate spike-frequency adaptation (larger than in RE neurons) (Fig 1H). In the case of depolarizing stimuli, characterized by negligible adaptation and regular firing activity, TC neurons exhibit an effective increase of activity (Fig 1F) according to the increasing external input and compatibly with the refractory period, where neuron is not allowed to fire. Therefore the firing activity increases proportionally with larger external sensory inputs. This is consistent with the strong input-output relation in the tonic mode (Fig 1E and 1F), in contrast with the bursting mode where there is no direct link between the EPSP and spike generation, which thus corresponds to a weak input-output correlation [38].

Overall these results show that the aeIF models properly capture the two firing modes (depolarization-driven and hyperpolarization-driven) for both TC and RE neurons. In the following we will show the transition between the two modes for TC neurons due to external inputs, and how recurrent activity drives a transition at the network level from stimulus-insensitive to stimulus-sensitive behavior.

Two-neuron loops

Before moving to large, structured networks we carefully analyzed the properties of the mutual interaction between TC and RE neurons, with two main objectives. First, we aimed at ensuring that our model is able to reproduce the observed global synchronization even in presence of biological variability leading for instance to heterogeneous intrinsic frequencies. Second, we intended to identify the more relevant connections in the generation of synchronized spindle patterns. Specifically, we studied different simple two-neuron loops formed by TC-RE and RE-RE neurons, and examined how self-sustained oscillatory patterns originated in these networks are modulated by synaptic strengths regulating the internal recurrent activity. We also studied the effect of GABA temporal decay dynamics on the frequency of oscillation, and the input-driven oscillatory pattern of a TC-RE loop. This analysis is informative towards the building of a full network.

We first built a minimal model of two bidirectionally coupled neurons, a RE neuron and a TC neuron (Fig 2A). Activating this RE-TC loop for 50 ms leads to oscillations that persist stably after the stimulus termination (Fig 2B). These oscillations are due to the rebound bursting properties of the TC relay cell, which is mutually connected with the RE neuron: the TC neuron provides depolarizing input to the RE neuron, which displays bursting activity (with two spikes per burst, see inset in Fig 2B) that generates strong hyperpolarization, followed by rebound firing activity in TC neurons. Consequently, in this configuration the RE neuron fires in response to depolarizing currents, while the TC neuron fires only in response to hyperpolarizing inputs. The time scale of the ISI interval (≈100 ms), which defines the spindle oscillation frequency, is determined by the TC-RE interplay, in particular by the time the TC neurons needs to recover from the RE-spikes induced hyperpolarization (Fig 2B).

thumbnail
Fig 2. Dynamical properties of two-neuron loops.

(A) Scheme of a two-neuron TC-RE loop. (B) Membrane voltage traces of the TC and RE neurons generated by this minimal TC-RE loop. (C) Interspike interval (ISI) distribution of the TC-RE loop as a function of the synaptic strength gTCRE. The value of gRETC is appropriately set to 550 μS in order to support self-sustained activity, while gTCRE varies between 10 μS and 60 μS. RE and TC ISI distributions are shown in the top and bottom plots, respectively. (D) ISI distribution of a TC-RE loop as a function of the synaptic strength gRETC. The value of gTCRE is chosen equal to 32 μS to reproduce the two-spike bursting dynamical regime of panel B while gRETC varies between 200 μS and 800 μS. RE and TC ISI distributions are shown in the top and bottom plots, respectively. (E) Scheme of a minimal purely reticular RE-RE loop. (F) ISI distribution of this loop as a function of the synaptic strength gRERE. gRERE varies between 200 μS and 800 μS. (G) Scheme of an input-driven two-neuron TC-RE loop. (H) ISI distribution of this loop as a function of external sensory input strength. RE and TC ISI distributions are shown in the top and bottom plots, respectively. The synaptic strengths are respectively: gRETC = 550 μS, gTCRE = 32 μS and gEXTTC = 1 μS.

https://doi.org/10.1371/journal.pone.0161934.g002

Next we investigated how these oscillatory patterns vary as a function of the synaptic strength gTCRE, keeping gRETC to a reference value of 550 μS. By increasing gTCRE, both the TC and RE neurons oscillate with higher frequencies, as can be seen from the decrease of the inter-spike interval (ISI) in Fig 2C (bottom). Stronger synaptic strengths enhance the firing activity of the RE neuron, which fires in advance along the oscillation cycle and thus leads the TC neuron to spike at an earlier phase. The net effect is an increase in the oscillation frequency. The RE neuron (Fig 2C, top) displays bursting activity in response to depolarizing input above a threshold value of gTCRE = 29 μS. It oscillates at around 11 Hz (inter-burst ISI ≈90 ms) with two spikes per burst with an intra-burst ISI ≈5 ms. By increasing the synaptic strength gTCRE, the neuron passes a second threshold gTCRE = 40 μS and presents three spikes per burst (three ISIs are present), eventually entering a regime in which the ISI approaches the intrinsic refractory period of the neuron (2.5 ms, see Material and Methods section).

Subsequently we performed the complementary analysis by fixing gTCRE to 32 μS (which led to two-spike bursting in the preceding analysis) and varying gRETC. Fig 2D shows that as gRETC is increased, the TC neuron oscillates with a gradually increasing frequency that stabilizes around 10.5 Hz (Fig 2D, bottom), while the RE neuron displays bursting activity with the same inter-burst ISI as the TC neuron and an intra-burst ISI of ≈3 ms (two-spikes-per-burst scenario of previous analysis) (Fig 2D. top). Note that the brief hyperpolarization induced in the TC cell by the firing of a single RE cell is able to trigger only one rebound spike, and consequently the number of spikes/burst in the RE cell remains constant. This is consistent with the results reported in Ref. [1], where spindle activity required at least a four-neuron network (see next Section).

Next we explored the dynamics of a purely GABAergic reticular RE-RE loop (Fig 2E) as a function of the synaptic strength gRERE. As Fig 2F shows, the RE neurons present a sustained strong and adapting bursting activity (corresponding to a wide range of intra-burst ISI) and for increasing values of the synaptic strength, the inter-burst ISI decreases. Importantly, unlike in previous studies, here the decreasing inter-burst ISI does not entail an increase in oscillation frequency, since here bursts last much longer (with more than 10 spikes per burst). This result shows that RE-RE synapses strengthen the rebound bursting properties and can be expected to enhance the bursting activity in a larger network.

In the simple TC-RE loop motif, the oscillation frequency can be tuned by the GABA decay time constant. For instance, by varying τdecay from 5 to 35 ms in the minimal model of Fig 2D, the frequency of the two neurons oscillates between ∼ 25 and 6 Hz (S1 Fig). This leads to corresponding changes in the ISI distributions (S2 and S3 Figs), without qualitative variations with respect to the behavior shown in Fig 2. The key role played by the GABA decay time constant in modulating the frequency of spindle oscillations is qualitatively similar to the way it affects gamma-range oscillation frequencies in LIF networks where rebound oscillations are not present [39].

After investigating the properties of stand-alone RE-TC loops, we moved to analyze an input-driven loop in which the TC neuron receives an external sensory input modeled as a Poisson pulse train with increasing rate (Fig 2G). We only considered inputs to TC, mimicking the sensory stimuli coming from the retina or the peripheral nervous system. We set reference values of gTCRE = 32 μS and gRETC = 550 μS and GABA τdecay = 20ms, for which the spontaneous activity (in the absence of external input) corresponds to low-frequency bursting with two spikes per burst. The value of GABA τdecay is lower in the full population model. When we increased the external input rate (Fig 2H) the ISI distribution was significantly different from the one observed in the absence of external stimulus (Fig 2D): both neurons show a strong variation in the bursting frequency due to the external stimulus, and the ISI displays a large variance due to the introduction of noise. In particular, the excitatory input on TC neurons decreases the ISI duration and determines the time scale of the oscillation cycle, by depolarizing the neuron and hence shortening the time needed for TC neurons to recover from the hyperpolarization due to inhibitory inputs (spindle frequency is now ≈15 Hz, while in Fig 2B is around ≈10 Hz). On the other hand, and consistently with Fig 2D, the RE neuron is in bursting mode for all values of external input, with the ISI approaching the refractory period.

Four-neuron motifs

As a last step before moving to the full network, we investigated several four-neuron motifs, made of two RE and two TC neurons, to understand what are the structural connectivity features more suitable to explain large oscillatory synchronization phenomena (spindle oscillations) during self-sustained activity in the bursting regime, even in presence of heterogeneity between neurons that leads to different oscillation frequencies. Previous work has shown [26] that aeIF models are able to reproduce this self-sustained oscillatory behavior in the form of periodic bursting, and that the minimal circuit reproducing the phenomenon is a circuit of two TC and two RE neurons fully connected with each other, with the exception of TC-TC connections, which are not present in the thalamus [40]. As in the case of the two-neuron loop, bursting is mainly due to the rebound bursting properties of TC cells and RE cells (Fig 3F) [1], and the oscillation frequency depends on the GABA temporal decay constant.

thumbnail
Fig 3. Four-neuron motifs in the form of coupled pairs of TC-RE loops.

The two TC-RE oscillators are bidirectionally coupled through (A) TC-RE connections, (B) RE-TC connections, (C) RE-RE connections, and (D) all three connections. (E) Frequency of the power spectral peak and (F) phase coherence at that frequency for the four different motifs. The power spectral density and phase coherence during self-sustained activity were averaged across 50 trials for random values of the GABA decay time (see text). GABA rise time and AMPA rise and decay times are set constant (see Materials and Methods section). When the corresponding connections exist in the motifs, the synaptic strengths are respectively: gRETC = 550 μS, gTCRE = 32 μS, and gRERE = 20 μS.

https://doi.org/10.1371/journal.pone.0161934.g003

Since oscillation frequency might slightly vary between loops, we checked the conditions for the onset of coherent oscillations. In particular, we studied different couplings between pairs of two-neuron TC-RE loops (which are equivalent to two bidirectionally coupled oscillators) with different intrinsic oscillation frequencies, and analyzed which coupling configuration leads more readily to oscillatory spindle patterns by examining the power spectrum of TC neurons and the phase coherence between them. Fig 3 shows the schemes of the different circuits explored depending on the coupling links being considered: TC-RE connections (Fig 3A), RE-TC connections (Fig 3B), RE-RE connections (Fig 3C) and all three types of connections (Fig 3D). For each circuit, we calculated the power spectral density and phase coherence between the two loops (see Materials and Methods section) by using the activity of TC neurons. The phase coherence is calculated by averaging 50 trials during self-sustained activity, each with a different GABA τdecay drawn from a Gaussian distribution with mean 20 ms and standard deviation 5 ms, which leads to variability in the frequencies of the two TC-RE loops being coupled.

Fig 3E shows the frequency at which the power spectrum of the TC neuron activity has its maximum, and Fig 3F the corresponding phase coherence at that frequency. The horizontal dashed red lines represent the corresponding values in the case of uncoupled loops. In the uncoupled case, the oscillation frequency is ≈10.4Hz and the loops are weakly synchronized (the phase coherence being ≈0.12). In configurations B and D, the two TC-RE loops strongly synchronize with a zero-lag phase (corresponding time lag is ≈0, not shown) with respect to the uncoupled case, while the loops are poorly zero-lag synchronized when only RE-RE connections are present. Therefore this result supports the idea that spindle generation is mainly due to an interplay between TC and RE cells [19, 20], which is enhanced by RE-RE connections.

Full thalamic network

We finally extended the size of the network to 500 neurons to capture the dynamics of a complex thalamic structure. Following experimental indications [4143], each RE projects on average to four TC or RE neurons, while TC neurons have on average one connection with RE neurons. Starting from these numbers we considered two network configurations, in order to investigate how the spindle oscillations are affected by network architecture during self-sustained activity. The first configuration was a purely random one (Fig 4A) with rewiring probability RP = 1 (see Methods), while the second favored RE-RE clustering with rewiring probability RP = 0.25 (Fig 4B). This amount of clustering can be expected from networks of spatially distributed neurons such as the one being considered here, since neighboring neurons should be more likely connected to each other than those that are further apart [4447].

thumbnail
Fig 4. Spindle self-sustained activity generated by a full network of TC-RE neurons depending on RE-RE clustering.

(A) Connectivity matrix of a random TC-RE network. The presynaptic neurons are represented in the x axis and the postsynaptic neurons in the y axis. The network is made of 500 neurons, of which the first 250 are RE neurons and the remaining ones are TC neurons. (B) Connectivity matrix in the presence of RE-RE clustering (rewiring probability RP = 0.25) (C) Membrane voltage dynamics of a couple of arbitrarily chosen TC and RE neurons in the case of random network. (D) Membrane voltage dynamics of a couple of arbitrarily chosen TC and RE neurons in the presence of clustering: evidence of typical spindle oscillations. (E,F) Distribution of inter-spike intervals (along the horizontal axis, color-coded and normalized to unit area) as a function of the rewiring probability (along the vertical axis) for RE (panel E) and TC (panel F) neurons. The synaptic strengths are respectively: gRETC = 300 μS, gTCRE = 200 μS and gRERE = 300 μS.

https://doi.org/10.1371/journal.pone.0161934.g004

We found that in the random network (Fig 4A), temporally irregular bursting is dominant (Fig 4C). On the other hand, in the presence of RE-RE clustering (Fig 4B) the network shows quite regular and synchronized self-sustained spindle oscillations at 8 Hz (Fig 4D). In order to characterize and quantify the bursting regular state (or spindle rhythm) and distinguish it from irregular bursting, we studied the inter-burst interval distribution (in particular the probability of a peak of ISI distribution at the already identified spindle timescale of 50 ms) as a function of the rewiring probability RP of the architecture (see Methods for details). Our results, shown in Fig 4E, reveal that fully regular networks (RP = 0, each neuron projects regularly to a fixed number of adjacent neurons) cannot support regular bursting activity and are often almost silent (with a firing rate of around 0.4 spikes/s, results not shown). At the other extreme, fully random networks (RP = 1) show sustained activity with temporally irregular bursting of TC and RE neurons. Between these two conditions, there is an optimal rewiring probability (RP ∼ 0.25) showing a relatively large ISI peak corresponding to frequency ∼8.5 Hz. The fraction of neurons displaying a large inter-burst ISI peak decreases substantially for increasing rewiring probability, namely when going towards fully random networks. Intuitively, given that connections between thalamic circuits are local but sparse [4143], excitatory synapses are very sparse and they are more effective when they impinge on small clusters of RE-RE neurons, enhancing and modulating the oscillatory spindle rhythm.

Given the results obtained above, we decided to study a network with the critical degree of clustering (obtained for RP = 0.25), and shift our focus from self sustained activity to constant external sensory input of different intensities impinging on TC neurons. We tested if by increasing the external input on these neurons the network showed a transition from bursting to tonic mode, which could be associated with the switch from sleep to awake state [3032]. The relation between input and output in the bursting mode is weak [38], but we expected to see that TC neurons (projecting to the cortex) increase their firing rate with the input strength when the network goes from bursting to tonic mode. Fig 5A shows the firing rate of TC (red) and RE (blue) neurons for increasing external sensory input on TC neurons. The case of input S = 0 spikes/s corresponds to the self-sustained condition discussed above. By increasing the input strength, the network displays a transition in the firing rate of TC neurons at around S = 50 spikes/s, after which the response of the thalamus increases sub-linearly with the external input. We interpret this as an indication of the switch from a purely bursting mode to a temporally irregular state. Note that the driver of this transition is the response of the recurrent activity to the external sensory input, since we did not change the intrinsic parameters of the model.

thumbnail
Fig 5. Bursting and tonic modes displayed by a TC-RE network with RE-RE clustering as a function of external input on TC neurons.

(A) Firing rate of TC (red) and RE (blue) neurons as a function of external driving input impinging on TC neurons. (B,C) ISI distribution as a function of external driving input on TC neurons of RE (B) and TC (C) neurons. (D) Mutual information between the set of increasing external stimulus (0-150 spikes/s) and the neural response given by the firing rate of TC and RE neurons. Different external sensory inputs are considered for the two regimes, following panel A: 0-50 spikes/s for the bursting mode and 60-150 spikes/s for the tonic mode. The white dashed line in the bar plots refers to significance threshold (p < 0.05, bootstrap test). The measures are averaged over 100 trials for each external stimulus. (E,F) Adaptation variable w of RE (E) and TC (F) neurons (color coded) as a function of the external input on TC neurons, averaged across 100 trials for each external stimulus. (G) Number of positive w values (depolarizing events, green) and negative w values (rebound events, black) of TC neurons. (H) Coefficient of variation of the ISI for both TC and RE cells as the input rate on TC neurons increases. The synaptic strengths are respectively: gRETC = 300 μS, gTCRE = 200 μS and gRERE = 300 μS.

https://doi.org/10.1371/journal.pone.0161934.g005

In order to explore this scenario further, we calculated the ISI distribution of RE and TC neurons by averaging over 100 trials for each different stimulus S. The RE neurons are the most insensitive to increasing external input, as can be seen in Fig 5B. On the other hand the fraction of TC neurons displaying a large inter-burst ISI decreased as the stimulus intensity surpasses a critical value (going from region S1 to region S2 in Fig 5A), and a corresponding increase of the intra-burst ISI peak approaching the refractory period (2.5 ms, see Materials and Methods section). We classified this as a further signature of a transition between a bursting mode and an irregular firing regime.

Next we calculated the information about the stimuli carried by the firing rates of the TC and RE neurons in the two different regimes. To that end we used the mutual information (see Methods), which quantifies the reduction of the uncertainty in predicting the applied stimulus given a single observation of the triggered response. In this case we considered a rate code, i.e. we selected as response the average firing rate over the whole stimulation [48]. Fig 5C compares I(S1; FR) and I(S2; FR) between the firing rates of TC (red) and RE (blue) neurons and the set of stimuli S1 and S2, where S1 ranges between 0 and 50 spikes/s, while S2 varies from 60 to 150 spikes/s, corresponding to the two dynamical regimes of Fig 5A. The figure clearly shows that in the bursting mode both the RE and TC neurons carry a lower information (0.13 bit, p < 0.05: bootstrap test), in comparison with the information encoded by TC neurons during the tonic mode (≈0.7 bit, p < 0.05: bootstrap test). RE neurons during the tonic mode do not encode significant information, in fact their firing rate decreases with respect to the bursting regime and after that remains constant for all inputs. These results show that the information about the stimulus that the thalamus carries (and is then potentially able to convey to the cortex) is much higher in the tonic mode, since in that regime spontaneous activity is enhanced and this contributes to keeping a strong relation between input and output and thus to minimizing rectification of the response [38].

In order to further interpret this transition, we examined the nature of each TC and RE spike by checking the sign of the adaptation variable w at the spiking time of each neuron. A positive value of w indicates that neuron fires via a depolarizing input (see Fig 1), while if negative we classify it is as a rebound spike. Fig 5E shows that RE neurons spike mostly due to a rebound in response to hyperpolarizing inputs (coming only from internal RE-RE clustered connections) for all the range of sensory input over TC neurons. TC neurons, in turn, also fire mainly in response to incoming hyperpolarizing currents (in this case coming from RE neurons) during the burst mode (Fig 5F), and after the transition from bursting to tonic mode a fraction of the spikes occur in response to depolarizing external inputs. Thus the transition occurring at around S = 50 spikes/s, shown in Fig 5A, underlies a shift in the spiking mechanism profile. This is confirmed in Fig 5G, which shows a quantitative estimation of the effective number of excitatory-driven spikes (blue) and inhibitory-rebound spikes (red) as the external input increases. Note that even though we divided the input rate in two regimes for the sake of simplicity, the lack of discontinuity in the increase of depolarization-driven spikes in this panel hints toward a smooth transition between the two firing modes. Finally we quantify the degree of temporal regularity of these two dynamic regimes by calculating the coefficient of variation (CV) of the interspike intervals (ISI) as a function of the input rate on TC neurons, consistently with the other analysis. The CV is calculated as the ratio between the standard deviation and the mean of the ISIs, averaged over all cells the network. In the S1 range the CV is approximately 0.3, indicating a very regular firing typical of bursting modes. In the S2 range, the CV increases with the input rate up to a value 0.6, showing an increasing irregularity of the firing activity, due to the coexistence of the tonic and bursting modes.

So far we have considered a thalamic network receiving an external sensory input impinging on TC neurons. We finally included a corticothalamic input [49] projecting to RE neurons. Fig 6A shows that the transition dynamics is not altered by the addition of a constant input from the cortex, which results only on an increase of the firing rate for both kind of neurons. The onset and the increase of depolarization spikes occur for similar levels of inputs (Fig 6B). The amount of information carried by RE and TC neurons in the two different regimes is relatively unaltered (Fig 6C and 6D), supporting the hypothesis that the information carried by projecting neurons during the tonic mode is higher than in the bursting mode. Interestingly, by increasing the amplitude of the cortical input on RE (from 1000 to 2000 spikes/s), the information encoded by TC neurons is increased for the tonic mode (from 0.6 to 0.66 bit, p < 0.05, bootstrap test) (Fig 6D). This result highlights the role of the intrinsic rebound bursting properties of TC neurons, which are essential in the generation of the spindle rhythm. They could also reinforce the role of corticothalamic feedback in information processing, for instance by recruiting TC neurons through inhibition and thus modulating TC firing rate [49]. To show the relevance of the rebound bursting properties of TC neurons, we plotted in S4 Fig the firing rate of TC and RE neurons, the ISI distribution and the w distribution at a fixed rate of external sensory input on TC neurons (150 spikes/s), for different levels of cortical input.

thumbnail
Fig 6. Bursting and tonic modes displayed by the TC-RE network with RE-RE clustering as a function of external input on TC neurons for different corticothalamic inputs.

(A) Firing rate of TC (red) and RE (blue) neurons as a function of the external driving input impinging on TC neurons for different corticothalamic input amplitudes. (B) Number of positive (depolarizing, red) and negative (rebound, blue) w values of TC spikes for different corticothalamic inputs. The w values are averaged across 100 trials for each external stimulus. (C,D) Mutual information carried by the firing rate of TC (red) and RE (blue) neurons with a cortico-thalamic input of (C) 1000 spikes/s and (D) 2000 spikes/s, calculated between the set of increasing sensory stimuli (10 − 150 spikes/s) and the neural response given by the firing rate. The white dashed lines in the bars refer to the significance threshold (p < 0.05, bootstrap test). Measures are averaged over 100 trials for each external stimulus. The synaptic strengths are respectively: gRETC = 300 μS, gTCRE = 200 μS and gRERE = 300 μS.

https://doi.org/10.1371/journal.pone.0161934.g006

Discussion

We have presented an adaptive exponential integrate-and-fire (aeIF) network model that is able to reproduce spindle oscillations and the transition between a stimulus-insensitive and a stimulus-sensitive state of the thalamus. Coherently to what was shown experimentally through direct optogenetic stimulation [18], in our model spindle oscillations are generated by RE activation leading to TC bursts as rebound from inhibition. Our simulations suggest that (i) these oscillations are stable for a specific range of RE-RE connection clustering, (ii) for external stimuli below a given threshold the network is in a purely rebound-bursting state insensitive to external stimuli, while when this threshold is crossed there is a non-zero contribution of the spikes due to depolarization, and this makes the TC neurons (and not the RE neurons) of the network sensitive to the stimulus intensity coherently with experimental observation.

Advantages and limitations of aeIF models

Choosing a simple model for the single neurons allowed us to focus on capturing the network effects. This choice also opens a number of interesting perspectives: due to their relative simplicity, IF models can be tackled analytically [28, 50, 51], and have substantially expanded the search for basic canonical computations [52]. Finally, most primary sensory cortex network models are built on IF neurons [29, 53], and hence aeIF neurons seem a more coherent choice to build models of corticothalamic interactions [54]. In our model, the switch from inhibitory-rebound-driven activity to depolarization-driven firing is proposed to represent a switch from sleep to awake state [7]. The information analysis shown in Figs 5 and 6 reveals the separation between a stimulus-independent state (sleep) and a stimulus-sensitive state (wakefulness). We did not directly deal with the role of thalamus, and in particular RE neuronal activity, in attention, to which a wealth of works have been devoted [9, 10] after the seminal intuition of Crick [2]. To compare our model results with these experimental observations we should (1) contrast different states inside the awake regime, and (2) take into account the temporal structure of the TC spike trains rather than their rate alone. This is certainly feasible on the ground of the results presented here, but is beyond the scope of this paper. We emphasize that our model is based on single-neuron models that are much simpler than those used previously. Although this has a number of advantages as discussed above, some features of thalamic behavior that are captured by more detailed models are not reproduced by our model. For instance, our spindle oscillations constitute a stable state, both in small and large TC-RE networks, and do not reproduce the wax-and-wane dynamics that has been observed experimentally [55], and which has been reproduced by more detailed models that take explicitly into account the dynamics of hyperpolarization-activated cation currents [56].

TC-RE loop studies

A recent computational paper [25] investigated the role of TC-RE interactions from a perspective complementary to the one discussed in this paper, using a Hodgkin-Huxley model much more detailed than the aeIF adopted here, and limiting the investigation only to minimal loops such as those we described here (Fig 2). Notwithstanding the higher realism of their model, the functional properties at the single-neuron level are similar to those described here (compare the two Figs 1 of the two works). Moreover, Willis and colleagues highlighted the fact that open-loops between TC and RE neurons might play a functional role in the thalamus, and indeed in our full network (Fig 4 and following) both open and closed TC-RE loops are taken into account. In a recent paper [21], Brown and collaborators stimulated optogenetically RE neurons, simultaneously recording from them. They found that the majority of those neurons (10/17) decreased significantly their firing rate, and only a minority of them (4/17) displayed a significant increase. At the same time they found that the activity of the TC neurons was inhibited, with functional consequences on the cortex. The interpretation of the authors was that a small increase in RE activity was sufficient to inhibit TC activity. Our model offers a simpler explanation: since most TC neurons fire due to hyperpolarization rebound, a decrease in RE activity can be associated to a decrease in TC firing (see Fig 6). Indeed, stimulating RE neurons has been shown to alter the temporal structure of TC neuron firing, without changing their average firing rate [18].

Perspectives

The present work focuses on the thalamus, taking into account only stable external inputs from the periphery to TC neurons or from the cortex to RE neurons. Preliminary analysis suggested that an accurate description of thalamocortical inputs and corticothalamic feedbacks required a separate study. In the future this network will be integrated in a full corticothalamic model comprising a primary visual cortex network (inspired by previous works [57, 58]). The following step will be to take into account (i) the layered structure of the cortex [29] and (ii) areas of the thalamus and the cortex associated to different sensory receptive fields and their interactions. Another interesting continuation of this work would be to contribute to the open challenge of modeling the Local Field Potential of the thalamus [59]. We recently showed [60] that an integrate-and-fire model like the one presented here can be combined with morphological data and transmembrane current simulation [61] to capture the LFP dynamics in a patch of cortex. Since morphological data are available for the thalamus, a similar procedure can be applied to the network introduced here, and would hopefully shed light on the way extracellular signals and neural activity are linked in this area, thus enhancing the possibility of experimental validations of the thalamic models. The potential applications of this work include the study of the consequences of deep brain stimulation (DBS). Thalamic DBS has been shown to contribute to the symptom mitigation of a variety of neural diseases including Parkinson [62] and Tourette’s syndrome [63]. However, the precise mechanisms of this mitigation are not completely clear, nor is the procedure to design specific trains of stimulations suited for different patients/conditions. Neural models are already exploited to test DBS patterns [64]. We think that a simple yet efficient model like the one presented here can valuably contribute to this field.

Materials and Methods

Computational model

We have used the adaptive exponential integrate-and-fire (aeIF) model [27], which is an evolution of a two-variable integrate-and-fire (IF) model proposed by Izhikevich [33], enriched by an exponential nonlinearity around the spike threshold, as in the exponential IF model of Fourcaud-Trocme et al. [34]. The combination of these two models leads to the aeIF formulated by Brette and Gerstner [27]. A detailed analysis of the dynamics of this system can be found in [28].

Single neuron model.

According to the aeIF model, the equations describing the evolution of membrane voltage of neurons in the thalamus are: (1) (2) The first equation describes the evolution of the membrane voltage: the capacitive current through the membrane with capacitance Cm = 1 nF equals the ionic currents, the adaptation current w and the input current I. The ionic currents are the ohmic leak current defined by the resting leak conductance gL = 0.05 μS and the resting voltage potential EL = −60 mV, and the exponential term which reproduces the Na+– current that is responsible for the generation of spikes. With this term we assume that the activation of Na+–channels is instantaneous (thus neglecting their activation), with Δ denoting the steepness of the exponential approach to threshold, taken equal to Δ = 2.5 mV, and VT = −50 mV is the threshold potential. The membrane time constant is τm = Cm/gL. When V is pushed over the threshold, the exponential term provides a positive feedback and a spike is emitted (occurring ideally at the time when V diverges towards infinity). The membrane potential is reset to Vr = −60 mV, and the adaptation variable w is increased by a value b (ww + b). After the spike, the neuron cannot spike again during a refractory period (2.5 ms).

The second equation describes the dynamics of w, with time constant τw = 600 ms. The parameter a (in μS) quantifies a conductance that mediates subthreshold adaptation, while the increment b (in nA) at each spike takes into account spike-triggering adaptation (it regulates the strength of adaptation). When the input current I to the neuron at rest reaches a critical value (1000 nA) the resting state is destabilized, leading to repetitive spiking for large regions of parameter space [45]. Without adaptation (a = b = 0) the model produces tonic spiking. Neurons in general can show a reduction in the firing frequency of their spike response if they are stimulated with a square pulse or step, known as spike frequency adaptation (SFA). With this model, an increase of a or b leads to SFA, characterized by a gradual increase in the inter-spike interval (ISI) until a steady-state spike frequency is reached. For b = 0, the model generates responses with a negative level of adaptation similar to the fast-spiking (FS) cells encountered in the cortex, often classified as inhibitory neurons. The strength of adaptation can be modulated by varying the parameter b, to get weakly adapting cells [45]. In what follows we choose values of a and b such that the range of values of w is on the order of those used in the literature [65].

In order to reproduce the peculiar properties of TC and RE when operating in bursting mode, we adopted specific values of a and b. When they are in their tonic mode, TC and RE neurons behave similarly to excitatory regular spiking (RS) neurons and inhibitory fast spiking (FS) neurons found in the cortex. With a = 0.4 μS, b = 0.02 μA, neurons display bursting activity in response to both depolarizing and hyperpolarizing stimuli typical of RE neurons. In contrast, with a = 0.2 μS, b = 0.0 μA, neurons display responses with moderate adaptation and strong rebound bursts, like TC neurons. RE and TC neurons can display different regimes (beyond bursting and tonic, fast spiking (FS), regular spiking (RS)) by tuning the parameters a and b [33, 45, 65]. In both cases, since a/gL > τm/τw we are in the parameter regime in which rebound firing is possible, as demonstrated by Touboul and Brette [28].

Thalamic network model.

The network is made of TC and RE cells, endowed with intrinsic properties and topographic connectivity specific to the thalamus [26]. Here we considered a network of 500 neurons, half of which are TC neurons and the other half being RE neurons. Given that thalamic interneurons do not contribute to the development of internal dynamics such as oscillations, they are neglected. Axonal projections within the thalamic circuitry are local but sparse. The excitatory projections from TC to RE had a connection probability of 1%, while RE to TC inhibitory projections were more dense, with a connection probability of 4%. The same density was assumed from inhibitory connections between RE cells. The structural connectivity is built starting from a ring network and then randomly rewiring with probability RP. This process allows us to control the clustering coefficient, which quantifies the connectedness or local connectivity of the network (in terms of the probability that two nodes that are connected to a common node are also connected between them). A rewiring probability equal to 0 implies a regular network with large clustering (provided the coupling extends beyond nearest neighbours), whereas a rewiring probability equal to 1 implies a completely random network with small clustering.

In Fig 4 we introduced different degrees of clustering (by tuning the rewiring probability RP), and according to the results we adopted RP = 0.25 for the continuation of the analysis. The network model was constructed based on this aeIF model, according to the following equations [26]: (3) (4) where Vi is the membrane potential of neuron i, and all parameters are as in Eqs (1) and (2), but were indexed to allow variations according to the cell type. The term ∑j gij(ViEj) accounts for the synaptic current coming from the neighboring neurons impinging on a neuronal cell, where gij is the conductance of the synapse from neuron j to neuron i (which can be zero), and Ej is the reversal potential of the synapse (Ej = 0 mV for excitatory synapses and −80 mV for inhibitory synapses). Synaptic conductances are described by: (5) where τdecay and τrise are the decay and rise synaptic time, respectively, and is constant and depends on the type of synapses and network (see Table 1). Once the presynaptic cell fires, gij exponentially increases up to a certain value, after which gij decays exponentially with a fixed time constant (5 ms for excitation and 10 ms for inhibition). Synaptic delays are equal to 1 ms. Different synaptic strengths are considered (see Table 2), depending on the network type. If different values are used, they are indicated in the captions of each figure. Note that the conductance values used here are higher than the ones observed experimentally. This is done to compensate for the unrealistically low amount of incoming inputs, due to the fact that we are considering small networks. This synaptic strength rescaling is a common practice in computational neuroscience [26, 66, 67]. In fact in order to design configurations with plausible synaptic conductance values it would be necessary to build unwieldy networks (with 10,000 to over 100,000 neurons), with large numbers of synapses per neuron (>500).

thumbnail
Table 1. Values of temporal rise and decay constants for RE and TC. These values are standard choices in the computational neuroscience literature [6870].

https://doi.org/10.1371/journal.pone.0161934.t001

thumbnail
Table 2. Values of synaptic strengths for a network of 500 neurons.

https://doi.org/10.1371/journal.pone.0161934.t002

To initiate activity, during the first 50 ms a number of randomly-chosen neurons were stimulated by an incoming current (with synaptic strength g = 40 μS), representing an heterogeneous Poisson train of excitatory presynaptic potential with an instantaneous event rate λ(t) that varies following an Ornstein-Uhlenbeck process: (6) where σ(t) is the standard deviation of the noise and is set to 0.6 spikes/s. τ is set to 16 ms, leading to a power spectrum for the λ time series that is approximately flat up to a cut-off frequency . η(t) is a Gaussian white noise of mean zero and intensity unity. In simulations in which we did not take into account external input after 50 ms, no input was given to the network, and thus the activity states described here are self-sustained with no external input or added noise. The only source of noise was the random connectivity. In simulations in which we took into account external sensory inputs, after 5 s of self-sustained activity we injected for 10 s homogeneous Poisson processes with rate comprised between 10 and 150 spikes/s.

Spectral analysis

We computed the power spectral density of LFPs using the Welch method: the signal is split up into 32768 point segments with 50% overlap. The overlapping segments are windowed with a Hamming window. The modified periodogram is calculated by computing the discrete Fourier Transform, and then calculating the square magnitude of the result. The modified periodograms are then averaged to obtain the PSD estimate, which reduces the variance of the individual power measurements. Spectral quantities and phase coherence are averaged over 50 trials.

Phase coherence

Phase coherence is calculated as in [71]: (7) where x and y denote the two signals, and Sxy(f, n) is the cross-spectrum between them. Since in each trial the cross spectral density is normalized by its amplitude, each term of the sum is a unit-length vector representation of the phase relation Δϕ(f, n). In other words, Δϕ(f, n) = ϕyϕx is the phase lag between the two signals at frequency f in the data segment n. Hence Cxy(f) quantifies how broad is the distribution of Δϕ(f, n) within the 2π-cycle. Averaging Δϕ(f, n) across all N data segments provides a mean angle Δϕ(f).

Mutual information

We calculate the mutual information I(S; FR) between the set of stimuli S given by the external Poisson inputs with different rates described above and the response FR, firing rate, as follows. Given that we are interested in how the specific neurons encode and carry information, in this case we select as response the average firing rate FR over the whole stimulation (other responses such as the power spectrum can be considered [72]). We consider as stimuli different inputs with increasing amplitude (from 0 to 150 spikes/s) impinging on TC neurons. We compute the information between the stimulus S and the response firing rate as: (8) where P(s) is the probability of having a stimulus s (equal to the inverse of the total number of different external firing rates, which act as stimuli), P(r) is the probability of observing a firing rate r across all trials in response to any stimulus, and P(r|s) is the probability of observing a firing rate r in response to a single stimulus s. I(S; FR) quantifies the reduction of uncertainty about the stimulus that can be gained from observing a single-trial neural response, measured in units of bits (1 bit means a reduction of uncertainty of a factor of two) [48]. This measure allows us to evaluate how well the firing rate r of both type of neurons encodes the stimulus s.

An important issue to be solved regarding the calculation of the mutual information is that it requires knowledge of the full stimulus-response probability distributions, and obviously these probabilities are calculated from a finite number of stimulus-response trials. This leads to the so-called limited sampling bias, which constitutes a systematic error in the estimate of information. We used the method described in [73] to estimate the bias of the information quantity and then we checked for the residual bias by applying a bootstrap procedure, in which mutual information is calculated when the stimuli and responses are paired at random. If the information quantity is not zero (as it should be in the case of non-finite samples), this is an indication of the bias, and the bootstrap estimate of this error should be removed from the mutual information. After applying these procedures, the information quantity estimation could be defined as significant. Several toolboxes provide different bias-correction techniques, which allow accurate estimates of information theoretic quantities from realistically collectable amounts of data [74, 75]. In order to accomplish those tasks, we used the Information Breakdown Toolbox (ibTB), a MATLAB toolbox implementing several information estimates and bias corrections [75].

Supporting Information

S1 Fig. Effect of the GABA decay time on the two-neuron TC-RE loop.

Interspike Interval (ISI) distribution of a TC-RE loop as a function of the GABA decay time for RE (A) and TC (B) neurons. The synaptic strengths are gRETC = 550 μS, gTCRE = 32 μS.

https://doi.org/10.1371/journal.pone.0161934.s001

(EPS)

S2 Fig. Effect of the synaptic strength on the two-neuron TC-RE loop.

Interspike Interval (ISI) distribution of a TC-RE loop as a function of the synaptic strength for RE (A,C) and TC (B,D) neurons. As in Fig 2C, in A and B the value of gRETC is appropriately set to 550 μS in order to support self-sustained activity, while gTCRE varies between 10 μS and 60 μS. In C and D, the value of gTCRE is chosen equal to 40 μS to reproduce the two-spike bursting dynamical regime, while gRETC varies between 200 μS and 800 μS. GABA decay τdecay is set equal to 10 ms in the two cases.

https://doi.org/10.1371/journal.pone.0161934.s002

(EPS)

S3 Fig. Effect of the synaptic strengths on the two-neuron RE-RE motif.

Interspike interval (ISI) distribution of a minimal purely reticular RE-RE motif as a function of the synaptic strength gRERE for three different values of the GABA decay time: (A) 5 ms, (B) 10 ms, (C) 20 ms. gRERE varies between 300 μS and 800 μS.

https://doi.org/10.1371/journal.pone.0161934.s003

(EPS)

S4 Fig. Influence of corticothalamic input on a full TC-RE network.

(A) Firing rate, (B) ISI distribution, and (C) distribution of the adaptation variable w of RE and TC neurons as a function of corticothalamic input. The external sensory input it set to 150 spikes/s. The synaptic strengths are respectively: gRETC = 300 μS, gTCRE = 200 μS and gRERE = 300 μS. Bars colors in panels (A) and (B) coincide with the lines colors in the other panels.

https://doi.org/10.1371/journal.pone.0161934.s004

(EPS)

Acknowledgments

We thank two anonymous reviewers for their constructive comments. This work was supported by the European Commission under ITN project NETT (FP7 contract 289146), by the Spanish Ministry of Economy and Competitiveness and FEDER (project FIS2015-66503-C3-1-P), and by the Generalitat de Catalunya (project 2014SGR0947). A.M. was supported by the NEBIAS European project (EUFP7-ICT-611687), by the PRIN/ HandBot Italian project (CUP: B81J12002680008; prot.: 20102YF2RY), and by the Italian Ministry of Foreign Affairs and International Cooperation, Directorate General for Country Promotion (Economy, Culture and Science) Unit for Scientific and Technological Cooperation, via the Italy-Sweden bilateral research project on “Brain network mechanisms for integration of natural tactile input patterns”. J.G.O. also acknowledges support from the ICREA Academia programme and from the “María de Maeztu” Programme for Units of Excellence in R&D (Spanish Ministry of Economy and Competitiveness, MDM-2014-0370).

Author Contributions

  1. Conceived and designed the experiments: AB JGO AM.
  2. Performed the experiments: AB.
  3. Analyzed the data: AB.
  4. Contributed reagents/materials/analysis tools: AB.
  5. Wrote the paper: AB JGO AM.

References

  1. 1. Destexhe a, Sejnowski TJ. Interactions between membrane conductances underlying thalamocortical slow-wave oscillations. Physiological reviews. 2003 Oct;83(4):1401–53. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2927823&tool=pmcentrez&rendertype=abstract. pmid:14506309
  2. 2. Crick F. Function of the thalamic reticular complex: the searchlight hypothesis. Proceedings of the National Academy of Sciences of the United States of America. 1984 Jul;81(14):4586–4590. Available from: http://www.pnas.org/content/81/14/4586.abstract. pmid:6589612
  3. 3. Sherman SM, Guillery RW. The role of the thalamus in the flow of information to the cortex. Philosophical transactions of the Royal Society of London Series B, Biological sciences. 2002;357(November):1695–1708. pmid:12626004
  4. 4. Reinhold K, Lien AD, Scanziani M. Distinct recurrent versus afferent dynamics in cortical visual processing. Nat Neurosci. 2015 Dec;18(12):1789–97. pmid:26502263
  5. 5. Llinás RR, Paré D. Of dreaming and wakefulness. Neuroscience. 1991;44(3):521–535. Available from: http://view.ncbi.nlm.nih.gov/pubmed/1754050. pmid:1754050
  6. 6. Steriade M, McCormick DA, Sejnowski TJ. Thalamocortical Oscillations in the Sleeping and Aroused Brain. Science. 1993;262:679–85. pmid:8235588
  7. 7. Dang-Vu TTT, Schabus M, Desseilles M, Albouy G, Boly M, Darsaud A, et al. Spontaneous neural activity during human slow wave sleep. Proceedings of the National Academy of Sciences of the United States of America. 2008 Sep;105(39):15160–15165. Available from: http://dx.doi.org/10.1073/pnas.0801819105. pmid:18815373
  8. 8. Guillery RW, Feig SL, Lozsadi DA. Paying attention to the thalamic reticular nucleus. Trends in Neurosciences. 1998;21:28–32.
  9. 9. McAlonan K, Cavanaugh J, Wurtz RH. Guarding the gateway to cortex with attention in visual thalamus. Nature. 2008 Nov;456(7220):391–394. Available from: http://dx.doi.org/10.1038/nature07382. pmid:18849967
  10. 10. Wimmer K, Compte A, Roxin A, Peixoto D, Renart A, de la Rocha J. Sensory integration dynamics in a hierarchical network explains choice probabilities in cortical area MT. Nature communications. 2015;6. Available from: http://view.ncbi.nlm.nih.gov/pubmed/25649611.
  11. 11. Llinás R, Jahnsen H. Electrophysiology of mammalian thalamic neurones in vitro. Nature. 1982 Jun;297(5865):406–408. Available from: http://www.nature.com/nature/journal/v297/n5865/pdf/297406a0.pdf. pmid:7078650
  12. 12. Jahnsen H, Llinás R. Electrophysiological properties of guinea-pig thalamic neurones: an in vitro study. The Journal of physiology. 1984 Apr;349:205–26. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1199334&tool=pmcentrez&rendertype=abstract. pmid:6737292
  13. 13. Andersen P, Eccles J. Inhibitory Phasing of Neuronal Discharge. Nature. 1962 Nov;196(4855):645–647. Available from: http://www.nature.com/nature/journal/v196/n4855/pdf/196645a0.pdf. pmid:14012799
  14. 14. Butts DA, Weng C, Jin J, Yeh CI, Lesica NA, Alonso JM, et al. Temporal precision in the neural code and the timescales of natural vision. Nature. 2007 Sep;449(7158):92–95. Available from: http://dx.doi.org/10.1038/nature06105. pmid:17805296
  15. 15. Land PW, Blas ALd, Reddy N. Immunocytochemical localization of GABAA receptors in rat somatosensory cortex and effects of tactile deprivation. Somatosensory & motor research. 1995;12(2):127–141.
  16. 16. Gilbert CD, Wiesel TN. Receptive field dynamics in adult primary visual cortex. Nature. 1992;356(6365):150–152. pmid:1545866
  17. 17. Pais-Vieira M, Kunicki C, Tseng PH, Martin J, Lebedev M, Nicolelis MA. Cortical and thalamic contributions to response dynamics across layers of the primary somatosensory cortex during tactile discrimination. Journal of neurophysiology. 2015;114(3):1652–1676. pmid:26180115
  18. 18. Halassa MM, Chen Z, Wimmer RD, Brunetti PM, Zhao S, Zikopoulos B, et al. State-dependent architecture of thalamic reticular subnetworks. Cell. 2014;158(4):808–821. pmid:25126786
  19. 19. McCormick DA. Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. Progress in neurobiology. 1992;39(4):337–388. pmid:1354387
  20. 20. Destexhe A, McCormick DA, Sejnowski TJ. A model for 8–10 Hz spindling in interconnected thalamic relay and reticularis neurons. Biophysical journal. 1993 Dec;65(6):2473–7. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1225988&tool=pmcentrez&rendertype=abstract. pmid:8312485
  21. 21. Lewis LD, Voigts J, Flores FJ, Schmitt LI, Wilson MA, Halassa MM, et al. Thalamic reticular nucleus induces fast and local modulation of arousal state. eLife. 2015;p. e08760. pmid:26460547
  22. 22. Destexhe A, Contreras D, Sejnowski TJ, Steriade M. A model of spindle rhythmicity in the isolated thalamic reticular nucleus. Journal of neurophysiology. 1994;72(2):803–818. pmid:7527077
  23. 23. Golomb D, Wang XJ, Rinzel J. Propagation of spindle waves in a thalamic slice model. Journal of Neurophysiology. 1996;75(2):750–769. pmid:8714650
  24. 24. Muller L, Destexhe A. Propagating waves in thalamus, cortex and the thalamocortical system: experiments and models. Journal of Physiology-Paris. 2012;106(5):222–238.
  25. 25. Willis AM, Slater BJ, Gribkova ED, Llano DA. Open-loop organization of thalamic reticular nucleus and dorsal thalamus: A computational model. Journal of neurophysiology. 2015;114(4):2353–2367. pmid:26289472
  26. 26. Destexhe A. Self-sustained asynchronous irregular states and Up-Down states in thalamic, cortical and thalamocortical networks of nonlinear integrate-and-fire neurons. Journal of computational neuroscience. 2009 Dec;27(3):493–506. Available from: http://www.ncbi.nlm.nih.gov/pubmed/19499317. pmid:19499317
  27. 27. Brette R, Gerstner W. Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity. Journal of Neurophysiology. 2005;94(5):3637–3642. pmid:16014787
  28. 28. Touboul J, Brette R. Dynamics and bifurcations of the adaptive exponential integrate-and-fire model. Biological cybernetics. 2008;99(4-5):319–334. pmid:19011921
  29. 29. Potjans TC, Diesmann M. The cell-type specific cortical microcircuit: relating structure and activity in a full-scale spiking network model. Cerebral Cortex. 2014;24(3):785–806. pmid:23203991
  30. 30. Steriade M, McCormick DA, Sejnowski TJ. Thalamocortical oscillations in the sleeping and aroused brain. Science (New York, NY). 1993 Oct;262(5134):679–685. Available from: http://dx.doi.org/10.1126/science.8235588.
  31. 31. Steriade M, McCarley RW. Brain Control of Wakefulness and Sleep. Boston, MA: Springer US; 2005. Available from: http://dx.doi.org/10.1007/b102230.
  32. 32. Livingstone MS, Hubel DH. Effects of sleep and arousal on the processing of visual information in the cat. Nature. 1981 Jun;291(5816):554–561. Available from: http://view.ncbi.nlm.nih.gov/pubmed/6165893. pmid:6165893
  33. 33. Izhikevich EM. Which model to use for cortical spiking neurons? Neural Networks, IEEE Transactions on. 2004 Sept;15(5):1063–1070.
  34. 34. Fourcaud-Trocme’ N, Hansel D, van Vreeswijk C, Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating inputs. The Journal of Neuroscience. 2003;23(5):11628–40.
  35. 35. Contreras D, Curro-Dossi R, Steriade M. Electrophysiological properties of cat reticular thalamic neurones in vivo. Journal of Physiology (Lond). 1993;470:273–94.
  36. 36. Domich L, Oakson G, Steriade M. Thalamic burst patterns in the naturally sleeping cat: a comparison between cortically projecting and reticularis neurones. J Neurophysiol. 1986;379:429–449.
  37. 37. Steriade M. Neuronal Substrates of Sleep and Epilepsy. Cambridge University Press; 2003. Available from: http://www.amazon.com/exec/obidos/redirect?tag=citeulike07-20&path=ASIN/0521817072.
  38. 38. Murray Sherman S. Tonic and burst firing: Dual modes of thalamocortical relay. Trends in Neurosciences. 2001;24(2):122–126.
  39. 39. Brunel N, Wang XJ. What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J Neurophysiol. 2003 Jul;90(1):415–30. Available from: http://www.ncbi.nlm.nih.gov/pubmed/12611969. pmid:12611969
  40. 40. Izhikevich EM, Edelman GM. Large-scale model of mammalian thalamocortical systems. Proceedings of the National Academy of Sciences. 2008 Mar;105(9):3593–3598. Available from: http://dx.doi.org/10.1073/pnas.0712231105.
  41. 41. FitzGibbon T, Tevah LV, Jervie-Sefton A. Connections between the reticular nucleus of the thalamus and pulvinar-lateralis posterior complex: A WGA-HRP study. Journal of Comparative Neurology. 1995;363:489–504. pmid:8847413
  42. 42. Minderhoud JM. An anatomical study of the efferent connections of the thalamic reticular nucleus. Experimental Brain Research. 1971;112:435–446. pmid:5579571
  43. 43. Kim U, Sanchez-Vives M, McCormick D. Functional dynamics of GABAergic inhibition in the thalamus. Science (New York, NY). 1997 October;278(5335):130–134. Available from: http://dx.doi.org/10.1126/science.278.5335.130.
  44. 44. Litwin-Kumar A, Doiron B. Slow dynamics and high variability in balanced cortical networks with clustered connections. Nature neuroscience. 2012;15(11):1498–1505. pmid:23001062
  45. 45. Ladenbauer J, Augustin M, Shiau L, Obermayer K. Impact of adaptation currents on synchronization of coupled exponential integrate-and-fire neurons. PLoS computational biology. 2012 Jan;8(4):e1002478. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3325187&tool=pmcentrez&rendertype=abstract. pmid:22511861
  46. 46. Litwin-Kumar A, Doiron B. Formation and maintenance of neuronal assemblies through synaptic plasticity. Nature communications. 2014;5. pmid:25395015
  47. 47. Klinshov VV, Teramae Jn, Nekorkin VI, Fukai T. Dense neuron clustering explains connectivity statistics in cortical microcircuits. PloS one. 2014;9(4):e94292. pmid:24732632
  48. 48. Ince R, Mazzoni A, Petersen RS, Panzeri S. Open source tools for the information theoretic analysis of neural data. Front Neurosci. 2010 Jan;4(1):62–70. Available from: http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=2891486&tool=pmcentrez&rendertype=abstract. pmid:20730105
  49. 49. Destexhe A. Modelling corticothalamic feedback and the gating of the thalamus by the cerebral cortex. Journal of Physiology Paris. 2000;94:391–410.
  50. 50. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. Journal of computational neuroscience. 2000;8(3):183–208. pmid:10809012
  51. 51. Barbieri F, Mazzoni A, Logothetis NK, Panzeri S, Brunel N. Stimulus dependence of local field potential spectra: experiment versus theory. The Journal of Neuroscience. 2014;34(44):14589–14605. pmid:25355213
  52. 52. Schwartz A. The promise of neurotechnology. Science. 2015;350(6256):11–11. pmid:26430090
  53. 53. Battaglia D, Hansel D. Synchronous chaos and broad band gamma rhythm in a minimal multi-layer model of primary visual cortex. PLoS computational biology. 2011 Oct;7(10):e1002176+. Available from: http://dx.doi.org/10.1371/journal.pcbi.1002176. pmid:21998568
  54. 54. Muller L, Destexhe A. Propagating waves in thalamus, cortex and the thalamocortical system: experiments and models. Journal of Physiology-Paris. 2012;106(5):222–238.
  55. 55. Bal T, McCormick DA. What stops synchronized thalamocortical oscillations? Neuron. 1996;17(2):297–308. pmid:8780653
  56. 56. Bazhenov M, Timofeev I, Steriade M, Sejnowski T. Self-sustained rhythmic activity in the thalamic reticular nucleus mediated by depolarizing GABAA receptor potentials. Nature neuroscience. 1999;2(2):168–174. pmid:10195202
  57. 57. Mazzoni A, Panzeri S, Logothetis NK, Brunel N. Encoding of naturalistic stimuli by local field potential spectra in networks of excitatory and inhibitory neurons. PLoS Comput Biol. 2008;4(12)(12):e1000239. pmid:19079571
  58. 58. Mazzoni A, Whittingstall K, Brunel N, Logothetis NK, Panzeri S. Understanding the relationships between spike rate and delta/gamma frequency bands of LFPs and EEGs using a local cortical network model. Neuroimage. 2010;52(3):956–972. pmid:20026218
  59. 59. Einevoll GT, Kayser C, Logothetis NK, Panzeri S. Modelling and analysis of local field potentials for studying the function of cortical circuits. Nature Reviews Neuroscience. 2013;14(11):770–785. pmid:24135696
  60. 60. Mazzoni A, H L, H C, A L, S P, T EG. Computing the local field potential (LFP) from integrate-and-fire models. PLoS Comput Biol. 2015;. pmid:26657024
  61. 61. Lindén H, Hagen E, Leski S, Norheim ES, Pettersen KH, Einevoll GT. LFPy: a tool for biophysical simulation of extracellular potentials generated by detailed model neurons. Frontiers in neuroinformatics. 2013;7.
  62. 62. Tasker RR. Deep brain stimulation is preferable to thalamotomy for tremor suppression. Surgical neurology. 1998;49(2):145–153. pmid:9457264
  63. 63. Servello D, Porta M, Sassi M, Brambilla A, Robertson MM. Deep brain stimulation in 18 patients with severe Gilles de la Tourette syndrome refractory to treatment: the surgery and stimulation. Journal of Neurology, Neurosurgery & Psychiatry. 2008;79(2):136–142.
  64. 64. McIntyre CC, Grill WM, Sherman DL, Thakor NV. Cellular effects of deep brain stimulation: model-based analysis of activation and inhibition. Journal of neurophysiology. 2004;91(4):1457–1469. pmid:14668299
  65. 65. Naud R, Marcille N, Clopath C, Gerstner W. Firing patterns in the adaptive exponential integrate-and-fire model. Biological Cybernetics. 2008;99(4-5):335–347. Available from: http://dx.doi.org/10.1007/s00422-008-0264-7. pmid:19011922
  66. 66. Kumar A, Schrader S, Aertsen A, Rotter S. The high-conductance state of cortical networks. Neural computation. 2008;20(1):1–43. pmid:18044999
  67. 67. El Boustani S, Pospischil M, Rudolph-Lilith M, Destexhe A. Activated cortical states: experiments, analyses and models. Journal of Physiology-Paris. 2007;101(1):99–109.
  68. 68. Destexhe A, Mainen ZF, Sejnowski TJ. Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. Journal of computational neuroscience. 1994;1(3):195–230. pmid:8792231
  69. 69. Destexhe A, Paré D. Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo. Journal of neurophysiology. 1999;81(4):1531–1547.
  70. 70. McCormick DA, Wang Z, Huguenard J. Neurotransmitter control of neocortical neuronal activity and excitability. Cerebral Cortex. 1993;3(5):387–398. pmid:7903176
  71. 71. Womelsdorf T, Schoffelen JM, Oostenveld R, Singer W, Desimone R, Engel AK, et al. Modulation of neuronal interactions through neuronal synchronization. Science. 2007 Jun;316(5831):1609–12. pmid:17569862
  72. 72. Barardi A, Sancristóbal B, Garcia-Ojalvo J. Phase-Coherence Transitions and Communication in the Gamma Range between Delay-Coupled Neuronal Populations. PLoS computational biology. 2014 Jul;10(7). Available from: http://view.ncbi.nlm.nih.gov/pubmed/25058021.
  73. 73. Panzeri S, Treves A. Analytical estimates of limited sampling biases in different information measures. Network. 1996;7(7):87–107.
  74. 74. Victor JD. Approaches to information-theoretic analysis of neural activity. Biol Theory. 2006;1:302–316.
  75. 75. Magri C, Whittingstall K, Singh V, Logothetis N, Panzeri S. A toolbox for the fast information analysis of multiple-site LFP, EEG and spike train recordings. BMC Neurosci. 2009;10(1):81. pmid:19607698