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

A novel direction to existing neural mass modeling technique is proposed where the commonly used “alpha function” for representing synaptic transmission is replaced by a kinetic framework of neurotransmitter and receptor dynamics. The aim is to underpin neuro-transmission dynamics associated with abnormal brain rhythms commonly observed in neurological and psychiatric disorders. An existing thalamocortical neural mass model is modified by using the kinetic framework for modeling synaptic transmission mediated by glutamatergic and GABA (gamma-aminobutyric-acid)-ergic receptors. The model output is compared qualitatively with existing literature on in vitro experimental studies of ferret thalamic slices, as well as on single-neuron-level model based studies of neuro-receptor and transmitter dynamics in the thalamocortical tissue. The results are consistent with these studies: the activation of ligand-gated GABA receptors is essential for generation of spindle waves in the model, while blocking this pathway leads to low-frequency synchronized oscillations such as observed in slow-wave sleep; the frequency of spindle oscillations increase with increased levels of post-synaptic membrane conductance for AMPA (alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic-acid) receptors, and blocking this pathway effects a quiescent model output. In terms of computational efficiency, the simulation time is improved by a factor of 10 compared to a similar neural mass model based on alpha functions. This implies a dramatic improvement in computational resources for large-scale network simulation using this model. Thus, the model provides a platform for correlating high-level brain oscillatory activity with low-level synaptic attributes, and makes a significant contribution toward advancements in current neural mass modeling paradigm as a potential computational tool to better the understanding of brain oscillations in sickness and in health.


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 4 -10 7 ), 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)(9)(10)(11)(12)(13). 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 neurotransmission 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. (1994Golomb et al. ( , 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") ionchannels (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(Destexhe et al., , 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-4isoxazolepropionic-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 postsynaptic 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.

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): wherew ∈ {e, i} represents pre-synaptic neuronal populations which make excitatory (e) and inhibitory (i) synapses on a post-synaptic neuronal population; τw is the time constant and FIGURE 1 | Block diagram of a single "neuronal mass" in current state-of-the-art neural mass models.
Hw is the amplitude of the synapse; E N w (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 P w is the "ensemble firing rate" of P and is defined by a sigmoid function where 2e 0 is the maximum firing rate of the population, s 0 is the threshold potential at which the neurons spike and ν is the sigmoid steepness parameter.

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-messengergated GABA-ergic receptors GABA A and GABA B , respectively. The input E N ξ (t),ξ ∈ {AMPA, GABA A , 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): where hξ(t) is the synaptic transmission function with τ ā ξ and τ b ξ as the rise and decay times, respectively; g N ξ denote the postsynaptic "ensemble" membrane conductance; Vξ is the reversal potential for the synapse mediated byξ; 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 P ξ (t) is as defined in Equation (5) and is the pre-synaptic firing rate input to other neuronal populations.

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): Let V χ , χ ∈ {1, 2} be the "ensemble" membrane potential of two pre-synaptic neuronal population that are afferent to the postsynaptic population P such that the synapses made by χ = 1 is mediated by a ligand-gated receptorξ 1 ∈ {AMPA, GABA A } while that made by χ = 2 is mediated by a secondary-messenger-gated receptorξ 2 ∈ {GABA B }. The concentration of neurotransmitters Equations (14-16)-the neurotransmitter T binds to the inactivated receptor R 0 ; a fraction of activated receptors R act as a catalyst to transform the G-protein from an inactivated form X 0 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).
[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ξ 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 αξ 1 and βξ 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ξ 2 χ is the fraction of activatedξ 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ξ 2 χ is the fraction of open ion channels caused by binding of X with independent binding sites; αξ 2 and βξ 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ξ ≡ξ 1 ∪ξ 2 due to a synapse from the pre-synaptic population χ is defined in Equation (17) where gξ and Vξ 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.

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: where¯ ∈ {ret, trn, trn} represent the afferent cell populations;Ῡ = {trn, trn} represent the efferent cell populations;η 1 ∈ {AMPA, GABA A },η 2 ∈ {GABA B },η ≡η 1 ∪η 2 ; Cūvw are connectivity parameters whereū ∈ {t, n} andv ∈ {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 intrapopulation afferent;w ∈ {e, i} represent an excitatory (e) or an inhibitory ( Golomb et al. (1996) and Suffczyński et al. (2004).
In Equation (20), both θ s and σ s act as bifurcation parameters in the model (see Bhattacharya et al., 2012). However, the emphasis here is on post-synaptic membrane attributes as in von Krosigk et al. (1993) and Golomb et al. (1996).
Thus these parameters (θ s = −35 and σ s = 2) are set by trial and error at values just before the model undergoes bifurcation from a "point-attractor" mode to a "limit-cycle" mode, based on a recent study where we observed rich model dynamics and power spectral behavior around the bifurcation point (Bhattacharya et al., 2013); T max = 1 mM (Destexhe et al., 1994 Bhattacharya et al. (2011b) for details], when the model output showed an increased frequency content within the theta (4-7 Hz) and alpha (8-13 Hz) bands. C tre and C nte are as in Bhattacharya et al. (2011b). All variables in the ODEs are initialized to an arbitrarily small value 0.0002.

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 AMPAkainate 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).
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 g gaba A TRN to TCR (Figures 5B,C) is observed. An increase in the parameter β gaba A affects the output only when the model is in a limit-cycle mode and counters the effect of increase in g gaba A TRN to TCR (Figure 5D)  does not affect the model output. For g gaba A TRN to TCR 0.5, which is the approximate bifurcation point (Figure 5B), increasing g gaba A TRN 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 g gaba A TRN to TRN causes increased synchronization within the spindling mode behavior of the model along with a decrease in the inter-spindle frequency. However, blocking g gaba A TRN 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 GABA A inhibitor either "abolished spindle waves or decreased within-spindle frequency," which correspond to the condition of either blocking or decreasing, respectively of g gaba A TRN 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 g gaba A , 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 g gaba A TRN 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 GABA A receptors in the TCR cell population (g gaba A TRN 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.
A quiescent state is observed corresponding to blocking either AMPA (g ampa = 0) or both GABA A (g gaba A TRN to TCR = 0) and GABA B (g gaba B TRN 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 GABA B pathway in our model was minimal-to sustain a non-quiescent model behavior with blockage of GABA A ; to sustain a high amplitude of limit-cycle oscillations in the model. Again, this

FIGURE 5 | (A)
The model output (corresponding to the "zoomed-in" plot in Figure 4G) when α ampa = 20, β ampa = 1, g ampa = 0.3). Retaining these parameter values, (B) increasing g gaba A n2r , 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. is in agreement with experimental studies (von Krosigk et al., 1993), where activation of GABA B receptors are reported as "not essential" for generating synchronized oscillations, while application of GABA B 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.

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 slowwave 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 modelbased (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 4 −10 7 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 2 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.