Causal Role of Thalamic Interneurons in Brain State Transitions: A Study Using a Neural Mass Model Implementing Synaptic Kinetics

1. Introduction

The thalamic interneurons (IN) are believed to play a fundamental role in linking retinal sensory input to visual perception by modulating thalamocortical alpha rhythms (8–13 Hz) recorded through electroencephalogram (EEG) from the occipital cortex (Lörincz et al., 2009). Alpha rhythms are traditionally known to represent an idling state of the brain when a subject is awake but relaxed with eyes closed. Subsequently, it emerged that alpha rhythms play a key role in controlling perception and are of significance during an awake cognitive state. A transition of EEG from the alpha to the theta (4–7 Hz) band in an awake state is associated with several neurological disorders such as neurogenic pain, tinnitus, Parkinson's disease, and is termed as thalamocortical dysrhythmia (TCD) (Sarnthein et al., 2003; Hughes and Crunelli, 2005; Llinas et al., 2005); a similar symptom in Alzheimer's disease is termed as “slowing” (a decrease of dominant frequency) of the alpha rhythms. On the other hand, alpha to theta band transition in a quiet resting state is a marker of a change of brain state from wakefulness to drowsiness. In the visual pathway, local field potential (LFP) recordings of alpha and theta rhythms from the thalamic lateral geniculate nucleus (LGN) show a high correlation with EEG recorded simultaneously from the occipital cortex (Lopes da Silva et al., 1973; Hughes and Crunelli, 2006). In this regard, synchronous oscillatory patterns showing waxing-and-waning of amplitude within the alpha band is a well known hallmark of EEG and LFP in an awake and “resting” state (i.e., devoid of sensory input or mental task), and have been a matter of extensive research; their generation is attributed to the feed-forward and feed-back connectivity between cell populations of the Thalamocortical Relay (TCR: the main carriers of sensory information to the cortex) and the Thalamic Reticular Nucleus (TRN: a thin sheet of inhibitory cells surrounding the thalamus, receiving “copies” of communications between the TCR and visual cortex) (Steriade et al., 1990). Thus, computational models simulating thalamocortical dynamics have focused on the TRN as the primary inhibitory influence on the TCR, and thereby, on the cortex (Destexhe et al., 1996; Golomb et al., 1996; Stam et al., 1999; Robinson et al., 2004; Grimbert and Faugeras, 2006; Bhattacharya et al., 2011a; Wang et al., 2014); the inhibitory influence of IN is largely ignored. A similar gap is seen in experimental research investigating the functional impact of the IN cells on the thalamocortical oscillations (Crunelli et al., 2006; Halassa et al., 2014), exceptions being some early research in Crunelli et al. (1988) and Zhu et al. (1999a). This is in spite of the IN constituting around 20–25% of the total number of cells in almost all thalamic nuclei processing sensory information in mammals; around 47% of the synaptic afferents of the the IN are from the information carrying spiking neurons of the retina (Sherman, 2004; Jones, 2007). Moreover, the critical role of the IN in the visual signal processing by the LGN and information transmission in the retino-geniculo-cortical pathway is now well established (Dublin and Cleland, 1977; Wang et al., 2007; Babadi et al., 2010; Saalmann and Kastner, 2011; Wang et al., 2011a,b; Pressler and Regehr, 2013; Bastos et al., 2014; Hirsch et al., 2015); also, their physiology and spiking characteristics are now understood fairly well (Pape and McCormick, 1995; Zhu et al., 1999a,b; Cox et al., 2003). Thus, it is surprising that the importance of the causality of IN on brain rhythms is underestimated in experimental research, perhaps due to the lack of appropriate technology (Zhu et al., 1999a,b) that prevented proper recordings of the IN population dynamics. The emphasis on the role of IN in brain rhythms was revived only recently when Lörincz et al. (2008), while studying the waking state alpha rhythm and their response to cortico-thalamic inputs, report the distinct inhibitory effect of IN over TCR; this validates the findings in early research on IN cell dynamics (Crunelli et al., 1988) where the authors report the significant GABA-ergic influence of IN on the TCR oscillatory dynamics in the rat LGN. In addition, Lörincz et al. (2008) report a minimal role of the TRN on the TCR dynamics in the awake state. This is in agreement with a recent study using a computational model on the role of thalamic cells in the disappearance of alpha rhythms in sleep-wake transitions (Bond et al., 2014), where we have reported a dramatic effect of the presence of the IN cells on the TRN response; the research presented here builds on this prior work.

Neural mass computational models (NMM) (Marreiros et al., 2008; Moran et al., 2013) are often used to simulate brain rhythms recorded in LFP and EEG that are believed to be generated through dynamic interaction between networks of meso-scale (10⁴−10⁷ neurons) neuronal populations. These models were conceptualized in the works of Wilson and Cowan (1973) and Freeman (1975), and later popularized by the classic work on alpha rhythm by Lopes da Silva et al. (1974) and Zetterberg et al. (1978). Subsequently, this classic model was extended in Jansen and Rit (1995) and Suffczyński (2000), and used extensively in model-based research of neuro-psychiatric disorders (Wendling et al., 2002; David and Friston, 2003; Suffczyński et al., 2004; Modolo et al., 2013; Wang et al., 2013; Taylor et al., 2014). In previous works, we have proposed an enhancement to state-of-the-art neural mass models by replacing the “alpha function” (Rall, 1967) with kinetic models of Glutamatergic and GABA-ergic synapses (Bhattacharya et al., 2012; Bhattacharya, 2013). The motivation for these research has been to take a step forward in building computational tools that can complement experimental research in understanding the underlying cellular mechanisms of anomalous EEG signals in neurological and psychiatric disorders. The main inspiration for this approach has been the work by Destexhe (1994), where the authors state the following when discussing the future benefits of kinetic modeling of synaptic processes (p. 223): “A considerable amount of experimental data is available from measurements of the average activity of populations of brain cells: recordings of electroencephalogram, LFPs, magnetoencephalograms, optical recordings, magnetic resonance images, etc. It would be interesting to attempt to establish a relationship between such global measurements and dynamics at the molecular level.” Our work on NMMs embedded with synaptic kinetics showed a high sensitivity to the neurotransmitter concentration, forward and reverse rates of reactions during synaptic transmission, and the membrane conductance of the cell populations. Besides enabling the correlation of lower-level synaptic attributes to population-level dynamics in the model, the approach provided a 10-fold decrease in computational times compared to classic NMMs. The study presented here uses a neural mass computational model of the LGN implementing kinetics of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) and γ-amino-butyric acid (GABA) neuro-receptor mediated synapses, and is an extension of the model in Bhattacharya (2013).

For a long time, the only role of the thalamus was thought to be relaying sensory information to the cortex. The significant advances in our current understanding of other important functions of the thalamic circuitry, especially its integral role in cortical dynamics, is based upon some crucial and pioneering experimental investigations of the thalamic cells in vitro (McCormick and Prince, 1987; von Krosigk et al., 1993). Based on de-corticated (disconnected from the cortex) thalamic slices from the LGN of mammals and rodents, these works established an in-depth understanding of the intrinsic dynamics of the intra-thalamic cell populations (McCormick and Pape, 1990; Steriade et al., 1993; Pape and McCormick, 1995; Bal et al., 1995) and its fundamental role in cortico-cortical communication (Sherman and Guillery, 2002). In this work, our objective is to gain a thorough understanding of the synaptic dynamics of the IN afferents and efferents and their role in the LGN dynamics. In fact, a pioneering work identifying the significant effect of IN on the TCR cell dynamics constitutes of an in vitro experimental study on LGN slices from the rat brain (Crunelli et al., 1988). The model presented in this work endeavors to simulate the dynamics of a de-corticated LGN, similar to these early experimental studies.

The model output dynamics are observed under conditions of deletion of the IN synaptic pathways in the LGN circuit and corresponding to a set of “base” parameter values in the model as well as with parametric deviations. Our motivation toward such a model simulation set up is two-fold: to investigate the functional significance of the IN population on LGN oscillatory dynamics, considering that these cells are ignored in current model-based studies of thalamocortical oscillations; to investigate dynamical transitions in the model output that may arise with a dysfunction of the IN cell circuitry, which in turn may be caused by some abnormal pathological conditions in the network. The model input is a white noise simulating background spiking activity of the retinal neurons with eyes closed i.e., absence of visual input. Our results show that: First, if the IN population is disconnected from the network, both TCR and TRN population outputs display synchronous oscillations with waxing-and-waning of amplitude and power spectral peak frequency within the alpha band at ≈ 11 Hz: the hallmarks of EEG and LFP in an awake and resting state with eyes closed. Similar alpha rhythmic synchronous oscillations are also seen for reduced glutamatergic retinal input to the IN population (i.e., when IN is included in the network but receiving reduced visual input). Second, a “smooth” transition of the power spectra from the alpha to theta band is effected by a progressive decrease of neurotransmitter concentration in the synaptic clefts. In contrast, if the inhibitory effect of the IN is reduced by decreasing its efferent synaptic connectivity to the TCR below a certain threshold, a “dramatic” (as opposed to “smooth”) transition from alpha to theta band is observed in the model output when the levels of neurotransmitter concentration are lower than the set base value; the TCR and TRN enter a state of synchrony, showing waxing-and-waning amplitude patterns within the theta band. Third, increasing the ensemble leakage conductance of the IN cell population hyperpolarizes its average membrane potential, and both TCR and TRN time-series show synchronous oscillations in the limit-cycle mode and within the theta band. Thus, our results make a strong case for the IN as an integral part of the LGN circuitry in the awake cognitive state, maintaining an overall homeostasis in the system by preventing dramatic oscillatory state transitions in the TCR and the TRN, which appear when specific sets of synaptic attributes linked to the IN are altered.

As a case study, the model is tested for simulating Steady State Visually Evoked Potentials (SSVEP)—brain signals corresponding to flickering visual input at a constant or slow-varying frequency that can be observed through EEG. The frequency of the SSVEP signal is easy to control through the visual input, thus making it popular as a tool to study both lower-level and higher-level vision (see Norcia et al., 2015 for a review). The primary advantage of SSVEP is that the response frequency follow that of the input stimulus, making it possible to distinguish between multiple stimuli at different frequencies. It is therefore not surprising that SSVEP is growing in popularity for applications in clinical neuroscience aimed at understanding brain diseases (see Vialatte et al., 2010 for a review). In addition, Brain-Computer Interface (BCI) research have also used SSVEP for advancing current state-of-the-art applications (Guger et al., 2012). Thus, it seems appropriate that a model simulating LGN dynamics be tested for SSVEP response corresponding to a cognitive brain state. A visual stimulus input to the model is simulated by superimposing the (above-mentioned) white random noise with periodic impulse trains at a single frequency within the range 5–50 Hz. For all frequencies in this range, the fundamental frequency of the model output reflects the frequency of the input impulse train and has distinct harmonic components, thus agreeing with existing literature on SSVEP characteristics (Vialatte et al., 2010; Norcia et al., 2015). Furthermore, the IN and TCR cells in the model show similar power spectra characteristics and their time-series are in-phase. However, the TRN time-series is in anti-phase (approximately) with TCR. There is a distinct change in the model behavior when the IN is disconnected from the circuit: the TRN and TCR time-series are now in-phase (approximately) and the resonant frequency of the circuit dominates the power spectra, thus implying a reduced effect of the periodic stimulus input.

The model structure, parameterization, and simulation methods are discussed in Section 2. Results are presented in Section 3 along with discussions on their implications and model-based predictions. We conclude in Section 4 and mention future directions that will build on this work.

2. Materials and Methods

A schematic of the model used in this work is shown in Figure 1. The synaptic layout of the model is based on the intra-thalamic connectivity presented in Sherman (2006) and Sherman (2007), which in turn are informed by experimental data obtained from the LGN of mammals and rodents and are as reported in Van Horn et al. (2000), Sherman and Guillery (2001) and Jones (2007). In addition, recurrent feedback in both IN and TRN (Huntsman et al., 1996) populations are observed in the experimental studies and are included in the model.

FIGURE 1

Figure 1. Schematic of the neural mass model of the thalamic Lateral Geniculate Nucleus (LGN) consisting of three cell populations viz. the Thalamocortical Relay (TCR), the Interneurons (IN) and the Thalamic Reticular Nucleus (TRN). The excitatory synapses are modeled with AMPA kinetics while the inhibitory synapses are modeled with GABA_A and GABA_B kinetics. The output of the TCR population is taken as the model output, while the model input is assumed to be mean activity of the retinal spiking neurons (RET). The model “base” (reference) parameter values are tuned for dominant alpha rhythms in the model output and are given in Tables 1, 2. The intra-thalamic population connectivities are based on experimental findings (Sherman and Guillery, 2001) and are discussed in Section 2.2. Both AMPA and GABA_A are ionotropic (fast) synapses and are modeled with one first order differential equation (Equation 2). The GABA_B is a metabotropic (slow) synapse and is modeled with two sets of first order differential equations (Equations 3–5).

The output of the TCR cells form the main source of sensory information to the visual cortex. Furthermore, in simultaneous studies on LGN and cortical outputs, the LFP from the TCR cells are observed to have a high coherence with EEG from the occipital scalp electrode (Lopes da Silva et al., 1974; Hughes and Crunelli, 2006; Bastos et al., 2014). Thus, in this model of the LGN, the time-series of the TCR cells is considered as the output and is hereafter referred to as the “model output.” Input to the model represent the mean activity of the retinal spiking neurons (the afferent population to the LGN) when the brain is in an awake state with eyes closed.

The retinal spiking neurons (ganglion cells) make excitatory (glutamatergic) synapses with the TCR and IN population of the LGN that are mediated by both fast ionotropic (iGluR) and slow metabotropic (mGluR) glutamate neuro-receptors (Cox et al., 2003; Wang et al., 2011b). However, Pape and McCormick (1995) report that in the presence of iGluR, mGluR have minimal effect on the IN membrane potential and their synaptic activities. This observation is subsequently confirmed by Govindaiah and Cox (2006). Hence, this pathway is often ignored in synaptic circuits of the LGN (Sherman, 2006). We follow the latter work, and consider only the ionotropic AMPA (iGluR) neuro-receptors mediated synaptic efferent pathways from the retina to both IN and TCR populations. The IN cells make feed-forward inhibitory (GABA-ergic) synapses on the TCR cells mediated by the ionotropic GABA_A neuro-receptors. The TRN cells receive excitatory synapses from the TCR population mediated by ionotropic AMPA neuro-receptors, and send inhibitory feedback to the TCR cells mediated by both ionotropic GABA_A and the metabotropic GABA_B neuro-receptors. In addition, the cells in both TRN and IN make feedback connections on their respective self populations, mediated by the ionotropic GABA_A neuro-receptors.

It is worth mentioning here that all thalamic nuclei that process sensory information are reported as having a similar architecture to that of the LGN (Sherman and Guillery, 2001; Saalmann and Kastner, 2011). Thus, the model in Figure 1 can also be used for simulating the thalamic dynamics corresponding to other sensory pathways. The mathematical framework for the model is mentioned in Section 2.1; model parameterization is discussed in Section 2.2; simulation methodologies are mentioned in Section 2.3. The reference values, referred to in this work as “base values,” of the model parameters are mentioned in Tables 1, 2.

TABLE 1

Table 1. (A) Data for the forward (α) and reverse (β) rates of synaptic transmission is according to the range mentioned in Golomb et al. (1996) and Destexhe et al. (1994).

TABLE 2

Table 2. Base values of the synaptic connectivity parameters C_uvw in Equation (6) are derived from experimental data on LGN of mammals and rodents (Van Horn et al., 2000; Sherman and Guillery, 2001; Jones, 2007) (see Section 2.2 for a brief overview).

2.1. Model Equations

The model in Figure 1 is defined by the set of first order differential equations in Equations (1)–(8). The neurotransmitter concentration in the synaptic cleft ([T]) is a function of the mean membrane potential of the pre-synaptic population, V_pre, and is simulated with a sigmoid defined in Equation (1):

$\begin{array}{cc}\begin{array}{ll}[T](V_{pre})=\frac{T_{max}}{1+e^{-\frac{V_{pre}-V_{thr}}{\sigma }}}, \hfill & \hfill \end{array}& (1)\end{array}$

where V_thr is the threshold voltage when the neurotransmitter concentration crosses the 50% of its maximum value T_max, and σ is the steepness parameter of the sigmoid. An increase of neurotransmitter concentration in the synaptic cleft can be simulated in the model by decreasing V_thr and increasing σ.

Equation (2) defines the dynamics of the ionotropic synapses in the model viz. those mediated by the AMPA and GABA_A neuro-receptors. The variable r defines the proportion of open ion-channels on the post-synaptic population caused by the binding of the glutamatergic and GABA-ergic neurotransmitters with the AMPA and GABA_A neuro-receptors respectively.

$\begin{array}{cc}\begin{array}{ll}\frac{dr(t)}{dt}=\alpha .\left[T\right]\left(V_{pre}\right).\left(1-r(t)\right)-\beta .r(t), \hfill & \hfill \end{array}& \left(2\right)\end{array}$

where α and β refer to the forward and reverse rates of chemical reactions respectively.

Equations (3)–(5) define the dynamics of the metabotropic synapses in the model viz. those mediated by the GABA_B neuro-receptors. They activate G-proteins which in turn act as the “secondary messengers” and initiate the opening of ion channels.

$\begin{array}{cc}\frac{dR(t)}{dt}={\alpha }_1.[T]\left(V_{pre}\right).\left(1-R(t)\right)-{\beta }_1.R(t)& \left(3\right)\end{array}$

$\begin{array}{cc}\begin{array}{c}\frac{d[X](t)}{dt}={\alpha }_2.R(t)-{\beta }_2.[X](t)\end{array}& (4)\end{array}$

$\begin{array}{cc}\begin{array}{c}r(t)=\frac{{[X]}^n(t)}{{[X]}^n(t)+K_d}\end{array}& \left(5\right)\end{array}$

where R is the fraction of activated GABA_B receptors, which acts as a catalyst in activating the secondary-messenger G-protein (guanine nucleotide binding proteins); [X] is the concentration of the activated G-protein; r is the fraction of open ion channels caused by binding of [X] with independent binding sites; α_{1, 2} and β_{1, 2} are the forward and reverse binding rate constants respectively; n is the number of bound receptor sites and K_d is the dissociation constant of binding of [X] with the ion channels.

The resulting post-synaptic current, I_psc(t), is defined in Equation (6):

$\begin{array}{cc}\begin{array}{ll}{I}_{psc}(t)=C_{uvw}.g_{max}^{syn}.r(t).\left(V_{psp}(t)-E_{rev}^{syn}\right), \hfill & \hfill \end{array}& \left(6\right)\end{array}$

where $g_{max}^{syn}$ and $E_{rev}^{syn}$ are the maximum conductance and reverse potential respectively and their values depend on the mediating synapse syn∈{AMPA, GABA_A, GABA_B}; V_psp is defined in Equation (7) and is the ensemble post-synaptic membrane potential; C_uvw is a normalized figure that represents the percentage of the synaptic contacts made on the post-synaptic cell population u by the pre-synaptic cell population v, and w represents the sign of the synapse i.e., excitatory or inhibitory.

$\begin{array}{cc}\begin{array}{ll}{\kappa }_m.\frac{d{V}_{psp}(t)}{dt}=-{{{\displaystyle \sum }}^{\text{}}}^{\text{}}{I}_{psc}(t)-I_{leak}(t), \hfill & \hfill \end{array}& (7)\end{array}$

where κ_m is the ensemble membrane capacitance of the post-synaptic cell population,

The parameter I_leak in Equation (7) is the ensemble membrane leak current of the post-synaptic cell population and is defined in Equation (8):

$\begin{array}{cc}\begin{array}{ll}{I}_{leak}(t)=g_{leak}\left(V_{psp}(t)-E_{leak}\right), \hfill & \hfill \end{array}& \left(8\right)\end{array}$

where g_leak and E_leak are conductance and reverse potential respectively corresponding to “non-specific” leak (Golomb et al., 1996; Suffczyński et al., 2004) in the ensemble membrane of the post synaptic cell population.

2.2. Model Parameterization

The model input is simulated by computer generated random noise that has a low variance and a white frequency spectrum. The resting state membrane potential of both excitatory cell populations of the model viz. the retinal cells and TCR are set to −65 mV. Thus, the mean of the retinal noisy input is −65 mV and the standard deviation is set to 2 mV² to reduce the stiffness of the solution (in Matlab) for the set of differential equations defining the model in Equations (1)–(8).

The base values of the synaptic connectivity parameters in the model are based on physiological data (Van Horn et al., 2000; Sherman and Guillery, 2001; Jones, 2007) and mentioned in Table 2 along with the corresponding nomenclature for specific connectivities used in this work. For the purposes of this work in simulating the LGN dynamics of a de-corticated thalamus (discussed in Section 1), we ignore all afferent cortico-thalamic inputs to the model. A brief overview of the literature survey on intra-thalamic cell connectivity is provided below:

• TCR afferents: Data from the dorsal cat LGN (LGNd) (Van Horn et al., 2000) suggest that the TCR receive ≈ 7.1% of their inputs from the retinal ganglion cells (C_tre), while ≈ 30.9% of their inputs are from inhibitory sources viz. IN and TRN. However, and to the best of our knowledge, there is no data available that distinguish between the afferent synaptic terminals from the IN and TRN. Thus, in this work, the GABA_A afferents from TRN and IN ($C_{tni}^a$ and C_tii respectively) and the GABA_B afferent from the TRN ($C_{tni}^b$) are tuned so that the sum total of all inhibitory afferents on the TCR is 30.9%. The remaining ≈ 62% of the connections are from the cortex as well as other sub-cortical sources, and are ignored in the present work for brevity.

• IN afferents: A study made on the cat LGNd in 1991 suggest that the IN cells receive around 25% synapses from the retinal spiking neurons, 37% from other inhibitory sources including themselves, while 26% synapses are from the cortex. However, according to a more recent study (Van Horn et al., 2000), these figures are reported as 47.4, 23.6, 29% respectively. On the other hand, data from the LGNd of a squirrel monkey (primate) indicates that the IN cells receive an equal proportion of each of the three categories of synaptic terminals (Jones, 2007). To maintain consistency with the source of data for TCR population, we follow the data by Van Horn et al. (2000) and set the retinal input (C_ire) and self-inhibitory (C_isi) connectivities in the IN as 47.4 and 23.6% respectively. It may be noted that the IN circuitry in the LGN has a unique triadic spiking arrangement (Sherman, 2004) consisting of dendrites that are not only post-synaptic to the retinal cells but are also pre-synaptic to the TCR cells as well as to themselves (dendro-dendritic synapses). These are referred to as the F2 terminals of the IN while the usual pre-synaptic axonal terminals are referred to as F1. However, the exact distribution of F1 and F2 terminals are not yet available from physiological studies; thus the above-mentioned figures for synaptic afferents to the IN refer to the combined numbers of both types of terminals.

Experimental observations of IN cell dynamics do mention an excitatory feed from the TCR to these cells (Crunelli et al., 1988; Zhu et al., 1999a; Lörincz et al., 2008), however, these were speculations based on cell behavior as opposed to cell physiology. On the other hand, experimental studies on IN physiology suggest two specific cell types in the LGN (Cox et al., 2003): the intra-layer IN cells that do not receive any afferents from the TCR cells; the inter-layer IN cells that do receive inhibitory feedback from the TCR cells, and are often thought to be “stray” cells of the TRN. In the present work, we consider the intra-layer IN cell population only; thus, as in our previous work (Bhattacharya et al., 2011b), the IN population do not receive any synaptic afferents from the TCR.

• TRN afferents: Both thalamocortical and corticothalamic synapses on the TRN sector corresponding to the rat LGNd and visual cortex are excitatory (glutamatergic) in nature and constitute ≈ 30–40% and ≈ 50% respectively of the total synapses; the remaining up to 25% of the synapses are from other inhibitory sources including neighboring intra-population cells (Jones, 2007). In our model, we maintain the TCR afferent connectivity (C_nte) as 35%, and self inhibitory connectivity (C_nsi) as 20%.

The parameter [T]_max in Equation (1) is well approximated by 1 mM (Destexhe et al., 1998), while the base values for V_thr and σ are obtained by trial simulation studies on the model (further elucidated in Section 3.1) and are set to −32 and 3.8 mV respectively. The capacitance κ_m is set at 1 μF/cm².

It is worth noting that while the model output time-series in Bhattacharya (2013) demonstrates rich dynamics corresponding to alterations in AMPA and GABA_A synaptic parameters and conforming to experimental observations, the period of oscillations in the model were below that of normal brain oscillations as seen in LFPs and EEGs. Appropriate modifications are made in this work to overcome the limitations of the previous work, and the parameters α, β and $g_{max}^{syn}$ in the Table 1A as well as g_leak in Table 1B reflect the modified unit scales.

The variable V_psp for TCR, TRN, and IN are initialized to the respective resting state values as in Table 1B; the variables “r, R, X” in Equations (2)–(5) are initialized to an arbitrarily small value 0.001 (and the ODE solutions do not show any dependency on the initial values). Please see the Table 1 legend for further information on parameter sources.

2.3. Simulation Methods

The ODEs are solved using the 4^th/5^th order Runge-Kutta-Fehlberg method (RKF45) in Matlab for a total duration of 40 s at a resolution of 1 ms. The output voltage time series is averaged over 20 simulations, where each simulation runs with a different seed for the noisy input. For frequency analysis, an epoch from 10 to 39 s of the output signal is sampled every 1 ms (1000 Hz) and bandpass filtered between 1 and 100 Hz with a Butterworth filter of order 10. The filtered signal is then transformed using 4-point FFT and power spectral density derived using the Welch periodogram. The power plots in Section 3 show the averaged power spectral density over 20 simulations. The bar plots show the total frequency content within each of the four frequency bands viz. delta (1–3.5 Hz), theta (3.75–7.5 Hz), alpha (7.75–13.5 Hz), and beta (13.75–20). Short Time Fourier Transform (STFT) is carried out on the averaged membrane potential of each cell population with a Hamming window of duration 1 s and overlap of 50%.

Simultaneous variation of parameters is an inherent feature in the brain as in all dynamical systems in nature. The kinetic modeling approach based NMM is computationally efficient in terms of both time and memory usage compared to the alpha-function based NMMs, as well as possessing greater biophysical plausibility. This approach is adopted in model presented here and has allowed for extensive simulation trials of the model involving simultaneous parameter variation in order to set the base parameter values in Table 1, as well as to study model behavior as discussed in the following text.

3. Results and Discussion

The results presented herewith make a comparative study on the effects of incorporating the IN in the LGN model with the case when it is excluded from the circuit. In Sections 3.1–3.4, we present and discuss parametric deviations that are observed to effect a state transition in the model output time-series and a shift in its power spectra. In Section 3.5, we examine a case study by simulating SSVEP in the model.

3.1. Causality of the Neurotransmitter Concentration: Setting Base Parameters

Figures 2A,B show an overview of the model output power within the alpha and theta bands respectively for simultaneous variation of σ in the range 3–4 mV at a resolution of 0.1 mV, and V_thr in the range −30 to −35 mV at a resolution of 1 mV, simulating a gradual increase of neurotransmitter concentration [T] (the reader may please refer to Equation 1). The power within the theta band has a left skew for lower values of both parameters; there is a drop in the overall band power for increasing values of V_thr. Power within the alpha band saturate for higher values of σ across all values of V_thr; however, the absolute power is less than that within the theta band for lower values of sigma. In general, for a fixed value of V_thr (σ), increasing (decreasing) σ (V_thr) causes an increase in [T] and effects a progressive shift of the model output power spectra from a dominant theta band power to a dominant alpha band power accompanied by a consistent increase in beta band power.

FIGURE 2

Figure 2. The model output power spectra for simultaneous variation of parameters V_thr and σ and within the (A) alpha and (B) theta frequency bands. (C) Histogram of the model output power within the theta and alpha bands at V_thr = −32 mV and for values of σ between 3 and 4 mV at a resolution of 0.2 mV. The power within the alpha band is greater than that within the theta band for σ ⩾ 3.6. For 3.6 ⩽ σ ⩽ 4, the alpha peak has a maximal difference with the theta peak at σ = 3.8 mV. Thus, the base values of V_thr and σ are set at −32 and 3.8 mV respectively.

As mentioned in Section 2.2, the base values of both parameters V_thr and σ are set by visual inspection of the results in Figure 2 such that the model output power is dominant within the alpha band. In Figure 2A, the image show minimal skewness for V_thr = −32 mV. In Figure 2C, the case for V_thr = −32 mV is shown: we note that at σ = 3.8, the alpha band power is greater than that of the theta band, and the difference between the two band powers is maximum for the tested range; thus, V_thr = −32 mV and σ = 3.8 mV are set as base parameter values in the model.

Discussion

The gradual alpha to theta band transition with reduced levels of neurotransmitter concentration in the model mimics the EEG marker of transition from a state of wakefulness to a state of drowsiness (Hughes and Crunelli, 2005). It is now well understood that the release of neurotransmitters in the synaptic cleft is mediated by calcium current; an in-depth study by Llinas et al. (1981) on squid giant axon suggests transmitter depletion in the cleft due to reduced pre-synaptic activity, which in turn lead to reduced post-synaptic activity. In the LGN, McCormick (1992) have demonstrated that thalamic slices in vitro display membrane depolarization and cessation of rhythmic oscillations on application of cholinergic and noradrenergic neuromodulatory inputs. This implies that conversely, depletion of these neuromodulators may lead to slow synchronous rhythms in the LGN. More recently, Zhao et al. (2002) have shown that dopaminergic neuro-receptors affect information transmission efficiency in the LGN by neuromodulation of the IN cells. Furthermore, their study suggest a stronger suppression of the TCR cells with blockage of GABA_A neuro-receptors. Our model-based study show an abrupt (as opposed to progressive) transition from alpha to theta rhythm when the IN feed-forward inhibition of the TCR cells are deleted in the circuit (see Section 3.2). Overall, it may be speculated that the state transition in the model may mimic the effect of IN cells due to neuromodulatory influences from extra-thalamic sources.

3.2. Causality of the Interneurons in the LGN Circuit

The time series and power spectra of the model cell populations when all parameters are at their base values are shown in Figure 3. The membrane potential of the TRN has a low peak-to-peak oscillatory envelope of ≈ 0.3 mV in comparison to ≈ 0.8 mV for TCR and IN. The power spectra of both TCR and IN are broad with peak frequencies at ≈ 8.5 and 13 Hz respectively. The corresponding STFT plots in Figure 4 (C: with IN) show that the power spectra of both populations span the theta, alpha and lower-beta (4–20 Hz) frequency bands. In comparison, the TRN has a sharp power spectral peak at around 7.5 Hz, and the STFT plot in Figure 4 (C: with IN) show a narrow power spectra within the theta to lower-alpha (4–10 Hz) region corresponding to base model parameter values.

FIGURE 3

Figure 3. Time-series (Left) and Power spectral density (Right) plots for the TCR, IN and TRN cell populations corresponding to base parameter values in the model. The dominant frequency of oscillation of both TCR and IN is within the alpha band (8–13 Hz) with peaks at around 8.5 and 13 Hz respectively. Furthermore, both these populations show a wide power spectra spanning the theta (4–7 Hz), alpha and lower-beta (14–20 Hz) bands. The TRN shows a comparatively narrow power spectra spanning the theta and lower-alpha (5–10 Hz) band with a dominant peak within the theta band at ≈ 7.5 Hz. The time-series show-case a representative sample from 15 to 25 s of the total simulation time of 40 s.

FIGURE 4

Figure 4. (A) The time-series plots show the phase-locked waxing-and-waning oscillatory patterns within the alpha band, (hallmarks of EEG and LFP in a quiet awake state with eyes closed, and resting) in both TCR and TRN, and in the absence of IN in the circuit. The average membrane potential from 15 to 25 s of the total simulation time of 40 s is show-cased and may be compared to the time-series plots in Figure 3, corresponding to when the IN is included in the circuit. (B) The bar plots for TCR (Left) and TRN (Right) show a comparison between when the IN is in the circuit and when the IN is disconnected: In the former case, the TCR and TRN have different power levels within the specified frequency bands. The dominant frequency of the TCR is within the alpha band, while the dominant frequency in the TRN is within the theta band. In the latter case, i.e., with the IN is removed from the circuit, the proportion of power within the delta, theta and beta bands in both TCR and TRN are negligible in comparison to the power within the alpha band. (C) The STFT plots in the lower panel (without IN) show a fairly stationary signal for both TCR and TRN with a peak at ≈ 12.5 Hz in the absence of IN. In comparison, when the IN is included in the circuit, the STFT plots (upper panel) for the TCR indicate a broad power spectra spanning the theta, alpha and lower-beta bands and is similar to that of the IN. However, STFT plots for the TRN output lies within the theta to lower-alpha region (4–10 Hz), implying suppressed dynamics in the presence of IN.

Next, the IN cell population is disconnected from the network by making C_tii = 0. This condition simulates the blocking of GABA-ergic feed-forward inhibition of the TCR due to neuromodulatory effects on the IN by several neurotransmitters such as acetylcholine (McCormick, 1992) and dopamine (Zhao et al., 2002). The time series output of both TCR and TRN in Figure 4A shows a distinct bifurcation with synchronized waxing-and-waning of amplitude within the alpha band frequency. The mean membrane potential of both populations show depolarization and increased peak-to-peak oscillation compared to their respective counterparts in Figure 3. The bar plots of both TCR and TRN in Figure 4B demonstrate the dramatic decrease in theta band power and the distinct dominance of alpha band power when the IN is removed from the circuit. The STFT plots for both TCR and TRN in Figure 4 (C: with IN) indicate a stationary power spectra with a dominant frequency within the alpha band. Furthermore, in the absence of IN, decrease of neurotransmitter concentration in the model effects a dramatic slowing of the model output: thus, while the time series still retain the amplitude-waxing-and-waning pattern in synchrony with the TRN, the dominant frequency of oscillation is within the theta band (not shown here; please refer to discussion in Section 3.1).

Discussion

The model results show a distinct suppression of TRN activity in the presence of IN in the circuit; the corresponding time series is low-amplitude with a broad power spectral content spanning the alpha to lower-beta regions, similar to EEG and LFP in awake cognitive state. Removal of IN effects feed-forward disinhibition of the TCR and generation of synchronous waxing-and-waning oscillations. Experimental studies have shown that the waxing-and-waning patterns in EEG, both in alpha rhythms (amplitude) as well as in sleep spindles (frequency), are the result of the feed-forward and -back interactions between the TCR and the TRN. As both patterns underpin a state of “low vigilance” (Saalmann and Kastner, 2011) i.e., devoid of task-related activity and reduced information transmission efficacy in the retino-geniculo-cortical pathway, it is suggested (Saalmann and Kastner, 2011) that a state of increased information transmission (as opposed to a state of low vigilance) is facilitated in the circuit by suppressing the inhibitory influence of the TRN. The model validates these experimental observations and further predicts that feed-forward inhibition from the IN play an important role in cognitive activities by suppressing the feedback inhibitory effects of the TRN cells.

Another interesting observation by Halassa et al. (2014) reports better task-related performance when the inhibitory efferents of the TRN are suppressed by application of an “external stimulus”—we hypothesize that the “external stimulus” may have simulated the indirect effect of IN on the TRN. Overall, the model predictions underpin a causal role of the IN in transitions between awake cognitive state to that of diminished cognition.

In addition, the model results show that common mechanisms in the LGN cause alpha to theta band transitions, thus agreeing with experimental observations suggesting that common neuronal mechanisms underlie EEG and LFP of alpha and theta rhythms (Hughes et al., 2004).

3.3. Varying the Connectivity Parameters

The synaptic connectivity parameters in the model are varied around their base values, the details of which are mentioned in Table 3 along with a brief overview about their effects on the model behavior. The results are presented and discussed below:

• When the self inhibition connectivity parameter of the IN (C_isi) is decreased, the TRN population output is suppressed significantly in spite having no direct connectivity with the IN cells. Overall power in the model output is decreased due to increased inhibition from the IN population, and the power spectra shifts right with relatively more power within the upper alpha (10–13 Hz) and lower beta bands (14–20 Hz) as shown in Figure 5B. Increasing the parameter values above the base value reduces the effect of IN on the TCR, and thereby on the TRN. A narrow peak within the alpha band in Figure 5B indicates an increased effect of the TRN population on the TCR cells.

• For C_isi ≤ 0.15 and as shown in Figure 5A, the overall power decreases with maximum frequency within the theta band, until C_isi≈ 0.05%, when the dominant frequency shifts to the delta band. When C_isi = 0, the dominant frequency is within the delta band.

Discussion

Cortico-cortical and cortico-thalamic modulation is well known to be vital for delta band oscillations corresponding to slow wave sleep and/or reduced cognitive state. In the model, the only instance of delta band oscillations are simulated for the case corresponding to synaptic deletion of the self-inhibitory afferent pathway of the IN that effectively disinhibits the IN, thus increasing the feed-forward inhibition of the TCR.

Based on these observations, we hypothesize that the cortical feedback to the IN cells target the self inhibitory mechanisms in the population, which in turn effects slow wave oscillations in the LGN. While it may not directly affect an alpha band oscillation, the cortical feedback may cause an increased self inhibition in IN cells, which in turn will lead to decreased effect on the TCR, thus allowing the TRN inhibition to dominate the output showing synchronous alpha rhythms.

• A progressive decrease of C_tii from its base value shows a progressive increase of power within both alpha and theta bands (Figure 6). For C_tii ⩽ 5, the TCR and TRN time-series display alpha rhythmic synchronous oscillations with amplitude waxing-and-waning. Increasing the parameter value did not have any significant effect and the model output showed a relatively flat spectrum within the alpha and lower beta bands, indicating an increased inhibition in the circuit; this is demonstrated in Figure 6 for C_tii = 23.

• Decreasing the retinal input connectivity for the IN (C_ire) effects a significant rise of the alpha band power in the TCR output until below a threshold, when the TRN dominates and the output show alpha band waxing-and-waning oscillations. Increasing C_ire from its base value does not affect the model output.

Discussion

In summary, the results show that a decrease in feed-forward inhibition in the LGN causes an increase in power within the alpha and theta bands, implying that at base parameter values, the inhibitory effect of TRN on TCR is suppressed by IN. We have discussed (above) the correspondence between cognitive states and suppression of TRN activity by the IN. Indeed, experimental evidence suggest that the IN feed-forward inhibition serves to enhance the “sensitivity to visual features” (Wang et al., 2011a), and thus increase the overall efficiency of retino-geniculo-cortical information transmission (Dublin and Cleland, 1977; Hirsch et al., 2015). The results support these experimental reports.

The model identifies a vital difference between the two GABA-ergic efferent pathways of the IN—synaptic variation in the efferent pathway to the TCR (C_tii) bears a linear relation to the power spectra, unlike in the self inhibitory efferent pathway (C_isi), where synaptic depletion/enhancement shows a non-linear relation with the power spectra dominant frequency. Thus, it seems likely that these two pathways form parts of different circuitry in the LGN: while the C_isi may be a part of the cortico-thalamic feedback mechanism (discussed above), the C_tii appears to be an intra-thalamic connectivity that aids retino-geniculo-cortical information transmission and regulates brain cognitive states. This hypothesis is further supported by identical effects of synaptic depletion in both the retino-geniculate (C_ire) and the intra-geniculate (C_tii) pathways, causing amplitude waxing-and-waning synchronous oscillations within the alpha band in both TCR and TRN, implying a reduced cognitive state.

• The GABA_A inhibitory feedback from the TRN to the TCR (${\mathbf{\text{C}}}_{\mathbf{\text{tni}}}^{\mathbf{\text{a}}}$) shows minimal effect on model output. However, there is a dramatic effect when the IN is disconnected from the network: The time-series of both TCR and TRN display bifurcation into a high amplitude alpha band limit cycle oscillations for $C_{tni}^a>$ its base value (Figure 7). The limit cycles disappear for $C_{tni}^a⩽15.45$ (base value), and the output signal displays synchronous waxing-and-waning alpha band oscillations (as is expected from our earlier observations with disconnecting the IN from the network).

• Decreasing the self inhibitory connectivity of TRN (C_nsi) from its base value of 20 cause a small but progressive increase of power within the theta band as shown in Figure 8A. However, removing IN from the network with a simultaneous decrease in parameter value induces a bifurcation of the time-series to a limit cycle mode for both TCR and TRN; the case for C_nsi = 10 is shown in Figure 8B, left panel. The corresponding STFT plots in Figure 8B, right panel show second harmonic of the dominant frequency within the alpha band. For both cases i.e., with and without IN, when C_nsi = 0, the dominant frequency of oscillation is within the theta band as seen in Figure 8A.

Discussion

The results are in agreement with the suppression of the TRN output in presence of the IN as discussed in Section 3.2: efferent and afferent synaptic pathways in the TRN in the model have minimal effect when the IN has the dominant inhibitory influence in the LGN circuit. Bifurcation to synchronous waxing-and-waning of amplitude as well as limit cycle oscillations are observed in both TCR and TRN when the IN is absent from the LGN circuitry. Limit cycle oscillations are often associated with abnormal brain behavior for example epilepsy (Suffczyński, 2000). Thus, the model predicts a homeostatic role of the IN in the LGN circuitry by controlling high amplitude synchronous oscillations.

In Saalmann and Kastner (2011), the authors speculate that a disinhibition mechanism adopted by the TRN (by increasing its self inhibitory connectivity) may regulate the appearance/absence of synchronous oscillations in the TCR. However, the model predicts the IN as the “master switch” that “enable” such mechanisms. When the IN synaptic pathway is sufficiently depleted, the IN is “disabled” in the circuitry, and both inhibitory efferent pathways of the TRN show significant role in TCR oscillations. In addition, the model predicts that with deletion of the self inhibitory mechanism in the TRN cells, the TCR output is “stalled” within theta band, an EEG biomarker of drowsy states or reduced cognition. This is contrary to the speculations made in Saalmann and Kastner (2011), where the authors suggest an increased self inhibition in the TRN population as a possible mechanism for facilitating useful thalamocortical information transmission. The disagreement may be due to the lack of a cortico-thalamic pathway in the model, and will need further investigation in future works with enhanced thalamo-cortico-thalamic circuitry.

TABLE 3

Table 3. Synaptic connectivity parameters that show sensitivity to the IN in the model are increased (>) and decreased (<) progressively with reference to their respective base values.

FIGURE 5

Figure 5. (A) Near zero values of the self-inhibitory connectivity of the IN population (C_isi) show dominant frequency within the theta and delta bands. (B) For increased values of C_isi from its base value (23.6), the power spectra is fairly flat with dominant frequency within the alpha and beta bands. For decreasing values of the parameter, there is an increase of power within the alpha band with a slight right-shift of the peak frequency.

FIGURE 6

Figure 6. Increase in power within the alpha band with diminishing values of the feed-forward inhibitory connectivity of the IN on TCR population (C_tii). For values less than around 5, the time-series shows synchronized oscillation of the TCR and TRN. For values higher than the base parameter value of 20, the power spectra is fairly flat with maximum power within the alpha and lower beta bands.

FIGURE 7

Figure 7. The TCR population output displays alpha band waxing-and-waning behavior in synchrony with the TRN population in the absence of IN in the circuit. Under this condition, if the ionotropic GABA_A connectivity from the TRN to the TCR ($C_{tni}^a$) is increased, the model time-series shows a bifurcation into a limit cycle mode. With the IN population incorporated in the LGN circuit, variations of $C_{tni}^a$ values do not have any significant effect on the output.

FIGURE 8

Figure 8. (A) Slowing of the power spectra with reduced value of the self-inhibitory pathway (C_nsi) of the TRN population. (B) In the absence of IN in the LGN circuit, the TCR and TRN population output are synchronized and oscillate in the limit cycle mode for reduced values of C_nsi. The figure shows the case when C_nsi = 10.

3.4. Effects of Leak Conductance

Increased or suppressed potassium leak currents are reported as facilitating hyperpolarization and depolarization respectively of cell membrane potential, and show evidence of the critical role in normal (Goldstein et al., 2001) and abnormal brain-state transitions (Gentiletti et al., in press). Neurotransmitters/neuromodulators such as serotonin and acetylcholine inhibit potassium leak currents to enhance neural excitability. On the other hand, anaesthetics are known to decrease excitability in muscles by increasing leakage conductance. To make a qualitative validation of these observations in our model, we vary the leakage conductances of all three cell populations in the model.

The g_leak for any one cell population is increased to 100 μS/cm² progressively, while the values of the same for the other two cell populations remain at their base values of 10 μS/cm². A simultaneous variation of neurotransmitter concentration parameter [T] is also made by varying σ. The results are summarized below:

• When g_leak = 100 μS/cm² for the IN cells, their mean membrane potential is hyperpolarized, causing a reduced effect on TCR and TRN, both of which show a depolarization. However, for reduced values of the neurotransmitter concentration (σ < 3.4), the TCR cells are hyperpolarized.

• When g_leak = 100 μS/cm² for TRN cells, their mean membrane potential is hyperpolarized and the TCR cells are depolarized, and a smooth alpha to theta transition occurs with decreasing neuro-transmitter concentration.

However, if the IN is removed from the circuit, increasing g_leak for the TRN population causes a depolarization in both TRN and TCR cells. The time series of both cell populations (not shown here) display synchronous alpha rhythmic waxing-and-waning oscillations while the power spectra display harmonics of the dominant frequencies within the alpha band.

• When g_leak = 100 μS/cm² for the TCR cells, both TCR and TRN are depolarized, with a larger beta content in the power spectra for increasing neurotransmitter concentration. Both time series and power spectra display smooth transition to lower frequency bands with decreasing neurotransmitter concentration. Removing the IN did not show any drastic change in the output characteristics.

Discussion

The results imply that for a model tuned to oscillate within the alpha band and consisting of all three LGN cell populations, hyperpolarization of the TCR is effected only when the leakage conductance of the IN population is increased, and under the condition of reduced transmitter concentration; we note that this is similar to the above-mentioned effects of anaesthetics (Goldstein et al., 2001). Furthermore, in their experimental investigation on interneuron physiological characteristics, Pape and McCormick (1995) mention that membrane hyperpolarization is affected by an increase in membrane potassium conductance, which in turn is effected by activating acetylcholine (ACh) receptors. All cell populations of the LGN are known to receive significant cholinergic inputs from the brainstem (Sherman, 2006), which play a dominant role in effecting alpha band oscillations in the LGN (Saalmann and Kastner, 2011). Thus, the effects of increased leak conductance in the model may be a simulation of the causality of brainstem cholinergic inputs to the LGN.

On the other hand, if the synaptic efferent from the IN to TCR is decreased while leak conductance of TRN is increased, both TCR and TRN are depolarized leading to abrupt state transitions to high amplitude synchronous oscillations with dominant alpha band frequency; decreased neurotransmitter concentration under these conditions generate synchronous waxing-and-waning within the theta band. This observation, once again, predicts a significant role of the IN in maintaining a homeostasis in the LGN; it minimizes any unwanted state transitions in the circuit that may be caused by increase of the leakage current.

Overall, our model-based study emphasize the role of leakage potassium currents in the LGN circuitry and identifies the need for further investigation and validation with experimental studies in the context of unwanted and abrupt transitions in the circuit.

3.5. Simulating Steady State Visually Evoked Potentials: A Case Study

As a case study, the LGN model tuned to simulate alpha rhythmic output in an awake state is tested for simulating SSVEP output corresponding to a periodic impulse train simulating visual stimulus. The white noise input to the model as used in Sections 3.1–3.4 is now superimposed by periodic impulse train of frequency f ∈ {5 − 50Hz}. The amplitude of the periodic impulses are set arbitrarily at a value of 10 mV.

Figure 9 shows the time series and power spectral density plots corresponding to an 8 Hz impulse train input:

FIGURE 9

Figure 9. The Steady state visually evoked potential (SSVEP) simulated in the model corresponding to an input periodic impulse train at 8 Hz superimposed on a white random noise. (A) A representative sample of 1 cycle of the total simulation time is shown for all three cell populations. For clarity, the input signal is represented by the impulse trains only i.e., the white noise component is not shown. When the IN is in the circuit (left panel), the TCR and IN population have similar waveform with peaks that are in-phase; both waveforms settle to a steady state prior to the next cycle. The TRN population peak lags that of the TCR and IN by 20 ms. When the IN is disconnected from the circuit (right panel), the TCR waveform is now similar to that of the TRN population, and peaks at a delay of 10 ms compared to when receiving inhibitory afferents from the IN (left panel). (B) The power spectral density plot (left panel) for both TCR and IN show peaks at the fundamental frequency corresponding to that of the periodic impulse train, with harmonics at integral multiples of the fundamental frequency. However, for TRN, the peak frequency is at ≈ 7.5 Hz. With IN feed-forward synaptic pathway deleted from the circuit (right panel), both TCR and TRN population dominant frequency is at ≈ 13 Hz. This do not resemble the input impulse train frequency, rather reflect the inherent circuit frequency with white noise input (see Figure 4). The magnitude of power in the dominant frequency for both TCR and TRN are much higher compared to when IN was in the circuit (the power spectra for IN is included just for ease of comparison). Power within harmonics for both TCR and TRN are negligible relative to that within the dominant frequency.

With all parameters at their respective base values, the TCR and IN time-series are in-phase (Figure 9A, left) with a power spectral density (Figure 9B, left) peaking at the fundamental frequency 8 Hz (periodic input stimulus frequency) for TCR and at the second harmonic 16 Hz for IN. Furthermore, both power spectra show harmonics at integer multiples of the fundamental frequency. In comparison, the TRN time series show a low amplitude of oscillation and peak at a delay of around 20 ms from that of IN and TCR. The power spectral peak for TRN is at 7.5 Hz. (We note that across all frequencies (f) in the tested range for this study, the power spectral peaks for the TRN population consistently lie in the range of ≈ 6–8 Hz when the IN is included in the LGN circuitry).

Next, the IN is removed from the circuit; the amplitude of oscillation for the TRN time-series increases significantly in (Figure 9A, right), which is reflected in the increased magnitude within the dominant frequency in the power spectra (Figure 9B, right). The TCR time-series pattern changes dramatically and is now similar to that of the TRN, albeit at a lead of around 10 ms. The power spectral plot also follows that of the TRN with a peak frequency at 13 Hz, and thus do not reflect the that of the input impulse train. In general, across all f within the tested range and when the IN is removed from the circuit, both TCR and TRN show dominant frequency within the range 12–14 Hz, thus over-riding the fundamental frequency of the input impulse train. For some f across the tested range, the power spectra of TCR population show harmonics. However, the power magnitudes in these harmonics are negligible in most cases; where the TRN population power spectra show harmonics, these are always at integral multiples of the peak frequency in the range 12–14 Hz.

Discussion

A study by Babadi et al. (2010) suggests that around 33% of TCR cells in the LGN receive “locked” feed-forward inhibitory inputs from the IN; “locked” (“non-locked”) refers to the case when TCR and IN receive inputs from same (different) retinal ganglion cells (RGC). Furthermore, output of locked TCR cells bear a high degree of correlation with RGC inputs and is speculated to be a mechanism for increasing LGN response precision (Blitz and Regehr, 2005). In the model, both TCR and IN receive the same retinal input for a particular simulation, and the output of both populations are perfectly in-phase; thus the model output may be thought to simulate locked TCR response observed in experimental studies. While we did not consider this aspect (i.e., locked/non-locked IN cells) in our model design, the validation with experimental studies imply a degree of robustness in the model; however, this needs to be tested thoroughly prior to making claims.

The appearance of alpha rhythmic dominant frequency in the model output power spectra underpin the dominant role of TRN with abolition of IN inhibitory input to the TCR and is consistent with our previous results in Sections 3.2 and 3.3. Moreover, the disappearance of the input impulse frequency in the TCR power spectra indicates a de-correlation with the input; rather, the TCR power spectra is now correlated with that of the TRN. This is an interesting phenomenon in the model considering that the TCR continues to receive direct retinal input. These observations possibly imply a greater significance of the indirect sensory pathway of the TCR via its afferent pathway from the IN. Along these lines, synchronous oscillations in the LGN are suggested as a mechanism for decorrelating retinal input from the thalamic output in order to reduce the efficiency of information transmission in the geniculo-cortical pathway (Saalmann and Kastner, 2011). The model validates these speculations and further predicts that: the dominant role of the TRN is facilitated by what appears to be a separate circuit mechanism that inhibits the IN cells beyond a threshold, so that these are in a “disabled” mode and have minimal inhibitory influence on the TCR cells.

Another interesting implication of the results are the known symptoms of fatigue during SSVEP-based BCI experiments (Cao et al., 2014). The increased level of fatigue in a subject undertaking the experiment can drastically reduce the performance, and the corresponding EEG shows an increased proportion of alpha and theta band power. In our study, deleting the IN pathway to the TCR cells effect a significant rise in the power within both alpha and theta bands. This raises speculation that neuromodulatory circuits inhibiting the IN may play an indirectly role in inducing fatigue corresponding to SSVEP-based BCI tasks.

Overall, the preliminary results on the model simulation of SSVEP validates experimental research reporting significant role of IN in modulating visual input that is relayed to the cortex by the TCR population (Lörincz et al., 2009; Wang et al., 2011b; Hirsch et al., 2015). We hypothesize that any anomaly in the IN circuitry will disturb the normal processing of visual information during an awake cognitive state.

4. Conclusion

We use a neural mass computational model of the thalamic LGN, implemented with synaptic kinetics, to understand the effects of the inhibitory interneuron (IN) population on the LGN dynamics. The IN population constitute around 20–25% of all sensory information carrying thalamic nuclei in mammals, and in the LGN of rats. Specifically in the LGN, the IN cells receive around 47% of their inputs from the retinal spiking neurons. It is thus not surprising that extensive research have investigated the functional significance of the IN afferent and efferent synaptic pathways in precise spatial and temporal transmission of visual information to the cortex (see Hirsch et al., 2015 for a review). In contrast, very little research has looked into the role of the feed-forward inhibitory pathway of the IN in normal and abnormal thalamocortical dynamics of health and disease, for example Crunelli et al. (1988), Zhu et al. (1999a) and Zhu et al. (1999b); these studies have emphasized both physiological and functional importance of the IN population in modulating LGN oscillatory activity. A similar trend is evident in computational model-based research on thalamocortical dynamics of health and disease that largely ignore the feed-forward inhibition by the IN; rather, the emphasis has been on the feedback inhibition by the TRN population. The interest in IN is revived in a recent set of experimental studies by Lörincz et al. (2009, 2008) emphasizing a vital role of the IN in modulating LGN oscillations in an awake cognitive state. In a previous work on a thalamo-cortico-thalamic model (Bhattacharya et al., 2011b), we have shown that the IN plays a role in “slowing” (left shift of the power spectral peak frequency) of alpha rhythms (8–13 Hz) recorded from the thalamo-cortical relay (TCR) cells and in the presence of cortical inputs, simulating an EEG biomarker of Alzheimer's disease. However, the non-triviality of the parameter space in the model puts a constraint on an in-depth understanding of the essential thalamic influence on cortical dynamics and vice-versa. Instead, we suggest that a “bottom-up” modeling approach, where the LGN cell dynamics corresponding to retinal input in both resting and awake conditions are studied independently of the cortical input, may better serve our goal of understanding the underlying cellular mechanisms of the IN circuitry and their role in the overall dynamics of the LGN. Indeed, the pioneering experimental studies that have provided the current fundamental knowledge about the thalamus were based on in vitro studies of LGN slices from mammals and rodents, when disconnected from the cortex. Here, based on insights from our prior model-based work, we follow the precedent set by the early experimental studies and model a de-corticated LGN responding to extrinsic input from the retinal pathway.

Several studies have shown that the LGN dynamics recorded from the TCR cells as LFP are correlated with an electroencephalogram (EEG) recorded from the occipital (visual) cortex. Thus, the time-series and power spectra from the TCR cells are considered as the “model output” and are validated qualitatively with reports of LFP (thalamus) and EEG (cortex) in existing literature. In Sections 3.1–3.4, the sensory input is taken as the background firing activity of the retinal ganglion cells in an eyes closed condition and is simulated with a “white” (flat power spectra) random noise. In Section 3.5, a case study is made to simulate Steady State Visually Evoked Potentials (SSVEP) in the model, a research paradigm that is increasing in popularity owing mainly to a direct correlation of the recorded brain response to a periodic input stimulus (see Norcia et al., 2015 for a review). Here, the model input is simulated with a periodic impulse train at a constant frequency in the range 5–50 Hz and superimposed on a white random noise, thus mimicking flashing LED lights provided as visual stimulus in SSVEP experimental studies.

The goal in this work has been to investigate the role of the afferent and efferent synaptic pathways of the IN population in effecting state transitions in the LGN circuitry. Toward this, systematic parametric variations with respect to their base values are explored in the model. The effect of synaptic depletion in the IN feed-forward pathway to the TCR is tested alongside the parametric variations, the objective being to draw a comparison with previous studies that ignore the IN population. Based on some preliminary results in Bond et al. (2014), our hypothesis in this work has been to understand possible underlying neuro-pathological changes that may “block” the IN inhibitory effects on the TCR cells, which may lead to state transitions in the LGN corresponding to both normal and pathological conditions. Indeed, experimental studies have shown the neuromodulatory effects of acetylcholine, noradrenaline (McCormick, 1992) and dopamine (Zhao et al., 2002) on the LGN that seem to “regulate” (under normal conditions) or affect (under abnormal pathological conditions) its synaptic efficacy and oscillatory patterns; specifically, dopamine has been implicated in affecting information transmission in the LGN by neuromodulation of IN. In this work, we have simulated the effects on neuromodulation by varying the neurotransmitter concentration in the synaptic clefts. However, for brevity in this study, we have maintained a single set of parameters for all neurotransmitter concentrations in the model (please refer to Equation 1); thus, any variation of parameter affects changes in all synaptic processes in the model.

The results validate several experimental reports and may be summarized thus: First, the IN population play a dominant role in efficient information transmission in the retino-geniculate pathway by suppressing the TRN feedback inhibitory effect on the TCR cells. Second, bifurcation to high amplitude synchronous oscillations are effected in the TCR cells with synaptic depletion and/or deletion in the feed-forward inhibitory afferent from the IN. A simultaneous decorrelation of the TCR output with retino-geniculate input is observed, thus simulating a transition from an awake state to a state of reduced cognition. Third, the IN plays a “homeostatic balancing” role in the LGN circuit, and its absence or reduced impact (owing to other related parametric deviation) may be speculated to aid unwanted periodic oscillations in the LGN that are often biomarkers of neuro-psychological disorders, for example “thalamocortical dysrhythmia” in schizophrenia or “slowing of alpha rhythms” in Alzheimer's disease. At the same time, a smooth transmission from alpha to theta band with a reduced inhibitory effect of IN simulates EEG of wake-sleep transitions. (The reader may please refer to Section 3 for detailed discussions, model validations and predictions).

The bidirectional connectivity between the thalamus and the cortex are well known to be fundamental to brain rhythms (Sherman, 2005; Destexhe, 2008). Specifically, the cortico-thalamic feedback is speculated to be a key factor in transition between a drowsy state showing theta rhythmic EEG to a deep sleep stage with delta rhythmic oscillations (Abeysuriya et al., 2014) and is consistent with studies on de-corticated thalamus Timofeev and Steriade (1996). Furthermore, physiological studies have attributed the delta band oscillations induced in the thalamus by cortico-thalamic inputs to a specific type of mGluR mediated synapse (Crunelli et al., 2006). However, the cortico-thalamic inputs are not modeled in this work; thus it is not surprising that there is a distinct lack of delta rhythms in the model output, the exception being the case corresponding to the total synaptic deletion in the self inhibitory pathway of the IN population. Our results have shown a non-linear relationship between the power spectral behavior of the model output and the IN self inhibitory mechanisms (the reader may please refer to Table 3, Figure 5); this non-linearity may imply a modulatory effect on the IN population. Along these lines, an experimental study (Lörincz et al., 2008) report that the impact of IN on thalamic oscillations could be identified only upon application of a cortico-thalamic stimulus. We hypothesize that the modulatory cortico-geniculate feedback to the IN targets its self inhibitory mechanism: a deletion of the pathway would lead to delta oscillations in the LGN. Ongoing work will extend our model by incorporating the thalamo-cortico-thalamic closed loop circuit to test this hypothesis. One other inconsistency in the model is that the GABA_B pathway does not seem to have a significant effect on the model output, and thus does not conform to experimental reports (von Krosigk et al., 1993). This calls for a further investigation of the parameters in this pathway, which is also expected to add to the dynamic repertoire of the model.

In conclusion, our model-based study implicates the IN population as a significant and vital constituent of the retino-geniculo-cortical pathway, regulating the state transitions of both TCR and TRN populations, and maintaining an overall homeostatic balance in the LGN circuitry in a normal awake state; any direct or indirect disruption to its synaptic mechanisms may cause unwanted brain rhythms that are EEG and LFP markers of neuro-psychiatric disorders.

Author Contributions

BB planned and designed the research and wrote the manuscript. TB and DT undertook parts of the technical simulation work during their 3rd year MEng project. LH co-supervised TB for the project work and introduced the concept of SSVEP in the model. SD provided consultation and guidance on validation of model results with EEG data and biological relevance.

Conflict of Interest Statement

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

Acknowledgments

The research is undertaken jointly by the School of Engineering and School of Psychology, University of Lincoln and is supported by the School of Psychology research funds. The authors are grateful to the reviewers for their valuable time and suggestions which helped improve the manuscript. BB is grateful to Piotr Suffczynski for useful discussions on the model during his academic visit to the University of Lincoln as a Distinguished Visiting Fellow, Royal Academy of Engineering, UK, during July–September, 2015.

References

Abeysuriya, R. G., Rennie, C. J., and Robinson, P. A. (2014). Prediction and verification of nonlinear sleep spindle. J. Theor. Biol. 344, 70–77. doi: 10.1016/j.jtbi.2013.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Babadi, B., Casti, A., Xiao, Y., Kaplan, E., and Paninski, L. (2010). A generalised linear model of the impact of direct and indirect inputs to the lateral geniculate nucleus. J. Vis. 10, 1–14. doi: 10.1167/10.10.22

CrossRef Full Text | Google Scholar

Bal, T., von Krosigk, M., and McCormick, D. A. (1995). Synaptic and membrane mechanisms underlying synchronized oscillations in the ferret lateral geniculate nucleus in vitro. J. Physiol. 483, 641–663. doi: 10.1113/jphysiol.1995.sp020612

PubMed Abstract | CrossRef Full Text | Google Scholar

Bastos, A. M., Briggs, F., Alitto, H. J., Mangun, G. R., and Usrey, W. M. (2014). Simultaneous recordings from the primary visual cortex and lateral geniculate nucleus reveal rhythmic interactions and a cortical source for gamma-band oscillations. J. Neurosci. 34, 7639–7644. doi: 10.1523/JNEUROSCI.4216-13.2014

PubMed Abstract | CrossRef Full Text

Bhattacharya, B. S. (2013). Implementing the cellular mechanisms of synaptic transmission in a neural mass model of the thalamo-cortico circuitry. Front. Comput. Neurosci. 7:81. doi: 10.3389/fncom.2013.00081

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, B. S., Coyle, D., and Maguire, L. P. (2011a). “Alpha and theta rhythm abnormality in Alzheimer's disease: a study using a computational model,” in Advances in Experimental Medicine and Biology, Vol. 718, eds C. Hernandez, J. Gomez, R. Sanz, I. Alexander, L. Smith, A. Hussain, and A. Chella (New York: Springer), 57–73.

Bhattacharya, B. S., Coyle, D., and Maguire, L. P. (2011b). A thalamo-cortico-thalamic neural mass model to study alpha rhythms in Alzheimer's disease. Neural Netw. 24, 631–645. doi: 10.1016/j.neunet.2011.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacharya, B. S., Coyle, D., Maguire, L. P., and Stewart, J. (2012). “Kinetic modelling of synaptic functions in the alpha rhythm neural mass model,” in ICANN 2012 Part I, Lecture Notes in Computer Science 7552, eds A. E. P. Villa, W. Duch, P. Erdi, F, Masulli and G. Palm (Verlag; Berlin; Heidelberg: Springer), 645–652. doi: 10.1007/978-3-642-33269-2_81

CrossRef Full Text

Blitz, D. M., and Regehr, W. G. (2005). Timing and specificity of feed-forward inhibition within the lgn. Neuron 45, 917–928. doi: 10.1016/j.neuron.2005.01.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Bond, T., Durrant, S., O'Hare, L., Turner, D., and Bhattacharya, B. S. (2014). Studying the effects of thalamic interneurons in a thalamocortical neural mass model. BMC Neurosci. 15(Suppl. 1):P219. doi: 10.1186/1471-2202-15-S1-P219

CrossRef Full Text | Google Scholar

Cao, T., Wan, F., Wong, C. M., da Cruz, J. N., and Hu, Y. (2014). Objective evaluation of fatigue by eeg spectral analysis in steady-state visual evoked potential-based brain-computer interfaces. Biomed. Eng. 13, 1–13. doi: 10.1186/1475-925x-13-28

PubMed Abstract | CrossRef Full Text | Google Scholar

Cox, C. L., Reichova, I., and Sherman, S. M. (2003). Functional synaptic contacts by intranuclear axon collaterals of thalamic relay neurons. J. Neurosci. 23, 7642–7646.

PubMed Abstract | Google Scholar

Crunelli, V., Cope, D. W., and Hughes, S. W. (2006). Thalamic t-type calcium channels and nrem sleep. Cell Calc. 40, 175–190. doi: 10.1016/j.ceca.2006.04.022

CrossRef Full Text | Google Scholar

Crunelli, V., Haby, M., Jassik-Gerschenfeld, D., Leresche, N., and Pirchio, M. (1988). Cl⁻-and k⁺-dependent inhibitory postsynaptic potentials evoked by interneurones of the rat lateral geniculate nucleus. J. Physiol. 399, 153–176. doi: 10.1113/jphysiol.1988.sp017073

PubMed Abstract | CrossRef Full Text | Google Scholar

David, O., and Friston, K. J. (2003). A neural mass model for MEG/EEG: coupling and neuronal dynamics. Neuroimage 20, 1743–1755. doi: 10.1016/j.neuroimage.2003.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A. (1994). Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. J. Comput. Neurosci. 1, 195–230. doi: 10.1007/BF00961734

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A. (2008). “Cortico-thalamic feedback: a key to explain absence seizures,” in Computational Neuroscience in Epilepsy, eds I. Soltez, and K. Staley (Amsterdam: Elsevier), 184–214.

Google Scholar

Destexhe, A., Bal, T., McCormick, D., and Sejnowski, T. J. (1996). Ionic mechanisms underlying synchronised oscillations and propagating waves in a model of ferret thalamic slices. J. Neurophysiol. 76, 2049–2070.

PubMed Abstract | Google Scholar

Destexhe, A., Mainen, Z., and Sejnowski, T. (1994). An efficient method for computing synaptic conductances based on a kinetic model of receptor binding. Neural Comput. 6, 14–18. doi: 10.1162/neco.1994.6.1.14

CrossRef Full Text | Google Scholar

Destexhe, A., Mainen, Z. F., and Sejnowski, T. (1998). “Kinetic models of synaptic transmission,” in Methods in Neuronal Modelling, eds C. Koch and I. Segev (Cambridge, MA: MIT Press), 1–25.

Dublin, M. E., and Cleland, B. G. (1977). Organisation of visual inputs to interneurons of lateral geniculate nucleus in cats. J. Neurophysiol. 40, 410–427.

PubMed Abstract | Google Scholar

Freeman, W. J. (1975). Mass Action in the Nervous System. New York, NY: Academic Press.

Google Scholar

Gentiletti, D., Suffczynski, P., Gnatkovsky, V., and de Curtis, M. (in press). Changes of ionic concentrations during seizure transitions - a modeling study. Int. J. Neur. Syst. doi: 10.1142/S0129065717500046

PubMed Abstract | CrossRef Full Text

Goldstein, S. A., Bockenhauer, D., O'Kelly, I., and Zilberberg, N. (2001). Potassium leak channels and the KCNK family of two-P-domain subunits. Nat. Rev. Neurosci. 2, 175–184. doi: 10.1038/35058574

PubMed Abstract | CrossRef Full Text | Google Scholar

Golomb, D., Wang, X.-J., and Rinzel, J. (1996). Propagation of spindle waves in a thalamic slice model. J. Neurophysiol. 75, 750–769.

PubMed Abstract | Google Scholar

Govindaiah, G., and Cox, C. L. (2006). Metabotropic glutamate receptors differentially regulate gabaergic inhibition in thalamus. J. Neurosci. 26, 13443–13453. doi: 10.1523/JNEUROSCI.3578-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimbert, F., and Faugeras, O. (2006). Bifurcation analysis of Jansen's neural mass model. Neural Comput. 18, 3052–3068. doi: 10.1162/neco.2006.18.12.3052

PubMed Abstract | CrossRef Full Text | Google Scholar

Guger, C., Allison, B. Z., Großwindhager, B., Prückl, R., Hintermüller, C., Kapeller, C., et al. (2012). How many people could use an ssvep bci? Front. Neurosci. 6:169. doi: 10.3389/fnins.2012.00169

PubMed Abstract | CrossRef Full Text | Google Scholar

Halassa, M. M., Chen, Z., Wimmer, R. D., Brunetti, P. M., Zhao, S., Zikopoulos, B., et al. (2014). State-dependent architectures of thalamic reticular subnetworks. Cell 158, 808–821. doi: 10.1016/j.cell.2014.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirsch, J. A., Wang, X., Sommer, F. T., and Martinez, L. M. (2015). How inhibitory circuits in the thalamus serve vision. Annu. Rev. Neurosci. 38, 309–329. doi: 10.1146/annurev-neuro-071013-014229

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, S. W., and Crunelli, V. (2005). Thalamic mechanisms of EEG alpha rhythms and their pathological implications. Neuroscientist 11, 357–372. doi: 10.1177/1073858405277450

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, S. W., and Crunelli, V. (2006). Just a phase they're going through: the complex interaction of intrinsic high-threshold bursting and gap junctions in the generation of thlaamic α and θ rhythms. Int. J. Psychophysiol. 64, 3–17. doi: 10.1016/j.ijpsycho.2006.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hughes, S. W., Lorincz, M., Cope, D. W., Blethyn, K. L., Kekesi, K. A., Parri, H. R., et al. (2004). Synchronised oscillations at α and θ frequencies in the lateral geniculate nucleus. Neuron 42, 253–268. doi: 10.1016/S0896-6273(04)00191-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Huntsman, M. M., Porcello, D. M., Homanics, G. E., DeLorey, T. M., and Huguenard, J. R. (1996). Reciprocal inhibitory connections and network synchrony in the mammalian thalamus. Science 283, 541–543. doi: 10.1126/science.283.5401.541

PubMed Abstract | CrossRef Full Text | Google Scholar

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybernet. 73, 357–366. doi: 10.1007/BF00199471

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, E. G. (2007). The Thalamus, Vol. I and II. 1st Edn., Cambridge, UK: Cambridge University Press.

Llinás, R., Steinberg, I. Z., and Walton, K. (1981). Relationship between presynaptic calcium current and postsynaptic potential in squid giant synapse. Biophys. J. 33, 323–352. doi: 10.1016/S0006-3495(81)84899-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Llinás, R., Urbano, F. J., Leznik, E., Ramírez, R. R., and van Marle, H. J. (2005). Rhythmic and dysrhythmic thalamocortical dynamics: Gaba systems and the edge effect. Trends Neurosci. 28, 325–333. doi: 10.1016/j.tins.2005.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H, Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of brain rhythmic activity. Kybernetic 15, 27–37. doi: 10.1007/BF00270757

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H., van Lierop, T. H., Schrijer, C. F., and van Leeuwen, W. S. (1973). Organisation of thalamic and cortical alpha rhythms: spectra and coherences. Electroencephalogr. Clin. Neurophysiol. 35, 627–639. doi: 10.1016/0013-4694(73)90216-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Lörincz, M. L., Crunelli, V., and Hughes, S. W. (2008). Cellular dynamics of cholinergically induced α (8–13 hz) rhythms in sensory thalamic nuclei In Vitro. J. Neurosci. 628, 660–671. doi: 10.1523/JNEUROSCI.4468-07.2008

CrossRef Full Text | Google Scholar

Lörincz, M. L., Kékesi, K. A., Juhász, G., Crunelli, V., and Hughes, S. W. (2009). Temporal framing of thalamic relay-model firing by phasic inhibition during the alpha rhythm. Neuron 63, 683–696. doi: 10.1016/j.neuron.2009.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Marreiros, A., Daunizeau, J., Kiebel, S., and Friston, K. (2008). Population dynamics: variance and the sigmoid activation function. Neuroimage 242, 147–157. doi: 10.1016/j.neuroimage.2008.04.239

CrossRef Full Text | Google Scholar

McCormick, D. A. (1992). Cellular mechanisms underlying cholinergic and noradrenergic modulation of neuronal firing mode in the cat and guinea pig dorsal lateral geniculate nucleus. J. Neurosci. 12, 178–289.

PubMed Abstract | Google Scholar

McCormick, D. A., and Pape, H.-C. (1990). Properties of a hyperpolarization-activated cation current and its role in rhythmic oscillation in thalamic relay neurones. J. Physiol. 431, 291–318. doi: 10.1113/jphysiol.1990.sp018331

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A., and Prince, D. A. (1987). Actions of acetylcholine in the guinea-pig and cat medial and lateral geniculate nuclei, In Vitro. J. Physiol. 392, 147–165. doi: 10.1113/jphysiol.1987.sp016774

PubMed Abstract | CrossRef Full Text | Google Scholar

Modolo, J., Thomas, A. W., and Legros, A. (2013). Neural mass modelling of powr-line magnetic fields effects on brain activity. Front. Comput. Neurosci. 7:34. doi: 10.3389/fncom.2013.00034

PubMed Abstract | CrossRef Full Text | Google Scholar

Moran, R., Pinotsis, D. A., and Friston, K. (2013). Neural masses and fields in dynamic causal modeling. Front. Comput. Neurosci. 7:57. doi: 10.3389/fncom.2013.00057

PubMed Abstract | CrossRef Full Text | Google Scholar

Norcia, A. M., Appelbaum, L. G., Ales, J. M., Cottereau, B. R., and Rossion, B. (2015). The steady-state visual evoked potential in vision research: a review. J. Vis. 15, 1–46. doi: 10.1167/15.6.4

PubMed Abstract | CrossRef Full Text | Google Scholar

Pape, H. C., and McCormick, D. A. (1995). Electrophysiological and pharmacological properties of interneurons in the cat dorsal lateral geniculate nucleus. Neuroscience 68, 1105–1125. doi: 10.1016/0306-4522(95)00205-W

PubMed Abstract | CrossRef Full Text | Google Scholar

Pressler, R. T., and Regehr, W. G. (2013). Metabotropic glutamate receptors drive global persistent inhibition in the visual thalamus. J. Neurosci. 33, 2494–2506. doi: 10.1523/JNEUROSCI.3458-12.2013

PubMed Abstract | CrossRef Full Text

Rall, W. (1967). Distinguishing theoretical synaptic potentials computed for different soma-dendritic distributions of synaptic inputs. J. Neurophysiol. 30, 1138–1168.

PubMed Abstract | Google Scholar

Robinson, P. A., Rennie, C. J., Rowe, D. L., and O'Connor, S. C. (2004). Estimation of multiscale neurophysiologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53–72. doi: 10.1002/hbm.20032

PubMed Abstract | CrossRef Full Text | Google Scholar

Saalmann, Y. B., and Kastner, S. (2011). Cognitive and perceptual functions of the visual thalamus. Neuron 71, 209–223. doi: 10.1016/j.neuron.2011.06.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarnthein, J., Morel, A., von Stein, A., and Jeanmonod, D. (2003). Thalamic theta field potentials and EEG: high thalamocortical coherence in patients with neurogenic pain, epilepsy and movement disorders. Thalamus Relat. Syst. 2, 231–238. doi: 10.1017/S1472928803000219

CrossRef Full Text | Google Scholar

Sherman, S. M. (2004). Interneurons and triadic circuitry of the thalamus. Trends Neurosci. 27, 670–675. doi: 10.1016/j.tins.2004.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M. (2005). Thalamic relays and cortical functioning. Prog. Brain Res. 149, 107–126. doi: 10.1016/S0079-6123(05)49009-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M. (2006). Thalamus. Scholarpedia 1, 1583. doi: 10.4249/scholarpedia.1583

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M. (2007). The thalamus is more than just a relay. Curr. Opin. Neurobiol. 17, 417–422. doi: 10.1016/j.conb.2007.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sherman, S. M., and Guillery, R. W. (2001). Exploring the Thalamus. 1st Edn., New York, NY: Academic Press.

Google Scholar

Sherman, S. M., and Guillery, R. W. (2002). Thalamic relay function and their role in cortico-cortical communication: generalizations from the visual system. Neuron 33, 163–175. doi: 10.1016/S0896-6273(01)00582-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Stam, C. J., Pijn, J. P., Suffczyński, P., and Lopes da Silva, F. H. (1999). Dynamics of the human alpha rhythm:evidence for non-linearity? Clin. Neurophysiol. 110, 1801–1813. doi: 10.1016/S1388-2457(99)00099-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Gloor, P., Llinás, R. R., Lopes de Silva, F. H., and Mesulam, M. M. (1990). Basic mechanisms of cerebral rhythmic activities. Electroencephalogr. Clin. Neurophysiol. 76, 481–508. doi: 10.1016/0013-4694(90)90001-Z

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., McCormick, D. A., and Sejnowski, T. J. (1993). Thalamocortical oscillations in the sleeping and aroused brain. Science 262, 679–685. doi: 10.1126/science.8235588

PubMed Abstract | CrossRef Full Text | Google Scholar

Suffczyński, P. (2000). Neural Dynamics Underlying Brain Thalamic Oscillations Investigated with Computational Models. Ph.D. thesis, Institute of experimental physics, University of Warsaw.

Suffczyński, P., Kalitzin, S., and Silva, F. L. D. (2004). Dynamics of non-convulsive epileptic phenomena modelled by a bistable neuronal network. Neuroscience 126, 467–484. doi: 10.1016/j.neuroscience.2004.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, P. N., Wang, Y., Goodfellow, M., Dauwels, J., Moeller, F., Stephani, U., et al. (2014). A computational study of stimulus driven epileptic seizure abatement. PLoS ONE 9:e114316. doi: 10.1371/journal.pone.0114316

PubMed Abstract | CrossRef Full Text | Google Scholar

Timofeev, I., and Steriade, M. (1996). Low frequency rhythms in the thalamus of intact-cortex and decorticated cats. J. Neurophysiol. 76, 4152–4168.

PubMed Abstract | Google Scholar

Van Horn, S. C., Erisir, A., and Sherman, S. M. (2000). Relative distribution of synapses in the A-laminae of the lateral geniculate nucleus of the cat. J. Comp. Neurol. 416, 509–520. doi: 10.1002/(SICI)1096-9861(20000124)416:4<509::AID-CNE7>3.0.CO;2-H

PubMed Abstract | CrossRef Full Text | Google Scholar

Vialatte, F.-B., Maurice, M., Dauwels, J., and Cichocki, A. (2010). Steady-state visually evoked potentials: focus on essential paradigms and future perspectives. Prog. Neurobiol. 90, 418–438. doi: 10.1016/j.pneurobio.2009.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

von Krosigk, M., Bal, T., and McCormick, D. A. (1993). Cellular mechanisms of a synchronised oscillation in the thalamus. Science 261, 361–364. doi: 10.1126/science.8392750

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Sommer, F. T., and Hirsch, J. A. (2011a). Inhibitory circuits for visual processing in thalamus. Curr. Opin. Neurobiol. 21, 726–733. doi: 10.1016/j.conb.2011.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Vaingankar, V., Soto Sanchez, C., Sommer, F. T., and Hirsch, J. A. (2011b). Thalamic interneurons and relay cells use complementary synaptic mechanisms for visual processing. Nat. Neurosci. 14, 224–231. doi: 10.1038/nn.2707

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Wei, Y., Vaingankar, V., Wang, Q., Koepsell, K., Sommer, F. T., et al. (2007). Feedforward excitation and inhibition evoke dual models of firing in the cat's visual thalamus during naturalistic viewing. Neuron 55, 465–478. doi: 10.1016/j.neuron.2007.06.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X.-J., Golomb, D., and Rinzel, J. (1995). Emergent spindle oscillations and intermittent burst firing in a thalamic model: specific neuronal mechanisms. Proc. Natl. Acad. Sci. U.S.A. 92, 5577–5581. doi: 10.1073/pnas.92.12.5577

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Goodfellow, M., Taylor, P. N., and Baier, G. (2014). Dynamic mechanisms of neocortical focal seizure onset. PLoS Comput. Biol. 10:e1003787. doi: 10.1371/journal.pcbi.1003787

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Taylor, P. N., and Baier, G. (2013). Computational modelling of micro-seizures and focal seizure onset. BMC Neurosci. 14(Suppl. 1):P14. doi: 10.1186/1471-2202-14-s1-p14

CrossRef Full Text | Google Scholar

Wendling, F., Bartolomei, F., Bellanger, J. J., and Chauvel, P. (2002). Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur. J. Neurosci. 15, 1499–1508. doi: 10.1046/j.1460-9568.2002.01985.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80. doi: 10.1007/BF00288786

PubMed Abstract | CrossRef Full Text | Google Scholar

Zetterberg, L. H., Kristiansson, L., and Mossberg, K. (1978). Performance of a model for a local neuron population. Biol. Cybernet. 31, 15–26. doi: 10.1007/BF00337367

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Kerscher, N., Eysel, U., and Funke, K. (2002). D₁ and D₂ receptor-mediated dopaminergic modulation of visual responses in cat dorsal lateral geniculate nucleus. J. Physiol. 539, 223–238. doi: 10.1113/jphysiol.2001.012721

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J. J., Lytton, W. W., Xue, J.-T., and Uhlrich, D. J. (1999a). An intrinsic osicllation in interneurons of the rat lateral geniculate nucleus. J. Neurophysiol. 81, 702–711.

PubMed Abstract | Google Scholar

Zhu, J. J., Uhlrich, D. J., and Lytton, W. W. (1999b). Burst firing in rat geniculate interneurons. Neuroscience 91, 1445–1460. doi: 10.1016/S0306-4522(98)00665-4

PubMed Abstract | CrossRef Full Text