Effects of Cholinergic Neuromodulation on Thalamocortical Rhythms During NREM Sleep: A Model Study

It has been suggested that cholinergic neurons shape the oscillatory activity of the thalamocortical (TC) network in behavioral and electrophysiological experiments. However, theoretical modeling demonstrating how cholinergic neuromodulation of thalamocortical rhythms during non-rapid eye movement (NREM) sleep might occur has been lacking. In this paper, we first develop a novel computational model (TC-ACH) by incorporating a cholinergic neuron population (CH) into the classical thalamo-cortical circuitry, where connections between populations are modeled in accordance with existing knowledge. The neurotransmitter acetylcholine (ACH) released by neurons in CH, which is able to change the discharge activity of thalamocortical neurons, is the primary focus of our work. Simulation results with our TC-ACH model reveal that the cholinergic projection activity is a key factor in modulating oscillation patterns in three ways: (1) transitions between different patterns of thalamocortical oscillations are dramatically modulated through diverse projection pathways; (2) the model expresses a stable spindle oscillation state with certain parameter settings for the cholinergic projection from CH to thalamus, and more spindles appear when the strength of cholinergic input from CH to thalamocortical neurons increases; (3) the duration of oscillation patterns during NREM sleep including K-complexes, spindles, and slow oscillations is longer when cholinergic input from CH to thalamocortical neurons becomes stronger. Our modeling results provide insights into the mechanisms by which the sleep state is controlled, and provide a theoretical basis for future experimental and clinical studies.

1. Introduction

Sleep plays a pivotal role in mental and physical health. During sleep, the brain alternates between two stages, rapid-eye movement (REM) and NREM. This alternation is reflected in the electrical rhythms generated by the thalamocortical system, evident in the electroencephalogram (EEG). Specifically, the thalamocortical rhythms in NREM sleep include K-complexes and slow oscillations (dominated by low-frequency [0.5, 2]Hz, high amplitude oscillations), as well as spindles (characterized by a waxing and waning waveform in the range [11, 16]Hz). In contrast, REM sleep exhibits low amplitude activity of higher frequency, which is similar with rhythms in wakefulness, including the alpha oscillation (Rasch and Born, 2013). Studies reveal that spindles and slow oscillations are very helpful in protecting and stabilizing sleep (Roth, 2009; Dang-Vu et al., 2010; Kim et al., 2012).

A number of results in behavioral and electrophysiological experiments suggest that cholinergic neurons, which are distributed in the basal forebrain and brainstem, help shape the oscillatory activity of the thalamocortical network (Bellingham and Funk, 2000; Kobayashi et al., 2003; Boutrel and Koob, 2004; Brown et al., 2012). Meanwhile, it has also been shown that thalamocortical processing is subject to the action of modulatory of acetylcholine (ACH) released by cholinergic neurons (McCormick, 1992). However, a theoretical framework for understanding the mechanisms of cholinergic neuromodulation on thalamocortical rhythms during NREM sleep is still lacking.

Neural mass models (NMMs), which describe the dynamics of neural populations (Wilson and Cowan, 1973; Lopes da Silva et al., 1974), have shown success in clarifying the mechanisms involved in the generation and regulation of oscillation patterns of the “wakefulness-sleep cycle.” For example, Robinson et al. (2002) developed a continuum model of large-scale brain electrical activity to understand the transition from the resting background state to the spike-wave state. Suffczynski et al. (2004) proposed a thalamocortical NMM to explain the relation between mechanisms that generate sleep spindles and those for that generate spike-wave activity. More recently, Weigenand et al. (2014) extended a cortical NMM to explore the mechanisms for generating K-complexes and slow oscillations. Cona et al. (2014) presented a new NMM to describe the sleeping thalamocortical system, where thalamic neurons exhibit two firing modes: bursting and tonic.

Although prior (NMM-based) studies have made progress in explaining sleep rhythms, the role of ACH modulation of thalamocortical activities, especially during sleep has not been thoroughly explored. In this work we model the details of cholinergic neuromodulation in the cortex and thalamus, and study the effects of cholinergic projections on thalamocortical rhythms during NREM sleep.

In this study, a new computational model (TC-ACH) is first proposed by integrating a cholinergic neuron population (CH) into classical thalamo-cortical circuitry, with connections between populations in TC-ACH constructed on the basis of the existing mechanisms. We concentrate on the neurotransmitter ACH released by cholinergic neurons, which impacts firing activity of thalamocortical neurons. The proposed model TC-ACH is used to investigate effects of ACH modulation on thalamocortical rhythms during NREM sleep in two ways: (1) measurement of effects of cholinergic modulation in the thalamus; (2) measurement of effects of cholinergic modulation in the thalamocortical system.

The rest of the paper is organized as follows. In section 2, the model framework and mathematical expression of TC-ACH are introduced. Section 3 presents the simulation results. A summary and discussion of our work are given in the last section.

2. Methods

2.1. Model Framework

This subsection introduces the topological structure and computational operators of TC-ACH systematically.

2.1.1. Topological Structure

Our computational model TC-ACH is composed of classical thalamo-cortical circuitry as well as a cholinergic neuron population, whose topological structure is illustrated in Figure 1.

FIGURE 1

Figure 1. Schematic representation of model TC-ACH.

The thalamo-cortical circuitry consists of two mutually interconnected modules: a cortical module and a thalamic module. Each module comprises two neural populations, where pyramidal cells (PY) and inhibitory interneurons (IN) are included in the cortical module; the thalamic module also includes thalamocortical cells (TC) and thalamic reticular cells (RE). Regarding connections between modules, excitatory neurons (PY or TC) in each module connect to both populations in the other module. Within each module, two populations project to each other, and there exist self-connections within PY, IN, and RE (Costa et al., 2016). Note that in the thalamo-cortical circuitry, the excitatory and inhibitory projections are mediated by AMPA and GABA receptors, respectively.

The cholinergic neuron population (CH) is located in the basal forebrain and brainstem. The neurons in CH release neurotransmitter (ACH) into the synaptic cleft, which changes the activity of thalamocortical neurons. In our proposed model, there exist three modulated cholinergic projections from CH to PY, RE, as well as the connection between TC and PY (represented by red dashed lines in Figure 1), and one mediated cholinergic projection from CH to RE (represented by blue the solid line in Figure 1). The coupling coefficients of these projections are denoted by K_PY, $K_{RE}^{Nrr}$, $K_{RE}^{Ca}$, K_TP, and $K_{RE}^I$, respectively. Here, all projections involve by ACH.

2.1.2. Computational Operators

The dynamical evolution in each population is implemented in two computational blocks.

The first computational block transforms the average membrane potential V(t) into the average firing rate Q(t), formulated by a sigmoid function with the form (Jansen et al., 1993)

$\begin{array}{ll}Q(t)=\frac{Q^{max}}{1+e^{-(V(t)-\theta )/\sigma }},& (1)\end{array}$

where Q^max, θ, σ represent the maximal firing rate, the firing threshold and the neural gain, respectively.

In the second computational block, the firing rate Q(t) is first transformed into the fraction of open channels r_ξ(t) by a convolution with an alpha function h_ξ(t), that is,

$\begin{array}{ll}{r}_{\xi }(t)=h_{\xi }(t)\otimes (N·Q(t)),& (2)\end{array}$

$\begin{array}{ll}{h}_{\xi }(t)={\gamma }_{\xi }^2·t{e}^{-{\gamma }_{\xi }t},t\ge 0.& (3)\end{array}$

Then the membrane potential V(t) is solved by

$\begin{array}{lll}\tau \dot{V}(t)& =& -I_L(t)-I_{AMPA}(t)-I_{GABA}(t)\\ & =& -{\overline{g}}_L·(V(t)-E_L)\\ & & -{\displaystyle \sum _{\xi }}{\overline{g}}_{\xi }{r}_{\xi }·(V(t)-E_{\xi }),& (4)\end{array}$

which is similar to the classical conductance-based form of Hodgkin and Huxley (1952) with one leak and two synaptic currents (say, ξ ∈ {AMPA, GABA}) (Weigenand et al., 2014). Moreover, for the convenience of calculating the convolution ⊗, r_ξ(t) can be equivalently obtained by solving the following second-order differential equation

$\begin{array}{ll}{\ddot{r}}_{\xi }(t)={\gamma }_{\xi }^2·(N·Q(t)-r_{\xi }(t))-2{\gamma }_{\xi }{\dot{r}}_{\xi }(t).& (5)\end{array}$

In (2)-(5), E_ξ and E_L denote the reversal potential of the synaptic current and leak current, respectively, ${\overline{g}}_{\xi }$ represents the synaptic input rate that scales r_ξ and E_ξ, g_L is the maximal conductivity of the leak current conductance, τ is the membrane time constant, N stands for the connectivity constant, and γ_ξ is the rate constant of synaptic response. The detailed dynamical evolution in each population is illustrated in Figure 2.

FIGURE 2

Figure 2. The dynamical evolution in each population in model TC-ACH.

2.2. Mathematical Expression

This subsection formulates the proposed model TC-ACH mathematically in the logic of five populations included in the model.

2.2.1. Cholinergic Neuron Population (CH)

We first consider the cholinergic neuronal population, which is an essential part of TC-ACH in studying the effects of cholinergic modulation on thalamocortical rhythms during NREM sleep. Due to the special role of ACH, here, we apply the concentration of ACH ([ACH]) rather than the firing rate applied in Equation (1) to complete the transformation from the membrane potential, with Destexhe et al. (1994)

$\begin{array}{ll}\left[ACH\right](t)=\frac{{\left[ACH\right]}^{max}}{1+e^{-(V_c(t)-{\theta }_c)/{\sigma }_c}},& (6)\end{array}$

where [ACH]^max is the maximum concentration of ACH.

Moreover, CH receives an external input current I_ext. Therefore, the membrane potential V_c(t) in CH obeys

$\begin{array}{ll}{\tau }_c{\dot{V}}_c(t)=-I_L^c(t)-I_{ext}(t),& (7)\end{array}$

here, I_ext has the form (Rudolph et al., 2004)

$\begin{array}{ll}{I}_{ext}=g_{ext}(t)(V_c(t)-E_c),& (8)\end{array}$

$\begin{array}{ll}\dot{g_{ext}}(t)=\frac{-(g_{ext}(t)-g_0(t))}{{\tau }_{ext}}+\sqrt{\frac{2{\sigma }_{ext}^2}{{\tau }_{ext}}}{\varphi }_C(t),& (9)\end{array}$

where g₀ is the average conductance, τ_ext is the time constant, σ_ext is the noise standard deviation (SD) value, and ϕ_C(t) denotes the independent Gaussian white noise process of unit SD and zero mean.

2.2.2. Populations PY and IN in the Cortical Module

According to Clearwater et al. (2008), we know that ACH acts on PY both pre-synaptically and post-synaptically over time scales from milliseconds to minutes.

In the post-synaptic case, ACH acts on PY through muscarinic receptors, which can cause a certain suppression of voltage-gated potassium channels (McCormick, 1992). Hence, the M-type potassium current $I_M^p$ generated by the inactivation of these channels will change the intrinsic excitability of PY. Motivated by this observation, the membrane potential V_p(t) in PY is then formulated by

$\begin{array}{ll}{\tau }_p{\dot{V}}_p(t)=-I_L^p(t)-I_{AMPA}^p(t)-I_{GABA}^p(t)-I_M^p(t),& (10)\end{array}$

where $I_L^p,I_{AMPA}^p,I_{GABA}^p$ are calculated by Equations (2)–(4). The newly added current $I_M^p$ is described as follows (Rich et al., 2018)

$\begin{array}{ll}{I}_M^p(t)=g_M(t)·(V_p(t)-E_M)z(t),& (11)\end{array}$

where z(t) represents the unitless gating variable of the ionic current conductance satisfying

$\begin{array}{ll}\dot{z}(t)=\frac{\overline{z}(t)-z(t)}{{\tau }_z}.& (12)\end{array}$

Here, $\overline{z}(t)=\frac{1}{1+e^{(-V_p(t)-39)/5}}$ stands for the voltage-sensitive steady-state activation function and τ_z is the time constant. Note that g_M in Equation (11) is not a constant anymore, but varies in term of [ACH], that is,

$\begin{array}{ll}{g}_M(t)=g_M^{*}+K_{PY}·\left[ACH\right](t),& (13)\end{array}$

where $g_M^{*}$ is a constant with the nominal value listed in Table 1.

TABLE 1

Table 1. Description and nominal values of parameters in model TC-ACH.

In the pre-synaptic case, ACH acts on PY to modulate the properties of thalamocortical synapses via nicotinic receptor, which can result in changing the thalamocortical synaptic connection strength from TC to PY (denoted by N_tp) (Gil et al., 1997). Furthermore, as summarized by Kimura (2000), the thalamocortical connection strength will increase when the value of [ACH] becomes larger. Hence, we apply the way in Clearwater et al. (2008) to reformulate the connection strength N in Equation (5) as

$\begin{array}{ll}{N}_{tp}(t)=(N_{tp}^{*}-N_{tp}^{\infty })e^{\frac{-K_{TP}·\left[ACH\right](t)}{\mu }}+N_{tp}^{\infty }.& (14)\end{array}$

Obviously, there have

$N_{tp}(t)=\{\begin{array}{cc}{N}_{tp}^{\infty },& \text{as }\left[ACH\right](t)\to \infty \text{;}\\ N_{tp}^{*},& \text{as }\left[ACH\right](t)\to 0\text{.}\end{array}$

where $N_{tp}^{\infty }$ and N^* represent the connectivity values under two extreme situations, and μ is a concentration constant.

Because in the cortical module IN contributes less than PY to the effects of cholinergic modulation on thalamocortical rhythms (Picciotto et al., 2012), in this paper, the dynamic evolution in IN is assumed to be not affected by [ACH]. That is to say, the mathematical formulation of IN obeys Equations (1)–(5).

2.2.3. Populations RE and TC in Thalamic Module

We see in Figure 1 that population RE appears to be a central hub in connecting the cortex and population CH. The electrical activity occurring in RE is important for regulating information transmission in the thalamocortical system and the shape of thalamocortical rhythms. Therefore, we first model the discharge activity of RE in this section.

On one hand, recent research indicates that thalamic reticular neurons can be directly excited by activating α7-containing nicotinic ACH receptor (Ni et al., 2016). Consequently, we introduce a new current I_ACH−n as an input of RE from CH (corresponding the blue solid line in Figure 1).

On the other hand, it has been shown that the potassium leak current and T-type calcium current are essential for the generation of spindle oscillation (Langdon et al., 2012). Based on this fact, two more currents I_LK and I_T are considered in the formulation of the membrane potential V_r(t) in RE simultaneously.

Hence, the final equation is formulated by

$\begin{array}{ll}{\tau }_r{\dot{V}}_r(t)& =-I_L^r(t)-I_{AMPA}^r(t)-I_{GABA}^r(t)-I_{ACH-n}(t)\\ & -C_m^{-1}{\tau }_r·(I_{LK}^r(t)-I_T^r(t)).& (15)\end{array}$

In Equation (15), $I_{ACH-n}^r$, $I_{LK}^r$ and $I_T^r$ are defined, respectively, as

$\begin{array}{ll}{I}_{ACH-n}(t)=g_{ACH-n}·(V_r(t)-E_{ACH-n})\omega (t)& (16)\end{array}$

where ω(t) is solved by

$\begin{array}{ll}\dot{\omega }(t)=\frac{{\omega }_{max}}{1+{(\frac{P_d}{K_{RE}^I·[ACH](t)})}^h},& (17)\end{array}$

$\begin{array}{ll}{I}_{LK}^r=g_{LK}^r·(V_r(t)-E_{LK}^r),& (18)\end{array}$

$\begin{array}{ll}{I}_T^r(t)=g_T^r(t)·(V_r(t)-E_{Ca}^r){(m(t))}^{2}h(t).& (19)\end{array}$

Here, ω denotes the proportion of open ion-channels caused by binding of ACH. g_ACH−n and E_ACH−n stand for the conductance and reversal potential, respectively. P_d and h represent the apparent dissociation constant and Hill coefficient of ACH binding to the receptors, and ω_max is the maximum conductance. m(t) and h(t) are activation and inactivation functions of T-type current. Details about Equations (16)–(19) can be seen in Baran et al. (2010), Sethuramanujam et al. (2016), and Destexhe et al. (1993). In addition, the calcium conductance $g_T^r$ in Equation (19) is recognized to be important for generating bursting oscillations in RE, whose value increases with increasing [ACH] (Fisher and Johnston, 1990). Specifically, it is expressed by

$\begin{array}{ll}{g}_T^r(t)={\overline{g}}_T^r\frac{K_{RE}^{Ca}·[ACH](t)}{K_{RE}^{Ca}·[ACH](t)+a}& (20)\end{array}$

where ${\overline{g}}_T^r$ and a are positive constants (Omori and Horiguchi, 2004).

In addition, it has been found that the GABAergic projection within RE is decreased by increasing [ACH] (Fisher and Johnston, 1990). This mechanism is then absorbed into defining the self-feedback connectivity of RE in our work and formulated by

$\begin{array}{ll}{N}_{rr}(t)=(-b·log{K}_{RE}^{Nrr}·[ACH](t)+c)·{\overline{N}}_{rr},& (21)\end{array}$

where ${\overline{N}}_{rr}$, b and c are positive constants (Omori and Horiguchi, 2004).

$\begin{array}{ll}{\tau }_p{\dot{V}}_p=-I_L^p-I_{AMPA}^p(r_{ep})-I_{GABA}^p(r_{gp})-C_m^{-1}{\tau }_p{I}_M& (22)\end{array}$

$\begin{array}{ll}{\tau }_i{\dot{V}}_i=-I_L^i-I_{AMPA}^i(r_{ei})-I_{GABA}^i(r_{gi})& (23)\end{array}$

$\begin{array}{lll}{\tau }_t{\dot{V}}_t& =& -I_L^t-I_{AMPA}^t(r_{et})-I_{GABA}^t(r_{gt})-C_m^{-1}{\tau }_t\\ & & (I_{LK}^t-I_T^t-I_h)& (24)\end{array}$

$\begin{array}{lll}{\tau }_r{\dot{V}}_r& =& -I_L^r-I_{AMPA}^r(r_{er})-I_{GABA}^r(r_{gr})-I_{ACH-n}\\ & & -C_m^{-1}{\tau }_r(I_{LK}^r-I_T^r)& (25)\end{array}$

$\begin{array}{ll}{\tau }_c{\dot{V}}_c=-I_L^c-I_{ext}& (26)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{ep}={\gamma }_e^2(N_{pp}{Q}_p+N_{tp}{\eta }_t+{\varphi }_P-r_{ep})-2{\gamma }_e{\dot{r}}_{ep}& (27)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{ei}={\gamma }_e^2(N_{pi}{Q}_p+N_{ti}{\eta }_t+{\varphi }_I-r_{ei})-2{\gamma }_e{\dot{r}}_{ei}& (28)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{et}={\gamma }_e^2(N_{pt}{\eta }_p+{\varphi }_T-r_{et})-2{\gamma }_e{\dot{r}}_{et}& (29)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{er}={\gamma }_e^2(N_{tr}{Q}_t+N_{pr}{\eta }_p-r_{er})-2{\gamma }_e{\dot{r}}_{er}& (30)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{gp}={\gamma }_g^2(N_{ip}{Q}_i-r_{gp})-2{\gamma }_g{\dot{r}}_{gp}& (31)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{gi}={\gamma }_g^2(N_{ii}{Q}_i-r_{gi})-2{\gamma }_g{\dot{r}}_{gi}& (32)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{gt}={\gamma }_r^2(N_{rt}{Q}_r-r_{gt})-2{\gamma }_r{\dot{r}}_{gt}& (33)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{gr}={\gamma }_r^2(N_{rr}{Q}_r-r_{gr})-2{\gamma }_r{\dot{r}}_{gr}& (34)\end{array}$

$\begin{array}{ll}\ddot{{\eta }_p}=v^2(Q_p-{\eta }_p)-2v{\dot{\eta }}_p& (35)\end{array}$

$\begin{array}{ll}\ddot{{\eta }_t}=v^2(Q_t-{\eta }_t)-2v{\dot{\eta }}_t& (36)\end{array}$

Next, we turn to modeling the discharge activity of TC. Different from RE, there exists one more h-type current in TC, which is responsible for the waxing and waning structure of spindle rhythms in thalamus. Hence, the membrane potential V_t(t) in TC is formulated by

$\begin{array}{ll}{\tau }_t{V}_t^{\dot{}}(t)& =-I_L^t(t)-I_{AMPA}^t(t)-I_{GABA}^t(t)-C_m^{-1}·{\tau }_t·(I_{LK}^t(t)\\ & -I_T^t(t)-I_h(t)),& (37)\end{array}$

where $I_{LK}^t,I_T^t$ are calculated by Equations (18)–(19), and I_h is described as Destexhe et al. (1996)

$\begin{array}{ll}{I}_h(t)=g_h·(V_t(t)-E_h)(m_{h1}(t)+g_{inc}·m_{h2}(t)).& (38)\end{array}$

Here, g_h and E_h represent the conductance and reversal potential, g_inc is the conductivity scaling. The details of functions m_h1(t) and m_h2(t) can be found in Destexhe et al. (1993).

2.2.4. Full Mathematical Expression of TC-ACH

There is, to emphasize, one key point in our model. On account of the long range afferent, there exist conduction delays between thalamus and cortex module. According to Costa et al. (2016), this delay is approximated as a convolution with the alpha function h(t). In this case, Equation (5) is then reformulated by

$\begin{array}{ll}{\ddot{r}}_{\xi }(t)={\gamma }_{\xi }^2(N·(h(t)\otimes Q(t))-r_{\xi }(t))-2{\gamma }_{\xi }{\dot{r}}_{\xi }(t)& (39)\end{array}$

during four transmissions “PY → RE,” “PY → TC,” “TC → PY,” and “TC → IN.”

The full mathematical expression of TC-ACH is then formulated by Equations (22)–(36) (see page 5). Note that the model output is V_p, which can be viewed as the simulated EEG signals.

Remark 1. The description and nominal values of all parameters included in the model TC-ACH are listed in Table 1, whose values are the nominal ones reported in Costa et al. (2016), Omori and Horiguchi (2004), Clearwater et al. (2008), Rich et al. (2018), Baran et al. (2010), Sethuramanujam et al. (2016), Bhattacharya et al. (2012), and Rudolph et al. (2004).

3. Results

In this section, we apply the proposed model TC-ACH to verify the effects of cholinergic modulation on thalamocortical rhythms during NREM in two ways: (1) measurement of effects of cholinergic modulation in the thalamus; (2) measurement of effects of cholinergic modulation in the thalamocortical system.

All of the numerical simulations are performed in MATLAB R2017b (MathWorks,USA), using a stochastic Runge-Kutta method of 4th order (Rößler, 2010) with a step size of 0.1ms.

3.1. The Oscillatory Phenomena Due to Cholinergic Modulation in Thalamus

The thalamus is believed to act as a “pacemaker” for thalamocortical rhythms, and is able to independently generate multiple brain rhythms during sleep (Hughes and Crunelli, 2005; Li et al., 2017). Hence, oscillatory phenomena due to cholinergic modulation in thalamus are first explored, where only the thalamic module and projections from CH to RE are considered. Here, the model output is V_t.

We first show different thalamic oscillations and transitions between them caused by cholinergic projections. To this end, we apply one-dimensional bifurcation analysis for several key parameters, represented by g_LK, $K_{RE}^I$, $K_{RE}^{Ca}$, and $K_{RE}^{Nrr}$. The bifurcation diagram is obtained by plotting the stable local minimum, as well as top three maximum values of V_t over changes in each of four parameters. All simulations are executed for 40 s and those minimum and maximum values are obtained from the latter stable 30 s of the time series.

Figures 3A–D shows bifurcation diagrams of V_t over changes in g_LK, $K_{RE}^I$, $K_{RE}^{Ca}$, and $K_{RE}^{Nrr}$, respectively. Seen from Figure 3A, it reveals that there are four different dynamical states including low firing (I), spindle (II), fast oscillation (III) and slow oscillation (IV) as g_LK varies in [0.018, 0.055]. Specifically, when the value of g_LK is extremely small, the model exhibits a low firing state and no oscillation behavior can be observed. As g_LK becomes a little larger, the model mainly experiences a spindle oscillation pattern for a period, in which multiple pairs of maximum and minimum values are found within each periodic complex. With further growth of g_LK, the model exhibits a fast oscillation pattern, in which only one pair of maximum and minimum values emerges within each periodic complex. When strong hyperpolarization through increasing g_LK attains, the model moves from a high frequency oscillation pattern into a slow oscillation pattern, in which multiple pairs of maximum and minimum values can also be observed within each periodic complex, but the distances between the top maximum and other maximum values become larger. Figure 3E illustrates the obtained thalamic oscillations, which correspond to four dynamical states (g_LK = 0.018, 0.03, 0.045, 0.052), respectively.

FIGURE 3

Figure 3. The bifurcation diagrams of V_t over changes in g_LK, $K_{RE}^I$, $K_{RE}^{Ca}$, and $K_{RE}^{Nrr}$. (A) Bifurcation diagram (varying g_LK). (B) Bifurcation diagram (varying $K_{RE}^I$). (C) Bifurcation diagram (varying $K_{RE}^{Ca}$). (D) Bifurcation diagram (varying $K_{RE}^{Nrr}$). (E) Four thalamic oscillations with respect to different values of g_LK.

A similar explanation can be given for the cases of $K_{RE}^I$, $K_{RE}^{Ca}$, and $K_{RE}^{Nrr}$, while only two dynamical states (i.e., II and III) are obtained. Moreover, the state transition of $K_{RE}^I$ is different from that of $K_{RE}^{Ca}$ and $K_{RE}^{Nrr}$. It is apparent that the model has a pattern transition from state III to an increasingly stable state II as $K_{RE}^I$ increases (see Figure 3B). By contrast, a transition from state II to state III with the increase of $K_{RE}^{Ca}$ (or $K_{RE}^{Nrr}$) can be observed in Figure 3C (or Figure 3D). That is to say, the model appears to demonstrate a stabler spindle oscillation state with the increasing value of $K_{RE}^I$, while the growth of $K_{RE}^{Ca}$ (or $K_{RE}^{Nrr}$) leads to reduced spindle rhythms. Prior work has established that spindle rhythms are helpful to protect sleep (Dang-Vu et al., 2010; Kim et al., 2012). Therefore, the obtained results implicate that cholinergic modulation in RE may help promote sleep (or arousal) states.

Next, we check whether our results can be generalized within a certain range of parameters. The above bifurcation analysis allows us to further distinguish different dynamical state regions in the two-parameter space (for example, see Figure 5A). Moreover, the power spectral analysis is applied to estimate the dominant frequency (Df) and second dominant frequency (Sec-Df) of neural oscillations from the time series of V_t. Figure 4 illustrates the box-plots of Df and Sec-Df with respect to three oscillations (fast, slow and spindle). It is clear from Figure 4 that the Sec-Df has better differentiation capability than Df. Therefore, the corresponding Sec-Df regions are drawn for each pair of parameters (for example, see Figure 5B) in the following analysis.

FIGURE 4

Figure 4. The box-plots of (A) Df and (B) Sec-Df with respect to three oscillations (fast, slow, and spindle).

FIGURE 5

Figure 5. The dynamical state regions (A,C,E) and Sec-Df regions (B,D,F) in three panels $(g_{LK},K_{RE}^I)$, $(g_{LK},K_{RE}^{Ca})$, $(g_{LK},K_{RE}^{N_{rr}})$.

Figures 5A,B illustrate the dynamical state regions and Sec-Df regions with 41 × 41 grids in the space of $g_{LK}\times K_{RE}^I\in \left[0.018,0.055\right]\times \left[1,9\right]$. As shown in Figure 5A, four different state regions are displayed, whose identification is same as above [i.e., low firing (I), spindle (II), fast oscillation (III) and slow oscillation (IV)]. It can be observed that along with $K_{RE}^I$ rising, the model generates more spindle oscillation pattern after lasting a certain fast oscillation period (see the black arrow). By combining the results of frequency analysis shown in Figure 5B, we can outline the spindle oscillation region that falls into the 11–14 Hz frequency range. It should be noted that, the obtained result is consistent with the frequency range [11, 16]Hz characteristic of sleep spindles.

Similar results are obtained in the spaces of $g_{LK}\times K_{RE}^{Ca}$ and $g_{LK}\times K_{RE}^{N_{rr}}$, which are illustrated in Figures 5C–F. We observe that increasing $K_{RE}^I$ promotes the generation of spindle oscillations in our model (see the black solid arrow in Figure 5A), while increasing $K_{RE}^{Ca}$ (or $K_{RE}^{N_{rr}}$) suppresses generation of spindle oscillations (see the black dash arrow in Figures 5C,E). These results are in line with the results mentioned above.

Additionally, in order to demonstrate the interactions among parameters $K_{RE}^I,K_{RE}^{Ca}$, and $K_{RE}^{Nrr}$, the dynamical state regions and Sec-Df regions for each pairwise combination of them are shown in Figure 6. As expected, both state and Sec-Df analysis in three panels ($K_{RE}^I$, $K_{RE}^{Ca}$), ($K_{RE}^I$, $K_{RE}^{Nrr}$), and ($K_{RE}^{Ca}$, $K_{RE}^{Nrr}$) provide the same evidence as the above description, that is, the spindle oscillation is generated as $K_{RE}^I$ increases (or as $K_{RE}^{Ca}$, $K_{RE}^{Nrr}$ decreases).

FIGURE 6

Figure 6. The dynamical state regions (A,C,E) and Sec-Df regions (B,D,F) in three panels ($K_{RE}^I$, $K_{RE}^{Ca}$), ($K_{RE}^I$, $K_{RE}^{Nrr}$), and ($K_{RE}^{Ca}$, $K_{RE}^{Nrr})$.

3.2. The Oscillatory Phenomena Due to Cholinergic Modulation in the Thalamocortical System

Here we concentrate on the effects of cholinergic modulation on thalamocortical rhythms during NREM sleep, including K-complexes, spindles, and slow oscillations. We fix the values of the four parameters to be g_LK = 0.034, $K_{RE}^I=5.2$, $K_{RE}^{Ca}=1.5$, and $K_{RE}^{Nrr}=2.1$. Based on the results obtained in the proceeding subsection, we know that spindle oscillations can be generated in the thalamus module with such settings.

We first show the dynamical behavior of our TC-ACH model with variation of K_PY, which is the strength of the cholinergic projection from CH to PY. In this case, we assume that there is no cholinergic projection from CH to “PY-TC connection” (K_TP = 0), and we apply linear ramps to increase or decrease the parameter value and observe the effects of continuous changes in K_PY on the model. Figure 7B illustrates the model output V_p as K_PY varies in [2.3, 5] linearly (see Figure 7A). In the beginning, when K_PY = 2.3, the model exhibits K-complex oscillations. Then we can find that a small ascent of K_PY drives a transition from K-complexes to slow oscillations. Furthermore, with the gradual increase of K_PY, the slow oscillation is kept for a short period until K_PY reaches its higher threshold. At this time, the slow oscillation is replaced by the α-like oscillation. As K_PY ramps down, the model transitions back into its original state gradually (that is, from α-like to slow oscillations, and then from slow oscillations to K-complexes). The enlarged graphs of three different oscillations in 20 s are displayed in Figures 7E–G.

FIGURE 7

Figure 7. The dynamical behavior of model TC-ACH with the variation of K_PY: (A) the variation of K_PY over time t; (B) the model output V_p; (C) the spectrogram of V_p; (D) the sample entropy of V_p; (E) the enlarged model output in [10,30]s (including K-complex); (F) the enlarged model output in [30,50]s (including slow oscillation); (G) the enlarged model output in [65,85]s (including α-like oscillation).

On the basis of the above observations, we conclude that K-complexes and slow oscillations can be triggered by weakening the strength of the cholinergic projection from CH to cortex, while the emergence of α-like activity, characteristic of REM, requires stronger cholinergic input. Here, if we roughly consider that K-complexes and slow oscillations mainly emerge during NREM sleep and α-like rhythms mostly appear during REM sleep, the obtained results in our work are consistent with the conclusion in Lena et al. (2005), that is, the concentration of ACH is lowest during NREM sleep and highest during wake and REM sleep.

The spectral analysis and the sample entropy extracted from the model output V_p also support this conclusion from another point of view. Figure 7C illustrates the spectrogram of V_p. During the first and last 40 s where K-complexes and slow oscillations are emerging, the spectral power of V_p mainly falls within the [0, 5]Hz range; while in the middle stage where α-like rhythms begin to appear, it increases to [0, 10]Hz. These results obviously conform with the frequency ranges of certain sleep rhythms. Besides, the sample entropy extracted from V_p also shows a similar tendency. In detail, it is larger in the α-like stage, but smaller in the K-complex and slow wave oscillation stage (see Figure 7D). Sample entropy measures signal complexity, hence this observation is consistent with the conclusion that the complexity of cortical rhythms are decreased as the deeper sleep state, but increased during REM and wake (Bruce et al., 2009; See and Liang, 2011).

Next, we show the dynamical behavior of constructed model TC-ACH with variation of K_TP, which is the strength of the cholinergic projection from CH to the connection between PY and TC. Here, K_PY shows the trend as in Figure 7, and K_TP is set to be 1.5, 3.5, 5.5, respectively, in simulations. Figure 8 illustrates the model output V_p under three different cases.

FIGURE 8

Figure 8. The dynamical behavior of model TC-ACH in the cases where K_TP = 1.5, K_TP = 3.5, K_TP = 5.5. (A1,B1,C1) The model output Vp; (A2,B2,C2) the enlarged model output in [20,40]s; (A3,B3,C3) the enlarged model output in [40,60]s; (A4,B4,C4) the enlarged model output in [65,85]s.

In the first case where the connectivity between TC and PY is regulated by a lower strength of cholinergic projection from CH (i.e., K_TP = 1.5), we observe in Figure 8A1 that some weaker spindle rhythms emerge when the value of K_PY is relatively small. It can be observed more clearly in Figure 8A2 where the enlarged graph exhibits spindle rhythms (with red circle) and K-complex oscillations in [20, 40] s. Meanwhile, in Figures 8A3,A4, we see that the duration of slow oscillation is much longer, and α-like activity is much less by comparing the results shown in Figure 7.

This conclusion can be further verified under other two cases. It can be seen in Figures 8B1–B4 that more spindle rhythms appear as K_PY gradually increases and the α-like activities have been largely replaced by slow waves in the case where K_TP = 3.5. Furthermore, when K_TP is set to be 5.5, the spindle rhythms run through the model output from the beginning to end. We can see from Figure 8C4 that there are still spindle rhythms when K_PY attains its maximum value. These observations reveal that the strength of cholinergic projections from CN to “PY-TC connection” play an important role in promoting the spindle rhythms and prolonging the duration of the NREM state in the thalamocortical system.

4. Summary and Discussion

In this paper, we first proposed a novel computational model (TC-ACH) by integrating a neuron population CH into classical thalamo-cortical circuitry (including populations PY, IN, TC, and RE). The connections between five populations are built in accordance with the established mechanisms. Our model considers the neurotransmitter ACH released by neurons in CH, which alters discharge activities of thalamocortical neurons. For simplicity, we represent these additions to the classic model by four cholinergic projections, where coupling coefficients $K_{RE}^{Ca},K_{RE}^{N_{rr}},K_{RE}^I$ represent different projections from CH to thalamus, and K_PY and K_TP represent the projections from CH to cortex and thalamocortical system, respectively. On the basis of established model framework of TC-ACH, the corresponding mathematical expression has been formulated in the logic of five populations systematically, where the average membrane potential V(t) is solved to simulate the rhythms generated by each population.

Next, we applied the developed model TC-ACH to study the effects of ACH modulation on thalamocortical rhythms during NREM sleep in two ways:

(1) Measurement of effects of cholinergic modulation in the thalamus. In this case, only the thalamic module and projections from CH to RE were considered. Simulation results suggest that cholinergic projection activity is a key factor in modulating oscillation patterns in the thalamic module. Specifically, the model appears to be a stabler spindle oscillation state with the increasing value of $K_{RE}^I$, while the growth of $K_{RE}^{Ca}$ (or $K_{RE}^{Nrr}$) leads to reduced spindle rhythms. Moreover, with variation of the potassium leak conductance g_LK, which is dramatically modulated by the concentration of ACH (McCormick, 1989, 1992), four different dynamical states including the low firing, spindle, fast oscillation, and slow oscillation can be obtained and transited.

(2) Measurement of effects of cholinergic modulation in the thalamocortical system. In this case, the dynamical behavior of our TC-ACH model was studied by varying K_PY and K_TP, respectively. Simulation results show that the K-complex and slow oscillations can be triggered by weakening strength of K_PY, while emergence of α-like activity requires stronger input. Furthermore, when there exists cholinergic input from CH to “PY-TC connection” (that is, K_TP ≠ 0), we found that the duration of oscillation patterns during NREM sleep including K-complexes, spindles and slow oscillations is longer. Additionally, when K_TP is relatively larger, more spindle rhythms appear and α-like activities are largely replaced by slow waves.

It should be noted that a number of electrophysiological experiments have investigated cholinergic modulation of cortex, thalamus or thalamocortical system, respectively (McCormick and Prince, 1986; McCormick, 1989, 1992; Clarke, 2004; Mesulam, 2004; Hasselmo and Giocomo, 2006; Beierlein, 2014). However, few studies correlate ACH modulation with rhythmic activities (Steriade et al., 1993; Steriade, 2004), where the model-based work to study in theory the effects of cholinergic modulation on thalamocortical rhythms during sleep is far less. More than that, what few existing studies are not comprehensive enough in studying the effects of ACH modulation on thalamocortical rhythms, especially during sleep. For example, the effects of ACH modulation are considered only on cortex or thalamus separately, but not on the whole thalamocortical system (Omori and Horiguchi, 2004; Li et al., 2017); the mechanism regarding ACH modulation is only modeled as a certain parameter, but not as a whole neuron population (Li et al., 2017); the model is constructed at the microscopic level, which cannot relate directly to thalamocortical rhythms at the mesoscopic level (Omori and Horiguchi, 2004). Therefore, in order to overcome such limitations, we have constructed a novel computational model (TC-ACH) by incorporating a cholinergic neuron population into the classical thalamo-cortical circuitry at the mesoscopic level. By thus doing, a deeper understanding of the role of cholinergic modulation on thalamocortical system will be got, and further, a critical insight into the mechanisms controlling sleep state may be found.

Besides five projections applied in the modeling (represented by $K_{RE}^{Ca},K_{RE}^{N_{rr}},K_{RE}^I$, K_PY, and K_TP), there still exist cholinergic modulators on other projections, such as corticoreticular and corticothalamic projections. Castro-Alamancos and Calcagnotto (2001) demonstrated that the corticothalamic activity can be filtered by cholinergic activation during arousal high-pass according to the experiments performed in vitro and in vivo. In addition, as indicated in Itier and Bertrand (2002), the cholinergic modulation on corticoreticular projection may induce the generation of more complex brain rhythms (such as spike and wave) during sleep. However, the quantitative description relating such cholinergic mechanisms to the generation and transition of typical rhythms during NREM sleep is very limited. Therefore, we hope these mechanisms could be further considered with more attempts, to study the cholinergic modulation of thalamocortical rhythms during sleep in different cases.

Another point to stress here is the synaptic connection ways in intra-RE. Previous electrophysiological studies indicated that the thalamic reticular neurons are functionally connected through chemical/electrical synapses (Sanchez-Vives et al., 1997; Landisman et al., 2002; Shu and McCormick, 2002; Long et al., 2004; Deleuze and Huguenard, 2006; Lam et al., 2006). A potentially paradigm-shifting question has been presented as to whether chemical synapses between thalamic reticular neurons are altogether absent in certain mammals or degenerate as a function of increasing age (Landisman et al., 2002; Cruikshank et al., 2010; Hou et al., 2016). However, there still have other works supporting the existence of GABAergic intrareticular synapses. For example, one of the most recent model-based study shows that intrareticular synapses, both chemical and electrical, manifest certain effects on the signal propagation and oscillation (Brown et al., 2019). In our modeling work, by means of considering the chemical synapses in intra-RE population, the measurement of effects of cholinergic modulation in thalamocortical system has been completed well. In contrast, the electrical synapses do not work effectively under the same situation. Consequently, we roughly hypothesize that GABAergic synapses may contribute much more than electrical synapses to the effects of cholinergic modulation in thalamocortical rhythms, at least at the mesoscopic level. But on the other hand, we certainly acknowledge that the electrical intrareticular synapses play important roles vis-à-vis thalamic signaling. Therefore, an important extension of the current work would be the modeling of electrical synapses between reticular neurons from the microscopic point of view, such as the model-based works in Pham and Haas (2018) and Brown et al. (2019).

An important limitation of our TC-ACH model is that it considers only ACH modulation. However, it is known that thalamocortical rhythms during sleep are also directly affected by other neuromodulators, such as noradrenalin (NE), serotonin (5-HT), histamine (HA), and dopamine (DA) from the hypothalamus and brainstem, whose concentrations vary over the night (Lena et al., 2005). Therefore, some apparent questions need to be answered: whether our method can be further developed to shed light on other sleep-related neuromodulators? If so, how to model the corresponding mechanisms and explore their effects on thalamocortical rhythms during sleep? Fortunately, based on the progressive mathematical description of sleep regulatory networks (Kumar et al., 2012; Booth et al., 2017), it may well be possible to carry our model further by constructing a new thalamocortical NMM which contains various sleep-related neuron populations (releasing NE, 5-HT, HA, DA); this is a topic we plan to pursue in future work. In addition, such model-based research can make a contribution to understanding sleep related pathological conditions, such as sleep-related epilepsy.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

QL, J-LS, S-HL, MW, and RZ designed and performed the research as well as wrote the paper.

Funding

This work was supported by the National Natural Science Foundation of China under Grant 61473223, the Innovative Talents Promotion Plan of Shaanxi Province under Grant 2018TD-016, the Key Industry Innovation Chain (Group) of Shaanxi Province under Grant 2019ZDLSF02-09-02. MW was supported during this research by the Glenn Foundation for Medical Research and the American Federation for Aging Research through a Breakthroughs in Gerontology Grant; through the American Academy of Sleep Medicine through an AASM Foundation Strategic Research Award; and by grants from the NIH (1R01NS102190, 1R01NS102574, 1R01NS107291, 1RF1AG064312).

Conflict of Interest

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

Baran, I., Iftime, A., and Popescu, A. (2010). Diffusion–convection effects on drug distribution at the cell membrane level in a patch-clamp setup. BioSystems 102, 134–147. doi: 10.1016/j.biosystems.2010.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Beierlein, M. (2014). Synaptic mechanisms underlying cholinergic control of thalamic reticular nucleus neurons. J. Physiol. 592, 4137–4145. doi: 10.1113/jphysiol.2014.277376

PubMed Abstract | CrossRef Full Text | Google Scholar

Bellingham, M. C., and Funk, G. D. (2000). Cholinergic modulation of respiratory brain-stem neurons and its function in sleep–wake state determination. Clin. Exp. Pharmacol. Physiol. 27, 132–137. doi: 10.1046/j.1440-1681.2000.03192.x

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 International Conference on Artificial Neural Networks, eds A. E. P. Villa (Berlin; Heidelberg: Springer), 645–652. doi: 10.1007/978-3-642-33269-2_81

CrossRef Full Text | Google Scholar

Booth, V., Xique, I., and Diniz Behn, C. G. (2017). One-dimensional map for the circadian modulation of sleep in a sleep-wake regulatory network model for human sleep. SIAM J. Appl. Dyn. Syst. 16, 1089–1112. doi: 10.1137/16M1071328

CrossRef Full Text | Google Scholar

Boutrel, B., and Koob, G. F. (2004). What keeps us awake: the neuropharmacology of stimulants and wakefulness promoting medications. Sleep 27, 1181–1194. doi: 10.1093/sleep/27.6.1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, J. W., Taheri, A., Kenyon, R. V., Berger-Wolf, T., and Llano, D. A. (2019). A computational model of intrathalamic signaling via open-loop thalamo-reticular-thalamic architectures. bioRXiv. 574178. doi: 10.1101/574178

CrossRef Full Text | Google Scholar

Brown, R. E., Basheer, R., McKenna, J. T., Strecker, R. E., and McCarley, R. W. (2012). Control of sleep and wakefulness. Physiol. Rev. 92, 1087–1187. doi: 10.1152/physrev.00032.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruce, E. N., Bruce, M. C., and Vennelaganti, S. (2009). Sample entropy tracks changes in EEG power spectrum with sleep state and aging. J. Clin. Neurophysiol. 26:257. doi: 10.1097/WNP.0b013e3181b2f1e3

CrossRef Full Text | Google Scholar

Castro-Alamancos, M. A., and Calcagnotto, M. E. (2001). High-pass filtering of corticothalamic activity by neuromodulators released in the thalamus during arousal: in vitro and in vivo. J. Neurophysiol. 85, 1489–1497. doi: 10.1152/jn.2001.85.4.1489

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, P. B. (2004). Nicotinic modulation of thalamocortical neurotransmission. Prog. Brain Res. 145, 253–260. doi: 10.1016/S0079-6123(03)45017-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Clearwater, J., Rennie, C., and Robinson, P. (2008). Mean field model of acetylcholine mediated dynamics in the thalamocortical system. J. Theor. Biol. 255, 287–298. doi: 10.1016/j.jtbi.2008.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cona, F., Lacanna, M., and Ursino, M. (2014). A thalamo-cortical neural mass model for the simulation of brain rhythms during sleep. J. Comput. Neurosci. 37, 125–148. doi: 10.1007/s10827-013-0493-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. S., Weigenand, A., Ngo, H.-V. V., Marshall, L., Born, J., Martinetz, T., et al. (2016). A thalamocortical neural mass model of the EEG during NREM sleep and its response to auditory stimulation. PLoS Comput. Biol. 12:e1005022. doi: 10.1371/journal.pcbi.1005022

CrossRef Full Text | Google Scholar

Cruikshank, S. J., Urabe, H., Nurmikko, A. V., and Connors, B. W. (2010). Pathway-specific feedforward circuits between thalamus and neocortex revealed by selective optical stimulation of axons. Neuron 65, 230–245. doi: 10.1016/j.neuron.2009.12.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Dang-Vu, T. T., McKinney, S. M., Buxton, O. M., Solet, J. M., and Ellenbogen, J. M. (2010). Spontaneous brain rhythms predict sleep stability in the face of noise. Curr. Biol. 20, R626–R627. doi: 10.1016/j.cub.2010.06.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Deleuze, C., and Huguenard, J. R. (2006). Distinct electrical and chemical connectivity maps in the thalamic reticular nucleus: potential roles in synchronization and sensation. J. Neurosci. 26, 8633–8645. doi: 10.1523/JNEUROSCI.2333-06.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Babloyantz, A., and Sejnowski, T. J. (1993). Ionic mechanisms for intrinsic slow oscillations in thalamic relay neurons. Biophys. J. 65, 1538–1552. doi: 10.1016/S0006-3495(93)81190-1

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Mainen, Z. F., and Sejnowski, T. J. (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

Fisher, R., and Johnston, D. (1990). Differential modulation of single voltage-gated calcium channels by cholinergic and adrenergic agonists in adult hippocampal neurons. J. Neurophysiol. 64, 1291–1302. doi: 10.1152/jn.1990.64.4.1291

PubMed Abstract | CrossRef Full Text | Google Scholar

Gil, Z., Connors, B. W., and Amitai, Y. (1997). Differential regulation of neocortical synapses by neuromodulators and activity. Neuron 19, 679–686. doi: 10.1016/S0896-6273(00)80380-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasselmo, M., and Giocomo, L. (2006). Cholinergic modulation of cortical function. J. Mol. Neurosci. 30, 133–135. doi: 10.1385/JMN:30:1:133

PubMed Abstract | CrossRef Full Text | Google Scholar

Hodgkin, A. L., and Huxley, A. F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500–544. doi: 10.1113/jphysiol.1952.sp004764

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, G., Smith, A. G., and Zhang, Z.-W. (2016). Lack of intrinsic gabaergic connections in the thalamic reticular nucleus of the mouse. J. Neurosci. 36, 7246–7252. doi: 10.1523/JNEUROSCI.0607-16.2016

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

Itier, V., and Bertrand, D. (2002). Mutations of the neuronal nicotinic acetylcholine receptors and their association with adnfle. Neurophysiol. Clin. 32, 99–107. doi: 10.1016/S0987-7053(02)00294-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Jansen, B. H., Zouridakis, G., and Brandt, M. E. (1993). A neurophysiologically-based mathematical model of flash visual evoked potentials. Biol. Cybernet. 68, 275–283. doi: 10.1007/BF00224863

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, A., Latchoumane, C., Lee, S., Kim, G. B., Cheong, E., Augustine, G. J., et al. (2012). Optogenetically induced sleep spindle rhythms alter sleep architectures in mice. Proc. Natl. Acad. Sci. U.S.A. 109, 20673–20678. doi: 10.1073/pnas.1217897109

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimura, F. (2000). Cholinergic modulation of cortical function: a hypothetical role in shifting the dynamics in cortical network. Neurosci. Res. 38, 19–26. doi: 10.1016/S0168-0102(00)00151-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobayashi, Y., Inoue, Y., and Isa, T. (2003). Cognitive function of the brainstem cholinergic system. Int. Congr. Ser. 1250, 397–406. doi: 10.1016/S0531-5131(03)00967-1

CrossRef Full Text | Google Scholar

Kumar, R., Bose, A., and Mallick, B. N. (2012). A mathematical model towards understanding the mechanism of neuronal regulation of wake-nrems-rems states. PLoS ONE 7:e42059. doi: 10.1371/journal.pone.0042059

PubMed Abstract | CrossRef Full Text | Google Scholar

Lam, Y.-W., Nelson, C. S., and Sherman, S. M. (2006). Mapping of the functional interconnections between thalamic reticular neurons using photostimulation. J. Neurophysiol. 96, 2593–2600. doi: 10.1152/jn.00555.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Landisman, C. E., Long, M. A., Beierlein, M., Deans, M. R., Paul, D. L., and Connors, B. W. (2002). Electrical synapses in the thalamic reticular nucleus. J. Neurosci. 22, 1002–1009. doi: 10.1523/JNEUROSCI.22-03-01002.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Langdon, A. J., Breakspear, M., and Coombes, S. (2012). Phase-locked cluster oscillations in periodically forced integrate-and-fire-or-burst neuronal populations. Phys. Rev. E 86:061903. doi: 10.1103/PhysRevE.86.061903

PubMed Abstract | CrossRef Full Text | Google Scholar

Lena, I., Parrot, S., Deschaux, O., Muffat-Joly, S., Sauvinet, V., Renaud, B., et al. (2005). Variations in extracellular levels of dopamine, noradrenaline, glutamate, and aspartate across the sleep–wake cycle in the medial prefrontal cortex and nucleus accumbens of freely moving rats. J. Neurosci. Res. 81, 891–899. doi: 10.1002/jnr.20602

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., Henriquez, C. S., and Fröhlich, F. (2017). Unified thalamic model generates multiple distinct oscillations with state-dependent entrainment by stimulation. PLoS Comput. Biol. 13:e1005797. doi: 10.1371/journal.pcbi.1005797

PubMed Abstract | CrossRef Full Text | Google Scholar

Long, M. A., Landisman, C. E., and Connors, B. W. (2004). Small clusters of electrically coupled neurons generate synchronous rhythms in the thalamic reticular nucleus. J. Neurosci. 24, 341–349. doi: 10.1523/JNEUROSCI.3358-03.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A. (1989). Cholinergic and noradrenergic modulation of thalamocortical processing. Trends Neurosci. 12, 215–221. doi: 10.1016/0166-2236(89)90125-2

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

McCormick, D. A., and Prince, D. A. (1986). Acetylcholine induces burst firing in thalamic reticular neurones by activating a potassium conductance. Nature 319:402. doi: 10.1038/319402a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mesulam, M.-M. (2004). The cholinergic innervation of the human cerebral cortex. Prog. Brain Res. 145, 67–78. doi: 10.1016/S0079-6123(03)45004-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, K.-M., Hou, X.-J., Yang, C.-H., Dong, P., Li, Y., Zhang, Y., et al. (2016). Selectively driving cholinergic fibers optically in the thalamic reticular nucleus promotes sleep. Elife 5:e10382. doi: 10.7554/eLife.10382

PubMed Abstract | CrossRef Full Text | Google Scholar

Omori, T., and Horiguchi, T. (2004). Dynamical state transition by neuromodulation due to acetylcholine in neural network model for oscillatory phenomena in thalamus. J. Phys. Soc. Japan 73, 3489–3494. doi: 10.1143/JPSJ.73.3489

CrossRef Full Text | Google Scholar

Pham, T., and Haas, J. S. (2018). Electrical synapses between inhibitory neurons shape the responses of principal neurons to transient inputs in the thalamus: a modeling study. Sci. Rep. 8:7763. doi: 10.1038/s41598-018-25956-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Picciotto, M. R., Higley, M. J., and Mineur, Y. S. (2012). Acetylcholine as a neuromodulator: cholinergic signaling shapes nervous system function and behavior. Neuron 76, 116–129. doi: 10.1016/j.neuron.2012.08.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Rasch, B., and Born, J. (2013). About sleep's role in memory. Physiol. Rev. 93, 681–766. doi: 10.1152/physrev.00032.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Rich, S., Zochowski, M., and Booth, V. (2018). Effects of neuromodulation on excitatory–inhibitory neural network dynamics depend on network connectivity structure. J. Nonlinear Sci. 8, 1–24. doi: 10.1007/s00332-017-9438-6

CrossRef Full Text | Google Scholar

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 | CrossRef Full Text | Google Scholar

Rößler, A. (2010). Runge–kutta methods for the strong approximation of solutions of stochastic differential equations. SIAM J. Num. Anal. 48, 922–952. doi: 10.1137/09076636X

CrossRef Full Text | Google Scholar

Roth, T. (2009). Slow wave sleep: does it matter? J. Clin. Sleep Med. 5(2 Suppl):S4.

PubMed Abstract | Google Scholar

Rudolph, M., Piwkowska, Z., Badoual, M., Bal, T., and Destexhe, A. (2004). A method to estimate synaptic conductances from membrane potential fluctuations. J. Neurophysiol. 91, 2884–2896. doi: 10.1152/jn.01223.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., Bal, T., and McCormick, D. A. (1997). Inhibitory interactions between perigeniculate gabaergic neurons. J. Neurosci. 17, 8894–8908. doi: 10.1523/JNEUROSCI.17-22-08894.1997

PubMed Abstract | CrossRef Full Text | Google Scholar

See, A. R., and Liang, C.-K. (2011). “A study on sleep eeg using sample entropy and power spectrum analysis,” in 2011 Defense Science Research Conference and Expo (DSR) (Singapore: IEEE), 1–4. doi: 10.1109/DSR.2011.6026802

CrossRef Full Text | Google Scholar

Sethuramanujam, S., McLaughlin, A. J., deRosenroll, G., Hoggarth, A., Schwab, D. J., and Awatramani, G. B. (2016). A central role for mixed acetylcholine/gaba transmission in direction coding in the retina. Neuron 90, 1243–1256. doi: 10.1016/j.neuron.2016.04.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Shu, Y., and McCormick, D. A. (2002). Inhibitory interactions between ferret thalamic reticular neurons. J. Neurophysiol. 87, 2571–2576. doi: 10.1152/jn.00850.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M. (2004). Acetylcholine systems and rhythmic activities during the waking–sleep cycle. Prog. Brain Res. 145, 179–196. doi: 10.1016/S0079-6123(03)45013-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Steriade, M., Contreras, D., Dossi, R. C., and Nunez, A. (1993). The slow (<1 hz) oscillation in reticular thalamic and thalamocortical neurons: scenario of sleep rhythm generation in interacting thalamic and neocortical networks. J. Neurosci. 13, 3284–3299. doi: 10.1523/JNEUROSCI.13-08-03284.1993

PubMed Abstract | CrossRef Full Text | Google Scholar

Suffczynski, P., Kalitzin, S., and Lopes Da Silva, F. (2004). Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience 126, 467–484. doi: 10.1016/j.neuroscience.2004.03.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Weigenand, A., Costa, M. S., Ngo, H.-V. V., Claussen, J. C., and Martinetz, T. (2014). Characterization of k-complexes and slow wave activity in a neural mass model. PLoS Comput. Biol. 10:e1003923. doi: 10.1371/journal.pcbi.1003923

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