Implementing the cellular mechanisms of synaptic transmission in a neural mass model of the thalamo-cortical circuitry

1. Introduction

Neural mass computational models mimicking synchronous behavior in populations of thalamocortical neurons are often used to study brain oscillations (David and Friston, 2003; Suffczyński et al., 2004; Breakspear et al., 2006; Sotero et al., 2007; Deco et al., 2008; Izhikevich and Edelman, 2008; Pons et al., 2010; Robinson et al., 2011; de Haan et al., 2012). The term “neural mass” was coined by Freeman (1975), while the neural mass modeling paradigm is based on the mathematical framework proposed by Wilson and Cowan (1973); each cell population in a neural mass model represents a neuronal “ensemble” of mesoscopic-scale (10⁴–10⁷), which are densely packed in space and work at the same temporal-scale, so that for all practical purposes, they can be mathematically treated as a single entity (Liljenström, 2012), whence “mass”. In a seminal work, da Silva et al. (1974) used a neural mass model of a simple thalamocortical circuitry to simulate EEG (Electroencephalography) alpha rhythms (8–13 Hz). Subsequently, this model has been the basis of several research (Zetterberg et al., 1978; Stam et al., 1999; Suffczyński, 2000; Bhattacharya et al., 2011a), albeit with modifications and enhancements; of special mention is the modification introduced by Jansen and Rit (1995) where the model is expressed as a set of ordinary differential equations (ODE). This modification, in turn, has been the basis of many significant research (Wendling et al., 2002; Grimbert and Faugeras, 2006; Ursino et al., 2010). However, the computational basis of the models remain the same—the conversion from firing rate to membrane potential by excitatory and inhibitory neurotransmitters is simulated by convolution of the input from a pre-synaptic neuronal mass with an exponential function, commonly known as the “alpha function”, proposed by Rall (1967). Although the alpha function is a fair estimate of the synaptic process (Bernard et al., 1994), it does not allow an insight into the underlying cellular mechanisms of synaptic transmission associated with abnormal brain oscillations—an aspect emphasized to be crucial as an aid to research in brain disorders (McCormick, 1992; Basar and Guntekin, 2008). The importance of understanding the neuro-transmission mechanisms in slow wave synchronized as well as spindle oscillations is also discussed in several relevant experimental studies (Steriade et al., 1993; von Krosigk et al., 1993). Moreover, correlating synaptic kinetics with brain oscillatory activity has the potential to aid neuropharmacological advances in treating the diseased brain (Aradi and Erdi, 2006). Along these lines, Destexhe et al. (1998) argue that the alpha function is inappropriate for representing post synaptic events other than the originally proposed post-synaptic potential in spiking neural networks; they propose a kinetic framework as a more biologically plausible method of modeling synaptic transmission compared to the alpha function (Destexhe et al., 2002). The ability of such a modeling framework to capture the physiological properties of synaptic transmission was demonstrated by fitting the model outputs to experimental data from hippocampal slices. Moreover, kinetic modeling is reported to be computationally efficient (Destexhe et al., 1994), a vital prerequisite in large-scale computational models. Subsequently, the kinetic models of neurotransmission was used in several single-neuronal-level model-based studies—to investigate thalamic oscillations (Destexhe et al., 1996) and corticothalamic influence on brain oscillatory activity (Destexhe, 2008); to investigate network synchrony (Breakspear et al., 2003); to simulate synchronous behavior observed during in vitro experimental studies on ferret thalamic slice by Wang and Rinzel (1992), Golomb et al. (1994, 1996) and Wang et al. (1995).

A significant modification to current neural mass modeling framework was proposed by Suffczyński et al. (2004) by applying single-neuronal-level model based techniques. Toward this, they proposed an “ensemble” representation of the membrane conductance and post-synaptic current in a neuronal mass model of the thalamocortical circuitry; an integrator is used to generate the “ensemble” post-synaptic membrane potential. In the work presented here, a similar approach is adopted to implement the kinetic framework of synaptic transmission in neural mass models—each post-synaptic attribute is assumed to be an “ensemble” representation corresponding to a “neuronal mass”. For brevity, only two-state (“open” and “closed”) ion-channels (Destexhe et al., 1998) are considered, the desensitized state is ignored. While two-state models are a significant simplification of the very complex nature of ion channel dynamics in biology, they have shown a remarkable fit to biological data compared to more-than-two-state models (Destexhe et al., 1998, 2002). This work aims to interface an abstraction of the ion channel dynamics, such as the two-state ion channel kinetic models, with an abstraction of the population level neuronal behavior, such as neural mass models. The goal is to enable the correlation of higher-level brain dynamics observed in EEG with cellular-level dynamics.

The work is presented thus: first, the kinetic framework for modeling AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic-acid) and GABA (γ-amino-butyric-acid) receptor mediated synapses is introduced in an existing thalamocortical neural mass model (section 2); second, a qualitative comparison of the model behavior with experimental studies on ferret thalamocortical tissue reported in von Krosigk et al. (1993) as well as to single-neuronal-level model based observations reported in Golomb et al. (1996; section 3) is presented; the lack of a quantitative study is mainly to avoid erratic conclusions as difference in model structure and simulation techniques are bound to induce mismatch in numerical results. The model behavior is observed to be consistent with these studies (von Krosigk et al., 1993; Golomb et al., 1996)—The post synaptic membrane conductance in both the thalamocortical relay (TCR) and thalamic reticular nucleus (TRN) cell population plays a role in effecting a bifurcation in model behavior from spindling mode [oscillations with the characteristic waxing-and-waning pattern seen in early stages of sleep (Steriade et al., 1993; Hughes et al., 2004) as well as in alpha rhythmic oscillations during resting brain state (da Silva et al., 1973)] to a limit-cycle mode (synchronized oscillations as seen in later stages of sleep or during absence seizures). The post-synaptic membrane conductance for both AMPA and GABA in the TRN cell population is responsible for sustaining and modulating spindle oscillations in the model output. Blocking the GABA-ergic synapses in the self-inhibitory loop of the TRN cell population effects a low-frequency synchronized oscillation in the model; this is aided by the secondary-messenger-gated GABA synapses in the TCR cell population. In addition, the reverse rate of transmitter binding plays a role in increasing or decreasing the frequency of synchronized oscillations, besides functioning as a bifurcation parameter, an observation that has not been reported in experimental studies. A comparison of the simulation time of the model with previous research using neural mass models based on alpha functions show a factor of 10 improvement in simulation time. This is a dramatic improvement on computational efficiency and emphasizes the appropriateness of the model proposed herein toward building large-scale software models for investigating neuronal disorders. The observations from this study as well as issues related to the modeling approach are discussed in section 4.

2. Materials and Methods

2.1. From Alpha Function to Kinetic Model: a Brief Outline

A single neuronal mass structure as used commonly in neural mass models is shown in Figure 1 and is defined in Equations (1–5):

$\begin{array}{cc}{h}_{\overline{w}}(t)=\frac{H_{\overline{w}}}{{\tau }_{\overline{w}}}\text{tex}p(-t/{\tau }_{\overline{w}})& (1)\end{array}$

$\begin{array}{cc}{y}_N(t)={\displaystyle \sum h_{\overline{w}}(t)\otimes E_{\overline{w}}^N(t)}& (2)\end{array}$

$\begin{array}{cc}{\ddot{y}}_N(t)=\frac{H_{\overline{w}}}{{\tau }_{\overline{w}}}{E}_{\overline{w}}^N(t)-\frac{2}{{\tau }_{\overline{w}}}{\dot{y}}_N(t)-\frac{1}{{\tau }_{\overline{w}}^2}{y}_N(t)& (3)\end{array}$

$\begin{array}{cc}{V}_P(t)={\displaystyle \sum _{N\in \{1,2,3,\dots n\}}{C}_N.y_N(t)}& (4)\end{array}$

$\begin{array}{cc}{E}_{\overline{w}}^P(t)=S(V_P)=\frac{2e_0}{1+e^{\nu (s_0-V_P)}}& (5)\end{array}$

where $\overline{w}\in \{e,i\}$ represents pre-synaptic neuronal populations which make excitatory (e) and inhibitory (i) synapses on a post-synaptic neuronal population; ${\tau }_{\overline{w}}$ is the time constant and $H_{\overline{w}}$ is the amplitude of the synapse; $E_{\overline{w}}^N(t)$, N ∈ {1, 2 … n} is the firing frequency of an extrinsic or intrinsic cell population that is pre-synaptic to the population P; C_N is a percentage of the total number of synapses from all afferents to P; V_P is the “ensemble post-synaptic membrane potential”; $E_{\overline{w}}^P$ is the “ensemble firing rate” of P and is defined by a sigmoid function where 2e₀ is the maximum firing rate of the population, s₀ is the threshold potential at which the neurons spike and ν is the sigmoid steepness parameter.

FIGURE 1

Figure 1. Block diagram of a single “neuronal mass” in current state-of-the-art neural mass models.

2.1.1. A modified neural mass representation

In a recent work, Suffczyński et al. (2004) modified the neural mass representation of a cell population and introduced post-synaptic current mediated by the ligand-gated glutamatergic receptors AMPA, and the ligand- and secondary-messenger-gated GABA-ergic receptors GABA_A and GABA_B, respectively. The input $E_{\overline{\xi }}^N(t)$, $\overline{\xi }\in \{\text{AMPA},{\text{ GABA}}_A,{\text{ GABA}}_B\}$, is the firing rate of an excitatory (AMPA) or inhibitory (GABA_A and GABA_B) pre-synaptic neuronal population N ∈ {1, 2 …, n}. The model (Figure 2A) is defined in Equations (6–11):

$\begin{array}{cc}{h}_{\overline{\xi }}(t)=H_{\overline{\xi }}\left(\text{exp}\left(-t/{\tau }_{\overline{\xi }}^a\right)-exp\left(-t/{\tau }_{\overline{\xi }}^b\right)\right),{\tau }_{\overline{\xi }}^b>{\tau }_{\overline{\xi }}^a& (6)\end{array}$

$\begin{array}{cc}{g}_{\overline{\xi }}^N(t)={\displaystyle \sum h_{\overline{\xi }}(t)\otimes E_{\overline{\xi }}^N(t)}& (7)\end{array}$

$\begin{array}{cc}\begin{array}{l}{\ddot{g}}_{\overline{\xi }}^N(t)=\frac{1}{{\tau }_{\overline{\xi }}^a{\tau }_{\overline{\xi }}^b}[H_{\overline{\xi }}\left({\tau }_{\overline{\xi }}^a-{\tau }_{\overline{\xi }}^b\right)E_{\overline{\xi }}^N(t)-\left({\tau }_{\overline{\xi }}^a+{\tau }_{\overline{\xi }}^b\right){\dot{g}}_{\overline{\xi }}^N(t)\\ -g_{\overline{\xi }}^N(t)]\end{array}& (8)\end{array}$

$\begin{array}{cc}{I}_{\overline{\xi }}^N(t)=g_{\overline{\xi }}^N(t)\left(V^P(t)-V_{\overline{\xi }}\right)& (9)\end{array}$

$\begin{array}{cc}{\kappa }_m{\dot{V}}^P(t)=-{\displaystyle \sum _{N\in \{1,2,3,\dots ,n\}}{C}_N.I_{\overline{\xi }}^N(t)-I_{\lambda }(t)}& (10)\end{array}$

$\begin{array}{cc}{I}_{\lambda }(t)=g_{\lambda }\left(V^P(t)-V_{\lambda }\right)& (11)\end{array}$

where $h_{\overline{\xi }}(t)$ is the synaptic transmission function with ${\tau }_{\overline{\xi }}^a$ and ${\tau }_{\overline{\xi }}^b$ as the rise and decay times, respectively; $g_{\overline{\xi }}^N$ denote the post-synaptic “ensemble” membrane conductance; $V_{\overline{\xi }}$ is the reversal potential for the synapse mediated by $\overline{\xi }$; V^P is the ensemble post synaptic membrane potential of the population P due to PSC from all pre-synaptic cell populations N ∈ {1, 2 …, n}; κ_m is the ensemble membrane capacitance; C_N is the synaptic connectivity parameter; I_λ, g_λ and V_λ are the ensemble leakage current, conductance and reversal potential, respectively for P. The ensemble firing rate $E_{\overline{\xi }}^P(t)$ is as defined in Equation (5) and is the pre-synaptic firing rate input to other neuronal populations.

FIGURE 2

Figure 2. Block diagram of (A) Suffczyński et al. (2004)'s modification of the structure in Figure 1 by introducing “ensemble” representation of post-synaptic membrane conductance and current. (B) Neuronal mass structure with the kinetic framework implemented for modeling synaptic transmission as an alternative to the alpha function ($h_{\overline{w}}(t)$ in Figure 1). A diagrammatic representation of the ion-channel kinetics during synaptic transmission is presented in Figures 3A,B. (C) A thalamocortical circuitry implementing the modified neuronal mass representation in (B).

2.1.2. Introducing kinetic model of synapses in a neural mass representation

The single neuronal mass structure presented in Figures 1, 2A is modified by replacing the alpha function with kinetic models of AMPA, GABA_A, and GABA_B synapses; the enhanced representation (Figure 2B) is defined in Equations (12–19):

$\begin{array}{cc}{[T]}_{\chi }(V_{\chi })=\frac{T_{\text{max}}}{1+e^{-\frac{V_{\chi }-{\theta }_s}{{\sigma }_s}}}& (12)\end{array}$

$\begin{array}{cc}\frac{d{r}_{\chi }^{\overline{{\xi }_1}}(t)}{dt}={\alpha }^{\overline{{\xi }_1}}{[T]}_{\chi }\left(1-r_{\chi }^{\overline{{\xi }_1}}(t)\right)-{\beta }^{\overline{{\xi }_1}}{r}_{\chi }^{\overline{{\xi }_1}}(t)& (13)\end{array}$

$\begin{array}{cc}\frac{d{R}_{\chi }^{\overline{{\xi }_2}}(t)}{dt}={\alpha }^{\overline{{\xi }_2}}{[T]}_{\chi }\left(1-R_{\chi }^{\overline{{\xi }_2}}(t)\right)-{\beta }^{\overline{{\xi }_2}}{R}_{\chi }^{\overline{{\xi }_2}}(t)& (14)\end{array}$

$\begin{array}{cc}\frac{d[X](t)}{dt}={\alpha }^{\overline{{\xi }_2}}{R}_{\chi }^{\overline{{\xi }_2}}(t)-{\beta }^{\overline{{\xi }_2}}[X](t)& (15)\end{array}$

$\begin{array}{cc}{r}_{\chi }^{\overline{{\xi }_2}}((t))=\frac{{[X]}^n(t)}{{[X]}^n(t)+K_d}& (16)\end{array}$

$\begin{array}{cc}{I}_{\chi }^{\overline{\xi }}(t)=g^{\overline{\xi }}{r}_{\chi }^{\overline{\xi }}(t)\left(V_P(t)-V^{\overline{\xi }}\right)& (17)\end{array}$

$\begin{array}{cc}{\kappa }_m\frac{d{V}_P(t)}{dt}=-{\displaystyle \sum _{\chi \in \{1,2\}}{I}_{\chi }^{\overline{\xi }}(t).C_{\chi }-I_P^{\lambda }(t)}& (18)\end{array}$

$\begin{array}{cc}{I}_P^{\lambda }(t)=g_P^{\lambda }\left(V_P(t)-V_P^{\lambda }\right)& (19)\end{array}$

Let V_χ, χ ∈ {1, 2} be the “ensemble” membrane potential of two pre-synaptic neuronal population that are afferent to the post-synaptic population P such that the synapses made by χ = 1 is mediated by a ligand-gated receptor $\overline{{\xi }_1}\in \{\text{AMPA},{\text{ GABA}}_A\}$ while that made by χ = 2 is mediated by a secondary-messenger-gated receptor $\overline{{\xi }_2}\in \left\{{\text{GABA}}_B\right\}$. The concentration of neurotransmitters [T]_χ in the synaptic cleft is defined as a function of V_χ and is expressed by a sigmoid function (Equation 12) where T_max is the maximum neuronal concentration in the synaptic cleft and is well approximated by 1 mM (Destexhe et al., 1998), θ_s represents the threshold at which [T] = 0.5T_max and σ_s denote the steepness of the sigmoid. The proportion of open ion-channels due to the bound receptors ${\overline{\xi }}_1$ on the ensemble membrane of the post-synaptic cell population corresponding to the synapse made by the population χ = 1 is defined in Equation (13) where ${\alpha }^{\overline{{\xi }_1}}$ and ${\beta }^{\overline{{\xi }_1}}$ are the forward and backward rate constants, respectively for transmitter binding. The transition diagram is shown in Figure 3A. However, GABA_B mediated synapses, unlike AMPA and GABA_A synapses, activate G-proteins which in turn act as the “secondary messengers” and initiate the opening of ion channels. The process is defined in Equations (14–16) where $R_{\chi }^{\overline{{\xi }_2}}$ is the fraction of activated $\overline{{\xi }_2}$ 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_{\chi }^{\overline{{\xi }_2}}$ is the fraction of open ion channels caused by binding of X with independent binding sites; ${\alpha }^{\overline{{\xi }_2}}$ and ${\beta }^{\overline{{\xi }_2}}$ are the binding rate constants; n is the number of bound receptor sites and K_d is the dissociation constant of binding of X with the ion channels. The transition diagram of this process is shown in Figure 3B. The resulting ensemble PSC mediated by the receptor $\overline{\xi }\equiv \overline{{\xi }_1}\cup \overline{{\xi }_2}$ due to a synapse from the pre-synaptic population χ is defined in Equation (17) where $g^{\overline{\xi }}$ and $V^{\overline{\xi }}$ are the maximum conductance and reverse potential, respectively corresponding to ξ mediated synapse. V_P (Equation 18) is the ensemble post-synaptic potential (PSP) of P, where κ_m is the ensemble membrane capacitance of P, C_χ, χ ∈ {1, 2} is the synaptic connectivity parameter. I^λ_P (Equation 19) is the ensemble leak current of the post-synaptic membrane, where g^λ_P and V^λ_P 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. In the following section, we implement this framework in a neural mass model of the thalamocortical circuitry.

FIGURE 3

Figure 3. The state transition diagrams for (A) AMPA and GABA_A neuro-receptor dynamics defined in Equation (13); α and β are rate of transitions between the two states and α is a function of the transmitter concentration in the synaptic cleft [T] defined in Equation (12); (B) GABA_B neuro-transmission as defined in Equations (14–16)—the neurotransmitter T binds to the inactivated receptor R₀; a fraction of activated receptors R act as a catalyst to transform the G-protein from an inactivated form X₀ to an activated form X, which binds at n independent sites to open a fraction of the ion channels. The desensitized state of the ion-channels are ignored in this work for brevity (see Destexhe et al., 1998 for a detailed comparison of kinetic models with more than two-states).

2.2. Neural Mass Model of a Thalamocortical Circuitry with Kinetic Synapses

The thalamocortical circuitry is shown in Figure 2C and consists of the two thalamic cell populations that communicate with the cortex viz. the TCR and TRN. The third group of cells viz. the Interneurons (IN) participate in intra-thalamic communications and are ignored here for brevity. The synaptic structure and connectivity are informed from experimental data based on the dorsal thalamic Lateral Geniculate Nucleus (LGNd) (Horn et al., 2000). The input to the model is assumed to be the ensemble membrane potential of pre-synaptic retinal cells (V_ret) in a resting state with no sensory input and is simulated using a Gaussian white noise (da Silva et al., 1973). The TCR cells make AMPA receptor mediated glutamatergic synapses on the TRN cells (other types of glutamatergic receptors are ignored in this work for brevity); the TRN cells make GABA-ergic synapses on the TCR cells mediated by both the ligand-gated GABA_A and the secondary-messenger-gated GABA_B receptors. Furthermore, the TRN cells make GABA_A receptor mediated synapses within the population. The model is defined in Equations (20–27); all variables and parameters in the model are assumed to be the ensemble representation corresponding to a neural mass:

$\begin{array}{cc}{[T]}_{\overline{\Psi }}(V_{\overline{\Psi }}(t))=\frac{T_{\text{max}}}{1+exp\left(-\frac{V_{\overline{\Psi }}(t)-{\theta }_s}{{\sigma }_s}\right)}& (20)\end{array}$

$\begin{array}{cc}\frac{d{r}_{\overline{\Psi }}^{{\overline{\eta }}_1}(t)}{dt}={\alpha }^{{\overline{\eta }}_1}{[T]}_{\overline{\Psi }}\left(1-r_{\overline{\Psi }}^{{\overline{\eta }}_1}(t)\right)-{\beta }^{{\overline{\eta }}_1}{r}_{\overline{\Psi }}^{{\overline{\eta }}_1}(t)& (21)\end{array}$

$\begin{array}{cc}\frac{d{R}_{\overline{\Psi }}^{\overline{{\eta }_2}}(t)}{dt}={\alpha }_1^{\overline{{\eta }_2}}{[T]}_{\overline{\Psi }}\left(1-R_{\overline{\Psi }}^{\overline{{\eta }_2}}(t)\right)-{\beta }_1^{\overline{{\eta }_2}}{R}_{\overline{\Psi }}^{\overline{{\eta }_2}}(t)& (22)\end{array}$

$\begin{array}{cc}\frac{d[X](t)}{dt}={\alpha }_2^{\overline{{\eta }_2}}{R}_{\overline{\Psi }}^{\overline{{\eta }_2}}(t)-{\beta }_2^{\overline{{\eta }_2}}[X](t)& (23)\end{array}$

$\begin{array}{cc}{r}_{\overline{\Psi }}^{\overline{{\eta }_2}}((t))=\frac{{[X]}^n(t)}{{[X]}^n(t)+K_d},& (24)\end{array}$

$\begin{array}{cc}{I}_{\overline{\Psi }}^{\overline{\eta }}(t)=g^{\overline{\eta }}{r}_{\overline{\Psi }}^{\overline{\eta }}(t)\left(V_{\overline{\Upsilon }}(t)-V^{\overline{\eta }}\right)& (25)\end{array}$

$\begin{array}{cc}{\kappa }_m\frac{d{V}_{\overline{\Upsilon }}(t)}{dt}=-{\displaystyle \sum _{\overline{\Psi }\in \{\text{ret},\text{ trn},\text{ trn}\}}{I}_{\overline{\Psi }}^{\overline{\eta }}(t).C_{\overline{u}\overline{v}\overline{w}}-I_{\overline{\Upsilon }}^{\lambda }(t)},& (26)\end{array}$

$\begin{array}{cc}{I}_{\overline{\Upsilon }}^{\lambda }(t)=g_{\overline{\Upsilon }}^{\lambda }(V_{\overline{\Upsilon }}(t)-V_{\overline{\Upsilon }}^{\lambda }),& (27)\end{array}$

where $\overline{\Psi }\in \{\text{ret},\text{ trn},\text{ trn}\}$ represent the afferent cell populations; $\overline{\Upsilon }=\{\text{trn},\text{ trn}\}$ represent the efferent cell populations; ${\overline{\eta }}_1\in \{\text{AMPA},{\text{ GABA}}_A\}$, ${\overline{\eta }}_2\in \left\{{\text{GABA}}_B\right\}$, $\overline{\eta }\equiv {\overline{\eta }}_1\cup {\overline{\eta }}_2$; $C_{\overline{u}\overline{v}\overline{w}}$ are connectivity parameters where $\overline{u}\in \{t,n\}$ and $\overline{v}\in \{r,t,n,s\}$ denote the post-synaptic and pre-synaptic cell populations, respectively of the retina (r), TCR (t), TRN (n), while s denote an intra-population afferent; $\overline{w}\in \{e,i\}$ represent an excitatory (e) or an inhibitory (i) synapse. All other parameter nomenclatures are as defined in section 2.1. The initial parameter values are mentioned in Table 1.

TABLE 1

Table 1. Initial values of the parameters defined in Equations (21–27).

3. Results

The ODEs are solved using the 4th/5th order Runge-Kutta-Fehlberg method (RKF45) in Matlab for a total duration of 600 s (10 min) at a resolution of 1 ms. The output voltage time series is averaged over 20 simulations, each simulation run with different seed for the noisy input. For frequency analysis, an epoch from 100–599 s of the output signal is sampled every 4 ms (250 Hz) and bandpass filtered between 3.5–14 Hz with a Butterworth filter of order 10. Short Time Fourier Transform (STFT) is done with a Hamming window of duration 10 s and overlap of 50%.

The model displays a point-attractor mode behavior (initial transient oscillations before settling down to a low amplitude noisy output, which reflects the noisy input of the model) corresponding to initial parameter values (Figure 4A). There is a behavioral transition in the model to a limit cycle mode with increasing values of β^ampa, which correlates with a decrease in the fraction of open ion channels in the post-synaptic ensemble membrane (Figures 4B,C). Varying α^ampa, on the other hand, does not affect the model behaviour (Figures 4D,E). A transition from the limit cycle mode to a spindling mode is effected in the model by increasing g^ampa, the post-synaptic membrane conductance for AMPA mediated synapses in both TCR and TRN cell population, and shown in Figures 4F,G. STFT of the output time series indicates the non-stationary behavior of the model (Figures 4H–K). A decrease and increase, respectively of the theta and alpha band components imply an overall increase in frequency with increasing values of g^ampa ≡ {g^ampa_TCR, g^ampa_TRN}, where g^ampa_TCR and g^ampa_TRN correspond to the incoming signal from the retina (to the TCR) and TCR (to the TRN), respectively in the model. These observations are consistent with similar reports of a transition in the state of the model output with increasing values of g^ampa in Golomb et al. (1996; pp. 756–757), accompanied by an abrupt increase in the ratio of the frequency of oscillation of the TCR and the TRN cell populations; we have not studied the latter aspect in this work. A more detailed study on the model presented herein where g^ampa_TCR and g^ampa_TRN are varied separately specify the g^ampa_TCR as the control parameter that causes a bifurcation in the model output from a limit cycle mode to the spindling mode with an increase in its value. On the other hand, the g^ampa_TRN does not effect any behavioral change in the model output, rather, it is effective in increasing the inter-spindle frequency with an increase in its value when the model is in a spindling mode. This observation implies that a change in AMPA receptor related attributes in the TRN plays a role in modulating thalamocortical spindle oscillations, which finds strong support in the experimental study by von Krosigk et al. (1993), where “activation of AMPA-kainate receptors on the PGN” (Peri-geniculate nucleus—the part of the TRN associated with the LGNd) is described as “critical to the generation of spindle waves”. Furthermore, this observation is in line with the TRN being widely implicated as being the key “ingredient” in the generation of thalamocortical spindle oscillations (McCormick, 1992; Steriade et al., 1993).

FIGURE 4

Figure 4. The model output time series with (A) all parameters at their initial values as in Table 1. (B) The model displays a bifurcation in output behavior when β^ampa is increased to 0.2; (C) a further increase in the parameter shows sustained oscillations with increased magnitude and decreased frequency. The frequency behavior may be observed qualitatively in the embedded line plots displaying the respective time-series in each plot for an arbitrarily selected period of 60 s. (D) An increase in α^ampa has a reverse effect resulting in reduced magnitude and increased frequency in the limit-cycle mode; (E) further increase in the parameter do not show any significant effect on the magnitude or the limit-cycle behavior while there is a slight increase in frequency. (F) Maintaining these modified values of α^ampa = 20 and β^ampa = 1, an increase in g^ampa brings a bifurcation in model behavior from a limit-cycle mode to a “spindling” mode. A “zoomed-in” plot from the 3rd minute to the 9th minute is shown; the initial transient oscillations are neglected. (G) The frequency of spindle oscillations increase with increasing values of g^ampa. This is indicated by a distinct (H,I) decrease (more blue pixels) in theta band components and (J, K) increase (more red, orange and yellow pixels) in alpha band components in the corresponding output time series plots. The abscissa in the figures denote (A–G) time (seconds) (H–K) time windows (seconds).

Varying the GABA-ergic synaptic attributes when the model is in a point-attractor mode does not show any change in model behavior. When the model is in a spindling mode (Figure 5A), increased synchronization within the limit cycle mode with increasing values of ggabaATRN to TCR (Figures 5B,C) is observed. An increase in the parameter βgabaA affects the output only when the model is in a limit-cycle mode and counters the effect of increase in ggabaATRN to TCR (Figure 5D). However, varying αgabaA does not affect the model output. For ggabaATRN to TCR ≥slant 0.5, which is the approximate bifurcation point (Figure 5B), increasing ggabaATRN to TRN causes the model to revert back to the spindling mode; the frequency of the inter-spindle oscillations increase with increasing values of the parameter (Figures 5E,F). This is also indicated by a decrease (Figures 5H,I) and increase (Figures 5J,K) of theta and alpha band components respectively in the STFT of the output time series. In other words, decreasing values of the parameter ggabaATRN to TRN causes increased synchronization within the spindling mode behavior of the model along with a decrease in the inter-spindle frequency. However, blocking ggabaATRN to TRN effects a switch in the model behavior to a very low-frequency oscillatory state (Figure 5G). These results are consistent with experimental findings (von Krosigk et al., 1993) where application of GABAA inhibitor either “abolished spindle waves or decreased within-spindle frequency,” which correspond to the condition of either blocking or decreasing, respectively of ggabaATRN to TRN in our model. Thus, the model implicate the intra-TRN synaptic activity to be a key factor in sustaining spindle oscillations in the thalamocortical circuitry, an observation which conforms to those made in Golomb et al. (1996; p. 755). Furthermore, a “frequency jump” with increasing ggabaA, and associated transition in model behavior is also reported in Golomb et al. (1996; see Figure 7) as a comparative study between the TCR and TRN cells. This is similar with the increase in frequency of the spindle oscillations corresponding to increasing ggabaATRN to TRN in the present model, although we have not done a comparative study with the TRN cell population behavior. However, the current study implicate the increased post-synaptic conductance for GABAA receptors in the TCR cell population (ggabaATRN to TCR) to play a significant role in effecting state-transition between spindle and slow-wave oscillations, an observation that is yet to find support from experimental or model-based studies.

FIGURE 5

Figure 5. (A) The model output (corresponding to the “zoomed-in” plot in Figure 4G) when αampa = 20, βampa = 1, gampa = 0.3). Retaining these parameter values, (B) increasing ggabaAn2r, where “n2r” denotes the TRN to TCR pathway, from its initial value effects a bifurcation in model behavior from a spindling mode to a limit-cycle mode, indicating highly synchronized oscillations in the thalamocortical circuitry. (C) Synchronization increases with increase in ggabaAn2r, indicated by increased magnitude and decreased frequency of oscillation. For ggabaAn2r ≥slant 0.5 [approximate point of bifurcation, shown in (B)], (D) increasing βgabaA effects an increase in frequency within the limit-cycle mode, while (E) increasing ggabaAn2n, where “n2n” denotes the self inhibitory pathway of the TRN, effects a transition from the limit-cycle mode to the spindling mode. (F) Further increase in ggabaAn2n causes an increase in the frequency of the spindle oscillations, indicated by (H,I) a decrease (more blue pixels) in theta band components and (J,K) increase (more red, orange and yellow pixels) in alpha band components in the corresponding output time series plots. (G) Blocking of ggabaAn2n shows a very low-frequency (≈0.03 Hz) synchronized oscillation whose magnitude decreases with time. The abscissa in the figures denote (A–G) time (seconds) (H–K) time windows (seconds).

A quiescent state is observed corresponding to blocking either AMPA (gampa = 0) or both GABAA (ggabaATRN to TCR = 0) and GABAB (ggabaBTRN to TCR = 0) mediated synapses in the TRN to TCR pathway (not shown). This is consistent with both experimental (von Krosigk et al., 1993) and model-based (Golomb et al., 1996) studies. The role of the synaptic parameters in the GABAB pathway in our model was minimal—to sustain a non-quiescent model behavior with blockage of GABAA; to sustain a high amplitude of limit-cycle oscillations in the model. Again, this is in agreement with experimental studies (von Krosigk et al., 1993), where activation of GABAB receptors are reported as “not essential” for generating synchronized oscillations, while application of GABAB antagonist abolished “evoked or spontaneous slowed oscillations.” The model in Golomb et al. (1996; see in Discussion p. 763) is also mentioned as being consistent with these experimental results.

In a recent work (Bhattacharya et al., 2012), a simple neural mass model implementing kinetic modeling for synaptic transmission is presented; the synaptic connectivity parameters in the model correlate directly to that of an alpha function based neural mass model [modified Alpha Rhythm model (modARm) from Bhattacharya et al. (2011a)]. The model behavior is studied corresponding to changes in the synaptic connectivity parameters as well as transmitter concentration related parameters, and a relevant comparison is made with the modARm. However, the model presented in this work has a larger set of synaptic connectivity parameters; model behavior corresponding to this parameter space and its usefulness in understanding neurological disorders will be the topic of a future work.

4. Discussion

The work presented here explores a novel approach toward correlating current neural mass model based studies with underlying cellular mechanisms during synaptic transmission. The aim is to underpin the synaptic correlates of abnormal brain oscillations in neurological and psychiatric disorders such as observed in Electroencephalogram (EEG). A kinetic framework for modeling AMPA and GABA receptor mediated synapses is implemented in an existing thalamocortical neural mass model consisting of an excitatory and an inhibitory neural mass, representing cell populations of the thalamocortical relay (TCR) and the thalamic reticular nucleus (TRN), respectively. Parameters in the model are assumed to be “ensemble” representations of the corresponding attributes in a single neuron. A preliminary observation is made on the model behavior by varying the parameters corresponding to the post-synaptic membrane conductance of the cell populations as well as the forward and reverse rates of synaptic reaction; of specific interest is the transition of the model behavior between the spindle oscillatory mode and the limit-cycle mode, the latter resembling the slow-wave (high-amplitude, low-frequency) synchronized oscillations that are signatures of absence seizures as well as slow-wave sleep. Furthermore, only the alpha (8–13 Hz) and theta (4–7 Hz) frequency bands of the output power spectra are studied here, as EEG alpha and theta bands are believed to have a strong correlation with thalamocortical oscillations (Hughes et al., 2004).

The results indicate that: (1) The post synaptic membrane conductance for both AMPA and GABA_A receptors in the TRN cell population play a role in sustaining spindle oscillations of the TCR cell population (the model output). (2) Blocking the GABA_A mediated synapses in the self-inhibitory feedback pathway of the TRN cell population effects synchronized oscillations with high amplitude and increased time-period of oscillation (≈0.03 Hz). (3) The post-synaptic membrane conductance for GABA_B in the TCR cell population does not play any role in generating or sustaining spindle oscillations, but is responsible for sustaining the slow-wave oscillations in the model associated with blocking of the intra-TRN GABA_A synapses. (4) Blocking both GABA_A and GABA_B or only the AMPA mediated synapses in the TCR cell population results in a quiescent model output. These findings are consistent with in vitro studies based on multiple unit recordings from ferret thalamic slices (von Krosigk et al., 1993) as well as single-neuron-level model based studies (Golomb et al., 1996). In addition, this study identifies—(a) the reverse rate of transmitter binding as an important attribute in effecting thalamocortical synchronized oscillations that can be induced in the model by increasing (decreasing) the fraction of open channels due to GABA_A (AMPA) mediated synapses in the TCR (TRN) cell populations; (b) the post-synaptic membrane conductance for GABA_A in the TCR cell population as a control parameter for effecting a behavioral transition in the model.

It may be noted that the above-mentioned observations are only a qualitative comparison with single-neuron-level model-based (Golomb et al., 1996) and experimental (von Krosigk et al., 1993) studies; a drawback of the current work is a lack of quantitative comparison with these studies. The neural mass model presented in this work is at a mesoscopic scale, representation of a population of ≈10⁴−10⁷ neurons, unlike that in Golomb et al. (1996), which is at single-neuronal-level. Similarly, the multiple unit recording based study in von Krosigk et al. (1993) observes neuronal behavior of either a single neuron or a population of <10² neurons. In addition, the modeling and simulation methods in the current work and that in Golomb et al. (1996) are not similar. Thus, a quantitative comparison of the current work with these studies may lead to erroneous conclusions. However, model validation with experimental data is a crucial criteria when investigating brain disorders. Along these lines, an ongoing work is investigating ways to validate the model presented herein with EEG data, and will be the topic of a future study.

The model structure in the current work is a considerably simplified representation of the thalamocortical circuitry. The role of the thalamocortical circuitry in generating slow wave brain oscillations is discussed at length in Steriade et al. (1993), based on in vivo and in vitro studies. More recently, three parameters in the thalamo-cortico-thalamic loop viz. the cortico-thalamic, thalamo-cortical and intra-thalamic pathways are specified in Breakspear et al. (2006) for generating instabilities in the thalamocortical circuitry, leading to synchronized oscillations such as seen during absence seizures. Furthermore, a non-linear dynamical analysis of the model is shown to predict seizure onset by validating with patient EEG data. In a previous research (Bhattacharya et al., 2011b), we have proposed a more elaborate alpha-function based neural mass model that have considered these vital pathways in the thalamocorticothalamic loop. Also, we have performed a non-linear dynamical analysis of a simple thalamocortical model based on alpha functions in Bhattacharya et al. (2013) to understand EEG power spectra abnormalities associated with several neurological disorders. Such research directions will be considered as an extended work based on the model presented herein.

It is worth mentioning here that biologically plausible parameterizations has been a major constraint in neural mass modeling of brain dynamics. This is largely due to insufficient experimental data, published or otherwise, as well as to a lack of “homogeneity” of published data from different experimental laboratories. The trend thus far has been to use biologically plausible data if and when available; otherwise, i.e., for parameter values that cannot be availed from experimental data, the models are tuned to estimated parameter values which provide a desirable output in context to the objectives of the research [the reader may refer to Robinson et al. (2004) for a model parameterizations related work and discussion]. Thus, the model in Breakspear et al. (2006) was based on neurophysiological parameters obtained from Robinson et al. (2002), which in turn are based on inverse parameterizations during model validation with EEG data from patients of epileptic seizures. The parameterizations of the model presented in this work is largely based on neurophysiological parameters obtained from experimental studies: the cellular-level parameters, including those of the synaptic kinetics, are based on in vitro studies and model-based studies of thalamocortical tissue by von Krosigk et al. (1993) and Golomb et al. (1996), respectively; the model connectivity parameters are based on experimental studies of the cat and rat thalamus obtained from Horn et al. (2000); Sherman and Guillery (2001). On the other hand, the extrinsic (retinal) input and neuro-transmitter concentration parameter values are adjusted to maintain a “dynamically active” model behavior (this is as opposed to a continuous “quiescent” state of the model corresponding to certain parameter values, and does not conform to biology). However, technological advances in the field of neuro-imaging during the last decade such as functional Magnetic Resonance Imaging (fMRI), Diffusion Tensor Imaging (DTI) and Transcanial Magnetic Stimulation (TMS) are paving the way for biologically-realistic mapping of parameter values in computational models; for example as in Izhikevich and Edelman (2008).

The observations made herein support the motivation toward this preliminary work, which is to correlate higher-level brain dynamics with underlying cellular-level synaptic mechanisms. It may be noted that in all our previous works using alpha function based neural mass models, the emphasis has been on studying the model behavior with varying values of synaptic connectivity parameters toward a meaningful mapping to Alzheimer disease-related EEG anomalies. However, such “synaptic parameter variation only” studies are highly constrained and do not make much sense when trying to understand generic brain-state conditions e.g., the sleep-awake cycle, or several other neurological and psychiatric disorders e.g., absence seizures, which rely heavily on various aspects of cellular dynamics in the thalamocortical circuitry. Rather, the emphasis of this work is on laying the ground-work for a more elaborate, and yet computationally efficient scheme, whereby large-scale computational models may be simulated to mimic brain rhythms, which can then be correlated to model parameters emulating cellular dynamics. The synaptic transmission kinetics and subsequent post-synaptic membrane parameters are some of the key constituents of brain signaling, and are affected significantly in various brain diseases. Clearly, the alpha-function based neural mass models are inadequate in dealing with research directions where model parameters can be mapped in a biologically plausible manner to synaptic attributes. In terms of computational efficiency, the time for simulating 20 trials with the model presented in this work takes 60 s; this may be contrasted with 600 s for simulating a similar model [the modified Alpha Rhythm model in Bhattacharya et al. (2011a)] based on alpha functions. This is a dramatic improvement in computational efficiency and highlight the plausibility of using the kinetic-model based neural mass modeling framework in simulating large-scale computational models toward mimicking real-time EEG signals. This in turn will provide a powerful tool for specifying cellular pathways that need be targeted for symptomatic alleviation of anomalous brain rhythms as well as to inform effective neuropharmacological research directions.

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.

References

Aradi, I., and Erdi, P. (2006). Computational neuropharmacology: dynamical approaches in drug discovery. Trends Pharmacol. Sci. 27, 240–243. doi: 10.1016/j.tips.2006.03.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Basar, E., and Guntekin, B. (2008). A review of brain oscillations in cognitive disorders and the role of neurotransmitters. Brain Res. Rev. 1235, 172–193. doi: 10.1016/j.brainres.2008.06.103

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bernard, C., Ge, Y., Stockley, E., Willis, J., and Wheal, H. (1994). Synaptic integration of NMDA and non-NMDA receptors in large neuronal network models solved by means of differential equations. Biol. Cybern. 70, 267–273. doi: 10.1007/BF00197607

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bhattacharya, B. S., Cakir, Y., Serap-Sengor, N., Maguire, L. P., and Coyle, D. (2013). Model-based bifurcation and power spectral analyses of thalamocortical alpha rhythm slowing in Alzheimer's disease. Neurocomputing 115, 11–22. doi: 10.1016/j.neucom.2012.10.023

CrossRef Full Text

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, et al. (New York, NY: 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 | Pubmed Full Text | CrossRef Full Text

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 E. P. V. Alessandro, D. Wlodzislaw, E. Peter, M. Francesco, and P. Gunther (Berlin; Heidelberg: Springer Verlag), 645–652.

Breakspear, M., Roberts, J., Terry, J., Rodrigues, S., Mahant, N., and Robinson, P. A. (2006). A unifying explanation of primary generalized seizures through nonlinear brain modelling and bifurcation analysis. Cereb. Cortex 16, 1296–1313. doi: 10.1093/cercor/bhj072

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Breakspear, M., Terry, J. R., and Friston, K. J. (2003). Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a nonlinear model of neuronal dynamics. Neurocomputing 52, 151–158. doi: 10.1016/S0925-2312(02)00740-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

da Silva, F. H. L., van Lierop, T. H. M. T., Schrijer, C. F., and van Leeuwen, W. S. (1973). Essential differences between alpha rhythms and barbiturate spindles: spectra and thalamo-cortical coherences. Electroencephalogr. Clin. Neurophysiol. 35, 641–645. doi: 10.1016/0013-4694(73)90217-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

de Haan, W., Mott, K., van Straaten, E. C. W., Scheltens, P., and Stam, C. J. (2012). Activity dependent degeneration explains hub vulnerability in Alzheimer's disease. PLoS Comput. Biol. 8:e1002582. doi: 10.1371/journal.pcbi.1002582

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4:e1000092. doi: 10.1371/journal.pcbi.1000092

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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. doi: 10.1016/B978-012373649-9.50016-8

CrossRef Full Text

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 | Pubmed Full Text

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

Destexhe, A., Mainen, Z., 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.

Destexhe, A., Mainen, Z., and Sejnowski, T. (2002). “Kinetic models for synaptic interactions,” in The Handbook of Brain Theory and Neural Networks, ed M. Arbib (Cambridge, MA: MIT Press), 1126–1130.

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

Golomb, D., Wang, X.-J., and Rinzel, J. (1994). Synchronization properties of spindle oscillations in a thalamic reticular nucleus model. J. Neurophysiol. 72, 1109–1126.

Pubmed Abstract | Pubmed Full Text

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 | Pubmed Full Text

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 | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

Izhikevich, E. M., and Edelman, G. M. (2008). Large-scale model of mammalian thalamocortical systems. Proc. Natl. Acad. Sci. U.S.A. 105, 3593–3598. doi: 10.1073/pnas.0712231105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liljenström, H. (2012). Mesoscopic brain dynamics. Scholarpedia 7:4601. doi: 10.4249/scholarpedia.4601

CrossRef Full Text

McCormick, D. A. (1992). Neurotransmitter actions in the thalamus and cerebral cortex and their role in neuromodulation of thalamocortical activity. Progr. Neurobiol. 39, 337–388. doi: 10.1016/0301-0082(92)90012-4

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pons, A. J., Cantero, J. L., Atienza, M., and Garcia-Ojalvo, J. (2010). Relating structural and functional anomalous connectivity in the ageing brain via neural mass modelling. Neuroimage 52, 848–861. doi: 10.1016/j.neuroimage.2009.12.105

Pubmed Abstract | Pubmed Full Text | 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 | Pubmed Full Text

Robinson, P., Rennie, C., and Rowe, D. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65:041924. doi: 10.1103/PhysRevE.65.041924

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Phillips, A. J. K., Fulcher, B. D., Puckeridge, M., and Roberts, J. A. (2011). Quantitative modelling of sleep dynamics. Philos. Trans. R. Soc. A 369, 3840–3854. doi: 10.1098/rsta.2011.0120

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Sotero, R. C., Tujillo-Barreto, N. J., and Iturria-Medina, Y. (2007). Realistically coupled neural mass models can generate EEG rhythms. Neural Comput. 19, 479–512. doi: 10.1162/neco.2007.19.2.478

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stam, C. J., Pijn, J., Suffczyński, P., and da Silva, F. (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 | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

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ńnski, P., Kalitzin, S., and Silva, F. H. 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 | Pubmed Full Text | CrossRef Full Text

Ursino, M., Cona, F., and Zavaglia, M. (2010). The generation of rhythms within a cortical region: analysis of a neural mass model. Neuroimage 52, 1080–1094. doi: 10.1016/j.neuroimage.2009.12.084

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

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. 92, 5577–5581. doi: 10.1073/pnas.92.12.5577

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, X.-J., and Rinzel, J. (1992). Alternating and synchronous rhythms in reciprocally inhibitory model neurons. Neural Comput. 4, 84–97. doi: 10.1162/neco.1992.4.1.84

CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

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 | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text