Computational Evidence for a Competitive Thalamocortical Model of Spikes and Spindle Activity in Rolandic Epilepsy

Rolandic epilepsy (RE) is the most common idiopathic focal childhood epilepsy syndrome, characterized by sleep-activated epileptiform spikes and seizures and cognitive deficits in school age children. Recent evidence suggests that this disease may be caused by disruptions to the Rolandic thalamocortical circuit, resulting in both an abundance of epileptiform spikes and a paucity of sleep spindles in the Rolandic cortex during non-rapid eye movement sleep (NREM); electrographic features linked to seizures and cognitive symptoms, respectively. The neuronal mechanisms that support the competitive shared thalamocortical circuitry between pathological epileptiform spikes and physiological sleep spindles are not well-understood. In this study we introduce a computational thalamocortical model for the sleep-activated epileptiform spikes observed in RE. The cellular and neuronal circuits of this model incorporate recent experimental observations in RE, and replicate the electrophysiological features of RE. Using this model, we demonstrate that: (1) epileptiform spikes can be triggered and promoted by either a reduced NMDA current or h-type current; and (2) changes in inhibitory transmission in the thalamic reticular nucleus mediates an antagonistic dynamic between epileptiform spikes and spindles. This work provides the first computational model that both recapitulates electrophysiological features and provides a mechanistic explanation for the thalamocortical switch between the pathological and physiological electrophysiological rhythms observed during NREM sleep in this common epileptic encephalopathy.

1. Introduction

Rolandic epilepsy (RE) is the most common idiopathic focal childhood epilepsy, characterized by a transient period of spontaneous seizures and cognitive deficits in school age children. Clinically RE is now recognized as a mild epileptic encephalopathy, characterized by the emergence of seizures and cognitive deficits during school age years (Lee et al., 2017; Ross et al., 2020). Electrographically, RE is characterized by distinctive high-voltage spikes in the Rolandic region. Both the seizures and epileptiform spikes in RE are most prominent during non-rapid eye movement (NREM) sleep (Pavlou et al., 2012). Although the genetic, electrophysiologic, and neurological features of RE have been studied in detail (Kellaway, 2000; Archer et al., 2003; Lemke et al., 2013; Mirandola et al., 2013; Lal et al., 2014; Li et al., 2017; Dryżałowski et al., 2018), the neural mechanisms underlying the stereotyped electrophysiologic features of this disease are poorly understood.

Converging evidence suggests that RE results from focal dysfunction in the Rolandic thalamocortical circuit. In support of this hypothesis, macrostructural and microstructural white matter abnormalities have been identified adjacent to the Rolandic cortex in RE (Ciumas et al., 2014; Ostrowski et al., 2019; Thorn et al., 2020). Thalamic abnormalities have also been observed in related developmental epilepsies characterized by NREM sleep-activated spikes (Fernández et al., 2012, 2017). Most recently, children with RE were reported to have focal spindle deficits in the Rolandic cortex that predicted their cognitive deficits (Kramer et al., 2021). Consistent with experimental observations in rodent models (Clementeperez et al., 2017) and cell culture (Beenhakker and Huguenard, 2009) that pathologic spikes competitively “hijack” the thalamocortical circuit that normally generates NREM sleep spindles, spindle rate anticorrelates with spike rate in RE patients (Kramer et al., 2021). As sleep spindles originate from a well-characterized thalamocortical circuit (Beenhakker and Huguenard, 2009), these findings provide a concrete pathophysiological model for RE. The potential neural mechanisms responsible for the competitive electrophysiological dynamics between the sleep activated spikes and spindles observed in RE remain unknown.

Neural computational modeling provides a tool to investigate the underlying dynamics of observed electrophysiological phenomena, test mechanistic hypotheses, and provide predictions for further experiments. Neural mass models (NMMs) describe the dynamics of a neural population through an averaged representation of specific cell types. A variety of NMMs have been developed to simulate several types of oscillatory and epileptiform patterns of brain activity in epilepsy, including high frequency oscillations, spike-wave complexes, and polyspike and wave discharges (Wendling et al., 2002; Suffczynski et al., 2004; Marten et al., 2009; Goodfellow et al., 2011; Helling et al., 2015). Similarly, NMMs have been constructed to reproduce typical rhythms during sleep, such as K-complexes, slow oscillations and spindles (Steynross et al., 2005; Cona et al., 2014; Weigenand et al., 2014; Costa et al., 2016a). Costa et al. developed a thalamocortical modeling framework to produce realistic time courses of EEG signals during NREM sleep, accurately replicating well-characterized NREM sleep architecture, including K-complexes and spindles (Costa et al., 2016b). However, for the sleep-activated epilepsies, such as RE, theoretical modeling to demonstrate the neural mechanisms underlying epileptiform spikes during NREM sleep, has not been explored.

In this study, we introduce a new neural mass model for Rolandic epilepsy (RE-NMM) incorporating the thalamocortical circuit underlying spikes and spindles in recent human and rodent studies (Clementeperez et al., 2017; Kramer et al., 2021) and genetic findings in RE (Lemke et al., 2013; Xiong and Zhou, 2017). This model extends the existing Costa model in Costa et al. (2016b) in three ways: (1) to target the dynamics of the Rolandic thalamocortical system, we include excitatory cells from the ventroposteriomedial (VPM) and ventroposteriolateral (VPL) thalamic nuclei and two types of inhibitory cells (parvocellular, PV; and somatostatin, SOM) from the reticular nucleus of the thalamus (TRN); (2) consistent with recent genetic studies in RE and idiopathic focal epilepsies (Lemke et al., 2013; Xiong and Zhou, 2017), an extra NMDA receptor-mediated current (NMDA current) is introduced in PV and SOM populations; (3) to reflect the role of ion channels mutations in idiopathic epilepsies (Poolos, 2004; Difrancesco and Difrancesco, 2015; Brennan et al., 2016), an h-type ionic current is included in the VPL and VPM populations.

We then use the model to propose the neuronal mechanisms that give rise to epileptiform spikes in NREM sleep in RE. Finally, we use the model to characterize the neuronal mechanisms responsible for the competitive dynamics observed between spikes and spindles in this disease. This work provides a new computational model to better study RE and related sleep-activated epilepsies and proposes testable hypotheses of the neuronal mechanisms underlying this common disease.

2. Methods

2.1. The Costa NMM

We start with the thalamocortical model in Costa et al. (2016b) which contains four neural populations: excitatory pyramidal cells (PY) and inhibitory interneurons (IN) in the cortical module, and excitatory thalamocortical cells (TC) and inhibitory thalamic reticular cells (TRN) in the thalamic module. The model's topological structure is shown in Figure 1.

FIGURE 1

Figure 1. The topological structure of Costa model.

In each population the dynamical evolution of neural activity is formulated by two computational operators: a sigmoid operator and a convolution operator. The sigmoid operator transforms the average membrane potential V(t) into the average firing rate Q(t), that is (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. The convolution operator maps the firing rate Q(t) into the proportion of open synaptic channels r_ξ(t). The membrane potential V(t) is then given by

$\begin{array}{l}\tau \dot{V(t)}=-I_L(t)-{\displaystyle \sum _{\xi }}{I}_{\xi }(t)\\ =-{\overline{g}}_L·(V(t)-E_L)-{\displaystyle \sum _{\xi }}{\overline{g}}_{\xi }{r}_{\xi }·(V(t)-E_{\xi }).& (2)\end{array}$

In (2), the proportion of open synaptic channels is given by

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

where h_ξ(t) is the impulse response function, denoted by

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

Here, ξ ∈ {e, i} denotes the type of synapse with e standing for AMPAergic synapses type for excitation and i for GABAergic type for inhibition. In Equations (2)–(4), E_L and E_ξ represent the reversal potential of the leak current and synaptic current, respectively, g_L is the maximal conductivity of the leak current conductance, ${\overline{g}}_{\xi }$ stands for the synaptic input rate that scales r_ξ, τ is the membrane time constant, N denotes the connectivity constant, and γ_ξ is the rate constant of synaptic response. The dynamics of each population can be illustrated as in Figure 2.

FIGURE 2

Figure 2. The dynamical evolution in each population in Costa model.

Note that several types of currents play important roles in the Costa model. Specifically, I_L and I_ξ in Equation (2) represent the leak and synaptic current, respectively. In addition, several intrinsic currents are also incorporated into the calculation of V(t), including the sodium dependent potassium current I_KNa, potassium leak current I_LK, h-type current I_h and T-type calcium current I_T.

2.2. Building the Rolandic NMM

2.2.1. Rolandic Thalamocortical Neural Circuit

Somatosensory thalamocortical neurons from the ventrobasal (VB) thalamus process sensory information delivered between the periphery and the cortex, and also modulate thalamocortical states including sleep and seizures (Kramer et al., 2021). The VB consists of the ventral posteromedial nucleus (VPM) and the ventral posterolateral nucleus (VPL) (Jones, 2007). VPL neurons relay somatosensory, proprioceptive, and nociceptive information from the body to the somatosensory cortex (Francis et al., 2008). The VPM transmits similar information for the face (Iavarone et al., 2019). Neurons in VPL and VPM receive glutamatergic input from the somatosensory cortex, and inhibitory input from the thalamic reticular nucleus (TRN).

The TRN serves as a gate for information flow between the cortex and thalamus, which is important for supporting sleep regulation and seizure disruption. The TRN is entirely composed of GABAergic neurons and is the main source of inhibition for the VB (Gentet and Ulrich, 2003). The TRN conversely receives excitatory input from relay cells of VB, as well as from pyramidal cortical neurons. The TRN includes primarily parvalbumin-expressing neurons (PV) (Csillik et al., 2005), but somatostatin-positive neurons (SOM) have also been reported (Ahrens et al., 2015). Recently, PV and SOM were shown to coexist in the human TRN and have distinct electrophysiological properties, segregate into different anatomical locations, and participate in predominantly non-overlapping anatomical pathways that distinctly modulate somatosensory thalamocortical circuits both in vitro and in vivo (Clementeperez et al., 2017).

2.2.2. Rolandic Thalamocortical Cellular Circuit

Recent evidence supports a genetic basis for Rolandic epilepsy (Lemke et al., 2013; Lal et al., 2014; Xiong and Zhou, 2017), with alterations of the gene encoding for α2 subunit of NMDARs (GRIN2A) as a major genetic risk factor for RE and related severe epileptic encephalopathies. The NMDAR is a glutamate-bound excitatory receptor with important roles in the discharge activity of neurons (Paoletti, 2011). Mutations in GRIN2A are thought to affect the NMDAR properties, including Mg²⁺ block levels and channel gating kinetics, giving rise to sleep-activated epileptic encephalopathies, including RE, continuous spike and wave of sleep with encephalopathy, and Landau Kleffner syndrome (Lemke et al., 2013; Reif et al., 2017; Xu and Luo, 2018). In addition, numerous studies have also suggested that other idiopathic epilepsies may be linked to a variety of ion channel mutations, which can lead to a change in the corresponding ionic current, such as the h-type current (Poolos, 2004; Difrancesco and Difrancesco, 2015; Brennan et al., 2016).

2.2.3. A Neural Mass Model of Rolandic Epilepsy (RE-NMM)

We introduce a new computational model RE-NMM based on the Costa NMM updated with features to reflect the somatosensory thalamocortical circuit and recent genetic findings in RE. First, we modify the Costa NMM to replicate features of RE by adding reticular nucleus cells and ventrobasal thalamus cells, populations present in that Rolandic thalamocortical circuit and implicated in the generation of spikes and spindles. These consist of one population in the cortical module (pyramidal cells, PY) and four populations in the thalamic module (PV, SOM, VPM, and VPL). Based on known neuroanatomical circuitry, the cortical PY population sends excitatory projection to inhibitory thalamic PV cells in the TRN and interacts mutually with excitatory thalamic VPM and VPL. Within the thalamic module, each pair of populations among PV, VPM, and VPL have excitatory or inhibitory projections; the SOM population in the TRN sends inhibitory projections to the VPL and the PV, and receives excitatory input from the VPL (Clementeperez et al., 2017). In addition, there are self-connections within PY, PV, and SOM populations, and external inputs to PY, SOM, and VPM. Note that the excitatory projections are mediated by AMPA or NMDA receptors, and inhibitory projections are mediated by GABA receptors in our model. The topological structure of our proposed RE-NMM is illustrated in Figure 3.

FIGURE 3

Figure 3. Schematic representation of model RE-NMM.

Motivated by the discovery that the excitatory NMDA receptors contribute to RE, we next add an additional NMDA current (I_NMDA) involving two inhibitory populations PV and SOM in thalamus, formulated by (Zador et al., 1990; Li and Cleland, 2017).

$\begin{array}{ll}{I}_{NMDA}=g_N·r·B(V_k(t))·(V_k(t)-E_N).& (5)\end{array}$

Here, k ∈ {PV, SOM}, V_k is the membrane potential of PV or SOM. The function B(V) implements the Mg²⁺ block for the NMDA current, denoted by

$\begin{array}{ll}B(V)=\frac{1}{1+\left[{Mg}^{2+}\right]e^{(-\mu ·0.062·V)}/3.57}.& (6)\end{array}$

To characterize the effects of GRIN2A mutations on NMDA receptor, a new parameter μ is added in the expression of B(V).

Finally, the h-type current in ventrobasal thalamus neurons is a potential regulation factor for spike activities. Therefore, both VPM and VPL are equipped with the h-type current I_h, as described in Destexhe et al. (1996)

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

where t ∈ {VPM, VPL}, V_t is the membrane potential of VPL or VPM. The value of the conductance g_h is assumed to be affected by ion channel mutations. E_h is the reversal potential, g_inc denotes the conductivity scaling. The details of functions m_h1(t) and m_h2(t) can be found in Destexhe et al. (1993).

We model the conduction delay between the cortical and thalamic modules as a convolution with an alpha function $\tilde{h}(t)$ (Costa et al., 2016b). That is, Equation (3) is reformulated as

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

during five transmissions “PY → PV”, “PY → VPM/VPL,” and “VPM/VPL → PY.”

Equations (9)–(29) present the full mathematical expression of the RE-NMM (see page 5–6). Note that the model output is V_p, which is viewed as the simulated EEG signal (Hocepied et al., 2013; Costa et al., 2016b). The parameters used in RE-NMM are listed in Table 1, whose values and ranges are set based on previous studies (Zador et al., 1990; Traub et al., 1991; Destexhe et al., 1993, 1996; Costa et al., 2016b; Clementeperez et al., 2017; Li and Cleland, 2017; Iavarone et al., 2019).

$\begin{array}{ll}{\tau }_p\dot{V_p}=-I_L^p-I_{AMPA}^p(r_e^p),& (9)\end{array}$

$\begin{array}{l}{\tau }_{pv}\dot{V_{pv}}=-I_L^{pv}-I_{AMPA}^{pv}(r_e^{pv})-I_{NMDA}^{pv}(r_{en}^{pv})-I_{GABA}^{pv}(r_i^{pv})\\ -C_m^{-1}{\tau }_{pv}(I_{LK}^{pv}+I_T^{pv}),& (10)\end{array}$

$\begin{array}{l}{\tau }_{som}\dot{V_{som}}=-I_L^{som}-I_{AMPA}^{som}(r_e^{som})-I_{NMDA}^{som}(r_{en}^{som})-I_{GABA}^{som}(r_i^{som})\\ -C_m^{-1}{\tau }_{som}(I_{LK}^{som}+I_T^{som}),& (11)\end{array}$

$\begin{array}{l}{\tau }_{vpm}\dot{V_{vpm}}=-I_L^{vpm}-I_{AMPA}^{vpm}(r_e^{vpm})-I_{GABA}^{vpm}(r_i^{vpm})\\ -C_m^{-1}{\tau }_{vpm}(I_{LK}^{vpm}+I_T^{vpm}+I_h^{vpm}),& (12)\end{array}$

$\begin{array}{l}{\tau }_{vpl}\dot{V_{vpl}}=-I_L^{vpl}-I_{AMPA}^{vpl}(r_e^{vpl})-I_{GABA}^{vpl}(r_i^{vpl})\\ -C_m^{-1}{\tau }_{vpl}(I_{LK}^{vpl}+I_T^{vpl}+I_h^{vpl}),& (13)\end{array}$

$\begin{array}{ll}{I}_{LK}^{\alpha }=g_{LK}·(V_{\alpha }-E_{LK}),\alpha \in \left\{pv,som,vpm,vpl\right\},& (14)\end{array}$

$\begin{array}{ll}{I}_T^{\beta }=g_T^{\beta }{m}^{2}h·(V_{\beta }-E_{Ca}),\beta \in \left\{pv,som,vpm,vpl\right\},& (15)\end{array}$

$\begin{array}{l}{\ddot{r}}_e^p={\gamma }_e^2(N_{pp}{Q}_p+N_{pt1}{\eta }_{t1}+N_{pt2}{\eta }_{t2}+{\varphi }_p-r_e^p)\\ -2{\gamma }_e{ṙ}_e^p,& (16)\end{array}$

$\begin{array}{l}{\ddot{r}}_e^{vpm}={\gamma }_e^2(N_{t1p}{\eta }_p+N_{t1t2}{Q}_{vpl}+{\varphi }_t-r_e^{vpm})\\ -2{\gamma }_e{ṙ}_e^{vpm},& (17)\end{array}$

$\begin{array}{ll}{\ddot{r}}_e^{vpl}={\gamma }_e^2(N_{t2p}{\eta }_p+N_{t2t1}{Q}_{vpm}-r_e^{vpl})-2{\gamma }_e{ṙ}_e^{vpl},& (18)\end{array}$

$\begin{array}{l}{\ddot{r}}_e^{pv}={\gamma }_e^2(N_{r1p}{\eta }_p+N_{r1t1}{Q}_{vpm}+N_{r1t2}{Q}_{vpl}-r_e^{pv})\\ -2{\gamma }_e{ṙ}_e^{pv},& (19)\end{array}$

$\begin{array}{ll}{\ddot{r}}_e^{som}={\gamma }_e^2(N_{r2t2}{Q}_{vpl}+{\varphi }_n-r_e^{som})-2{\gamma }_e{ṙ}_e^{som},& (20)\end{array}$

$\begin{array}{l}{\ddot{r}}_{en}^{pv}={\gamma }_g{\gamma }_n(N_{r1p}{\eta }_p+N_{r1t1}{Q}_{vpm}+N_{r1t2}{Q}_{vpl}-r_{en}^{pv})\\ -2{\gamma }_n{ṙ}_{en}^{pv},& (21)\end{array}$

$\begin{array}{ll}{\ddot{r}}_{en}^{som}={\gamma }_g{\gamma }_n(N_{r2t2}{Q}_{vpl}-r_{en}^{som})-2{\gamma }_n{ṙ}_{en}^{som},& (22)\end{array}$

$\begin{array}{ll}{\ddot{r}}_i^{vpm}={\gamma }_N^2(N_{t1r1}{Q}_{pv}-r_i^{vpm})-2{\gamma }_N{ṙ}_i^{vpm},& (23)\end{array}$

$\begin{array}{ll}{\ddot{r}}_i^{vpl}={\gamma }_N^2(N_{t2r1}{Q}_{pv}+N_{t2r2}{Q}_{som}-r_i^{vpl})-2{\gamma }_N{ṙ}_i^{vpl},& (24)\end{array}$

$\begin{array}{ll}{\ddot{r}}_i^{pv}={\gamma }_N^2(N_{rr1}{Q}_{pv}+N_{r2r1}{Q}_{som}-r_i^{pv})-2{\gamma }_N{ṙ}_i^{pv},& (25)\end{array}$

$\begin{array}{ll}{\ddot{r}}_i^{som}={\gamma }_N^2(N_{rr2}{Q}_{som}-r_i^{pv})-2{\gamma }_N{ṙ}_i^{som},& (26)\end{array}$

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

$\begin{array}{ll}{\ddot{\eta }}_{t1}=v^2(Q_{vpm}-{\eta }_{t1})-2v{\dot{\eta }}_{t1},& (28)\end{array}$

$\begin{array}{ll}{\ddot{\eta }}_{t2}=v^2(Q_{vpl}-{\eta }_{t2})-2v{\dot{\eta }}_{t2}.& (29)\end{array}$

TABLE 1

Table 1. The parameters and their values in RE-NMM.

2.3. Empirical Fitting of Model Parameters

To reflect the abundant spike activity common in RE, we utilized a dataset of EEGs collected from 5 children with active RE during NREM sleep (ages 4.9–14.7, 5M) with a minimum of 3 epileptiform spikes per minute during NREM sleep to constrain the parameters of the model. Clinical EEG data were referenced to the average reference. We then applied a second-order Butterworth filter with a high pass of 1 Hz, a low pass of 50 Hz, and a notch at 60 Hz to denoise RE-EEGs. To capture the focal Rolandic features, we utilized data from six centrotemporal channels: C3, C4, C5, C6, T3, and T4. Each EEG signal in the selected channels was segmented into 1-min RE-EEG epochs.

For spike detection in each RE-EEG, we applied an automated and validated spike detector, Persyst 13 (Persyst Development Corporation, San Diego) (Scheuer et al., 2017). Only epochs in which spike rates were equal to or larger than 3/min were included (n = 85), denoted by S_spike. Table 2 summarizes descriptive information about the clinical dataset.

TABLE 2

Table 2. Detail information of RE-EEGs.

We then extracted features from the EEG recordings to fit the parameters of our model. Here we calculated two features: the power spectra and the histogram of inter-spike intervals. To compute the power spectra, given an EEG signal S with sampling rate F_s, the power spectral of S is calculated using the multitaper method (Babadi and Brown, 2014), denoted by

$\begin{array}{ll}\text{P}=(P_1,P_2,\cdots ,P_N).& (30)\end{array}$

Here, P_i is the power spectral at the i frequency point with $N=\left[\frac{F_s}{2}\right]$. To compute the histogram of inter spike intervals we sorted the intervals into equal-sized bins M, denoted by

$\begin{array}{ll}\text{D}=(D_1,D_2,\cdots ,D_M).& (31)\end{array}$

Here, D_j denotes the probability of inter-spike intervals in the jth bin. Figure 4 illustrates these data analysis steps on an example clinical EEG epoch from one subject.

FIGURE 4

Figure 4. Example data analysis pipeline. (A) NREM sleep data segment from one clinical EEG. The red circles represent detected spikes. (B) The power spectra computed from a 60 s epoch. We note that in this epoch where spikes are abundant, the expected sigma bump reflecting sleep spindles is absent. (C) The histogram of inter-spike intervals computed from the same data epoch.

Next, we estimate the posterior distribution for the parameters of our RE-NMM using the method in Hartoyo et al. (2019). Let θ be a parameter of RE-NMM and π(θ) be the prior distribution of θ, then the posterior distribution

$\begin{array}{ll}p(\theta |\text{D},\text{P})=\frac{p(\text{D}|\theta )p(\text{P}|\theta )\pi (\theta )}{p(\text{D},\text{P})}& (32)\end{array}$

of θ is evaluated using the Markov chain Monte Carlo (MCMC) approach. In (32), we apply the product of distributions of P_i to compute p(P|θ) (Thomson and Haley, 2014; Hartoyo et al., 2019), that is,

$\begin{array}{ll}p(\text{P}|\theta )={(\frac{K^K{e}^{-K}}{\Gamma (K)})}^N{(\alpha )}^{NK}({\displaystyle \prod _{i=1}^N}\frac{P_i^{K-1}}{\stackrel{-}{P_i}{(\theta )}^K}).& (33)\end{array}$

The likelihood function p(D|θ) is defined as

$\begin{array}{ll}p(\text{D}|\theta )={\displaystyle \prod _{i=1}^H}\frac{1}{\sqrt{2\pi }{\sigma }^2}{e}^{\frac{{({\stackrel{-}{D}}_i(\theta )-D_i)}^2}{2{\sigma }^2}}& (34)\end{array}$

with the assumption that $\stackrel{-}{\text{D}}=D+ϵ$. In this work, the error ϵ is assumed to obey the Gaussian distribution [i.e., ϵ ~ N(0, σ²)]. Here, $\stackrel{-}{\text{P}}$ and $\stackrel{-}{\text{D}}$ stands for the extracted features from the simulated EEG signals (say, the model output). A detailed explanation of each symbol in Equations (33) and (34) can be found in Moussaoui et al. (2006), Hartoyo et al. (2019).

The parameter value can be estimated by selecting the point from the Markov chain with the largest value of the posterior distribution, that is, $\widehat{\theta }=\underset{\theta }{\text{arg max }}p(\theta |\text{D},\text{P})$.

All numerical simulations are performed using MATLAB R2017b using the stochastic Runge-Kutta method of 4th order (Rosler, 2010) with step size of 0.1 ms.

3. Results

3.1. Rolandic NREM Spike Rate Is Influenced by NMDA Receptor Activity and Ionic Channel Currents

We test two parameters to generate Rolandic spikes in our model: μ and g_h, reflecting the density of NMDA current activity and the h-type current primarily influenced by ionic channel activity. All other parameters in the model are set to their nominal values (see Table 1).

We find that a range of dynamical spike activities in the RE-NMM emerge by varying the values of μ and g_h respectively (for example, see Figure 5). Increasing μ from 0.5 to 2 results in the emergence of spike events (see Figures 5A,B). Decreasing g_h from 0.066 to 0.056 also increases spike rate (Figures 5C,D). Decreasing g_h results in an increased spike rate (Figure 5D). The histograms of inter-spike intervals from the simulated spikes (Figures 5B–D) approximates that observed in the clinical EEGs (Figure 4C).

FIGURE 5

Figure 5. The simulated EEGs (i.e., the model output V_p) as well as the histogram of inter-spike intervals under different situations. (A) μ = 0.5 and g_h = 0.066; (B) μ = 2 and g_h = 0.066; (C) μ = 0.5 and g_h = 0.056; (D) μ = 2 and g_h = 0.056.

To characterize how spike rate varies with adjustments to parameters μ and g_h, we simulated EEG signals across a range of plausible values and computed the corresponding spike rates. This allows us to observe the evolution of spikes in the parameter space of μ and g_h. Figure 6 shows the distribution of spike rates with 41 × 41 grids in the space of μ × g_h ∈ [0.1, 3] × [0.046, 0.066]. We observe that higher spike rates occur for larger values of μ and smaller values of g_h, and lower spike rates occur with smaller values of μ and larger values of g_h. Note that larger values of μ and smaller values of g_h correspond to lower NMDA and h-type currents respectively (see Figures 13A,C). Our simulations indicate that reducing either the NMDA current or the h-type current causes an increase in epileptiform spikes.

FIGURE 6

Figure 6. The evolution of spike rates in the parameter space μ × g_h.

3.2. Competition Between Rolandic Spikes and Spindles Is Induced by TRN Inhibition

We note that spindles can be generated by the Costa model and our RE-NMM by adjusting parameters g_T, g_h, and g_LK, which are the conductance of T-type Calcium current, h-type current and potassium leak current. The dominant frequency of our simulated spindles is around 11 Hz, consistent with the sigma frequency range of human sleep spindles. An example simulated time series with a detected spindle (Wamsley et al., 2012) and accompanying power spectra is shown in Figure 7.

FIGURE 7

Figure 7. The simulated EEG signal and its power spectra. (A) 30s EEG segment; (B) power spectra (the red labeled part corresponds to the spindles).

Based on our observations in 3.1 and previous results (Destexhe et al., 1996; Żygierewicz et al., 2001), we first fix the values of the four parameters to be μ = 2, g_h = 0.058, g_T = 2.9, and g_LK = 0.042, to generate spikes and spindle oscillations simultaneously in simulation signals. To test for a competitive interaction between spikes and spindles during NREM sleep in our model, we then focus on the inhibitory outputs from reticular nucleus of the thalamus, which are responsible for modulating spindles (Żygierewicz et al., 2001; Beenhakker and Huguenard, 2009). Here, we compare spike and spindle rates for three inhibitory projection strengths N_t1r1, N_t2r1, and N_t2r2 (from population PV to VPM, PV to VPL, and SOM to VPL).

We implement the simulation in parameter space N_t1r1 × N_t2r1 × N_t2r2 ∈ [1, 5.5] × [1, 2] × [1, 2] with 41 × 41 × 41 grids. For ease of visualization, spike rates and spindle rates are shown for every fifth of the 41 points in Figure 8. We find that the spike rates decrease (increase) as the inhibitory projection becomes stronger (weaker), while there is the opposite relationship for spindle rate under the same situations (Figures 8A,B). The monotonously decreasing trend of the linear fit (red line) reveals that spike rate decreases with rising spindle rate, and vice versa, namely they are anticorrelated (Figure 9). Figures 9B,C show details of two points A and B in Figure 9A (marked by the square). The obtained result shows that the competitive relationship between spikes and spindles can be induced by changes of inhibitory transmission occurring in thalamus. This is consistent with the competitive relationship observed between spikes and spindles in patients with Rolandic epilepsy (see Figure 5 in Kramer et al., 2021).

FIGURE 8

Figure 8. The evolution of spike rate (A) and spindle rate (B) in parameter space N_t1r1 × N_t2r1 × N_t2r2 ∈ [1, 5.5] × [1, 2] × [1, 2].

FIGURE 9

Figure 9. (A) The combinations of spike rate and spindle rate in the diagonal grids of the parameter space N_t1r1 × N_t2r1 × N_t2r2 ∈ [1, 5.5] × [1, 2] × [1, 2]. (B) 60s simulated EEG segment with high spindle rate and low spike rate (corresponding to the point A); (C) 60s simulated EEG segment with high spike rate and low spindle rate (corresponding to the point B) (the red labeled part represents the spindles, the green circle represents the spikes).

To evaluate the impact of the NMDA currents on spindle rate, we tested the μ values ranging from 0.1 to 3. We found that spindle rate remained near constant at all values (4.5 ± 2.3).

3.3. PV and SOM Inhibitory Thalamic Populations Differentially Modulate Spikes and Spindle Activity

We included both PV and SOM neurons in our thalamic inhibitory populations. These inhibitory thalamic neurons have different pre- and post-synaptic connectivity, in that the somatosensory cortex exclusively targets the PV neurons, whereas subcortical structures preferentially target the SOM cells, and PV and SOM neurons project to distinct thalamic relay nuclei. In our model, we observed that the PV neurons have much higher firing rates than the SOM neurons (Figure 10A), consistent with Clementeperez et al. (2017).

FIGURE 10

Figure 10. (A) PV cells (top) have increased firing compared to SOM cells (bottom). Boxplots of Kruskal-Wallis test statistic for the spindle rate (B) and spike rate (C) with respect to three cases (PV+SOM, only PV, only SOM).

To evaluate the contributions of the PV and SOM neurons to the model, we compared the output from the original model including both PV and SOM populations, to models with only the PV population and only the SOM inhibitory thalamic populations. We found that models including only SOM neurons had lower spindle rates (Figure 10B) and higher spike rates (Figure 10C), suggesting that the PV neurons, but not the SOM neurons, preferentially support sleep-spindles and reduce epileptiform spikes.

3.4. Validation of Mechanisms Underlying RE Using Real RE-EEGs

To observe how Rolandic spikes in real RE-EEGs evolve in the parameter space μ × g_h, we select 27 segments from dataset S_spike in which all possibilities of spike rates (from 3 to 29) are represented. Values of μ and g_h were estimated using an MCMC-based approach using the selected RE-EEG segments (Figure 11). We observe that RE-EEG segments with higher spike rate localize to the upper-left corner of the graph (corresponding to large $\widehat{\mu }$ and small ĝ_h), and RE-EEG segments at the bottom-right corner of the graph (corresponding to small $\widehat{\mu }$ and large ĝ_h) have relatively lower spike rates.

FIGURE 11

Figure 11. The spike rate distribution of 27 real RE-EEG segments in the estimated parameter space $\widehat{\mu }\times {ĝ}_h$.

Furthermore, we display more details implicated in Figure 11 taking one point as an example (marked by the circle). The real RE-EEG segment at this point is shown in Figure 12A, from“patient 1”. Figure 12B illustrates the simulated RE-EEG segment, which is obtained as the model output V_p setting values for μ, g_h to their optimized values that $\widehat{\mu }=2.3411$ and $\widehat{g_h}=0.0522$ and using nominal values for other parameters. The power spectra of simulated EEG is in good agreement with that computed from the real data (Figure 12C) and the distribution of inter-spike intervals (ISI) of both simulated EEG and real one are mainly concentrated at lower inter-spike intervals, reflecting the higher spike rate in simulated EEG and real RE-EEG (Figure 12D).

FIGURE 12

Figure 12. An example of the model fit to one real RE-EEG segment. (A) A 60s real RE-EEG segment from “patient 1”; (B) The simulated 60s EEG segment; (C,D) Comparison between the features extracted from the real EEG and the simulated one.

These results are plausible from a physiological point of view. Specifically, the larger the value of μ, the less the density of excitatory current I_NMDA (see Figure 13A), thus the weaker the inhibitory activity in the TRN (shown as the green dot in Figure 13B), the stronger the excitatory activity in the VB (shown as the red dot in Figure 13B); the reduced I_h (lower value of g_h, Figure 13C) causes VB to be more excitable, which sends more excitation to TRN, and causes the TRN spike rates to increase (Figure 13D),—leading to a higher spike rate (see Figure 6). Notably, the firing rates of population VB decrease first (Figure 13D) because the reduced I_h could result in the increase of the time course and amplitude of after hyperpolarization potential. These observations provide new insights into how alterations in the NMDA receptor and h-type currents can result in Rolandic spikes.

FIGURE 13

Figure 13. The density of I_NMDA (A) and firing rate of population TRN and VB (B) with the increase of μ; The density of I_h (C) and firing rate of TRN and VB (D) with the decrease of g_h.

4. Summary and Discussion

In this study we introduce a computational thalamocortical model of Rolandic epilepsy informed at both the neural circuit and cellular levels by experimental observations. We validate that the model generates spike dynamics comparable to that observed in human EEG. We then use the model to demonstrate that Rolandic spikes can be triggered and promoted by a reduced NMDA current to the inhibitory thalamic cells or h-type current in the excitatory thalamic cells; and that changes in inhibitory transmission in thalamus lead to a dynamic switch between epileptiform spikes and spindles in the shared thalamocortical circuit.

Our modeling results highlight the pacemaker role of TRN, the primary inhibitory nucleus in the thalamus, in producing spike discharges and spindles during sleep or seizure processes. Previous models- to study the spike-like activities and spindle rhythms simultaneously have primarily focused on absence epilepsy, not sleep activated syndromes. For example, Suffczynski et al. (2004) developed a thalamocortical NMM to explain the relation between mechanisms that generate spindle-like activity and those that generate spike-wave activity. Zhao and Robinson (2015) extended a thalamocortical model with bursting dynamics to explore the mechanisms underlying spike and wave seizures, as well as sleep spindles etc. Fan et al. (2017) considered a single-compartment thalamocortical model and spatially extended 3-compartment coupled network to explore the role of TRN in regulating spindles and spike-wave discharges. Knox et al. (2018) introduced a thalamocortical model to understand the mechanisms of the transformation of sleep spindles to spike and wave discharges.

Our study provides the first computational model that recapitulates thalamocortical circuit “competition” between spikes and spindles in NREM sleep. That spikes may “hijack” thalamocortical spindle circuits has been proposed in theoretical papers based on experimental observations (Beenhakker and Huguenard, 2009). We provide a mechanistic model to both test and explain these observations in RE. For example, our model suggests that decreased inhibition in the TRN may explain the onset of sleep-activated spikes and seizures in these patients. These findings are consistent with common empirical treatment approaches to use GABAergic medications to treat sleep activated spikes (Sánchez Fernández et al., 2013).

We found that reducing NMDA current and TRN output increases spike rate. In contrast, reducing TRN output decreases spindle rate, but changes to NMDA current do not impact spindle rate. Thus, changes of inhibitory outputs from the reticular nucleus can result in a competitive relationship between epileptiform spikes and sleep spindles, but changes to NMDA alone do not appear to impact this dynamic. As we tested only a subset of potential mechanisms several other potential factors, including changes to the excitatory output from thalamocortical neurons (Żygierewicz et al., 2001) or slow currents (Zhao and Robinson, 2015) could contribute to the generation of epileptiform spikes and spindles, the competitive dynamics, or their properties, such as the distribution of inter-spike intervals or the waxing-and-waning structure of the spindle oscillation.

Previous electrophysiological studies indicate that enhancing SK-channel activity promotes rhythmic bursting in TRN neurons (Wimmer et al., 2012) and the cyclical Ca-SK channel interaction may be necessary for spindle generation (Cueni et al., 2008). In contrast, we were able to produce spindles and replicate an antagonistic relationship between spikes and spindles by without including the SK-channel, but by tuning just three parameters reflecting the potassium leak current, h-type current and T-type calcium current (Destexhe et al., 1996; Costa et al., 2016b) suggesting that the SK type current is not required for these dynamics. However, SK channels may play an important role in more subtle dynamics of thalamocortical rhythms, such as bursting activity (Ritter-Makinson et al., 2019), that were not explored here.

5. Conclusion

Our study provides the first computational model that both recapitulates and provides a mechanistic explanation for the thalamocortical “competition” between epileptiform spikes and sleep spindles in the most common epileptic encephalopathy. These data provide hypotheses for empirical testing of the neural mechanisms underlying this disease and related sleep-activated epilepsie.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

QL, MW, RZ, and CC designed and performed the research as well as wrote the paper. All authors contributed to the article and approved the submitted version.

Funding

QL and RZ were supported by the National Natural Science Foundation of China under Grant 12071369, the Innovative Talents Promotion Plan of Shaanxi Province under Grant 2018TD-016, and the Key Industry Innovation Chain (Group) of Shaanxi Province under Grant 2019ZDLSF02-09-02. CC was supported by NIH NINDS R01NS115868. MW was supported by NIH (1R01NS102190, 1R01NS102574, 1R01NS107291, and 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

Ahrens, S., Jaramillo, S., Yu, K., Ghosh, S., Hwang, G., Paik, R., et al. (2015). ErbB4 regulation of a thalamic reticular nucleus circuit for sensory selection. Nat. Neurosci. 18, 104–111. doi: 10.1038/nn.3897

PubMed Abstract | CrossRef Full Text | Google Scholar

Archer, J. S., Briellman, R. S., Abbott, D. F., Syngeniotis, A., Wellard, R. M., and Jackson, G. D. (2003). Benign epilepsy with centro-temporal spikes: spike triggered fMRI shows somato-sensory cortex activity. Epilepsia 44, 200–204. doi: 10.1046/j.1528-1157.2003.02502.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Babadi, B., and Brown, E. N. (2014). A review of multitaper spectral analysis. IEEE Trans. Biomed. Eng. 61, 1555–1564. doi: 10.1109/TBME.2014.2311996

CrossRef Full Text | Google Scholar

Beenhakker, M. P., and Huguenard, J. R. (2009). Neurons that fire together also conspire together: is normal sleep circuitry hijacked to generate epilepsy? Neuron 62, 612–632. doi: 10.1016/j.neuron.2009.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Brennan, G. P., Baram, T. Z., and Poolos, N. P. (2016). Hyperpolarization-activated cyclic nucleotide-gated (HCN) channels in epilepsy. Cold Spring Harbor Perspect. Med. 6:a022384. doi: 10.1101/cshperspect.a022384

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciumas, C., Saignavongs, M., Ilski, F., Herbillon, V., Laurent, A., Lothe, A., et al. (2014). White matter development in children with benign childhood epilepsy with centro-temporal spikes. Brain 137, 1095–1106. doi: 10.1093/brain/awu039

PubMed Abstract | CrossRef Full Text | Google Scholar

Clementeperez, A., Makinson, S. R., Higashikubo, B., Brovarney, S., Cho, F. S., Urry, A., et al. (2017). Distinct thalamic reticular cell types differentially modulate normal and pathological cortical rhythms. Cell Rep. 19, 2130–2142. doi: 10.1016/j.celrep.2017.05.044

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., Born, J., Claussen, J. C., and Martinetz, T. (2016a). Modeling the effect of sleep regulation on a neural mass model. J. Comput. Neurosci. 41, 15–28. doi: 10.1007/s10827-016-0602-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. S., Weigenand, A., Ngo, H. V., Marshall, L., Born, J., Martinetz, T., et al. (2016b). 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Csillik, B., Mihaly, A., Krisztinpeva, B., Chadaide, Z., Samsam, M., Knyiharcsillik, E., et al. (2005). Gabaergic parvalbumin-immunoreactive large calyciform presynaptic complexes in the reticular nucleus of the rat thalamus. J. Chem. Neuroanat. 30, 17–26. doi: 10.1016/j.jchemneu.2005.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Cueni, L., Canepari, M., Luján, R., Emmenegger, Y., Watanabe, M., Bond, C. T., et al. (2008). T-type ca 2+ channels, sk2 channels and sercas gate sleep-related oscillations in thalamic dendrites. Nat. Neurosci. 11, 683–692. doi: 10.1038/nn.2124

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

CrossRef Full Text | Google Scholar

Difrancesco, J. C., and Difrancesco, D. (2015). Dysfunctional hcn ion channels in neurological diseases. Front. Cell. Neurosci. 9:71. doi: 10.3389/fncel.2015.00071

PubMed Abstract | CrossRef Full Text | Google Scholar

Dryżałowski, P., Jóźwiak, S., Franckiewicz, M., and Strzelecka, J. (2018). Benign epilepsy with centrotemporal spikes-current concepts of diagnosis and treatment. Neurologia i Neurochirurgia Polska 52, 677–689. doi: 10.1016/j.pjnns.2018.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, D., Liao, F., and Wang, Q. (2017). The pacemaker role of thalamic reticular nucleus in controlling spike-wave discharges and spindles. Chaos 27:073103. doi: 10.1063/1.4991869

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernández, I. S., Peters, J. M., Akhondi-Asl, A., Klehm, J., Warfield, S. K., and Loddenkemper, T. (2017). Reduced thalamic volume in patients with electrical status epilepticus in sleep. Epilepsy Res. 130, 74–80. doi: 10.1016/j.eplepsyres.2017.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernández, I. S., Takeoka, M., Tas, E., Peters, J., Prabhu, S., Stannard, K., et al. (2012). Early thalamic lesions in patients with sleep-potentiated epileptiform activity. Neurology 78, 1721–1727. doi: 10.1212/WNL.0b013e3182582ff8

PubMed Abstract | CrossRef Full Text | Google Scholar

Francis, J. T., Xu, S., and Chapin, J. K. (2008). Proprioceptive and cutaneous representations in the rat ventral posterolateral thalamus. J. Neurophysiol. 99, 2291–2304. doi: 10.1152/jn.01206.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gentet, L. J., and Ulrich, D. (2003). Strong, reliable and precise synaptic connections between thalamic relay cells and neurones of the nucleus reticularis in juvenile rats. J. Physiol. 546, 801–811. doi: 10.1113/jphysiol.2002.032730

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodfellow, M., Schindler, K., and Baier, G. (2011). Intermittent spike-wave dynamics in a heterogeneous, spatially extended neural mass model. Neuroimage 55, 920–932. doi: 10.1016/j.neuroimage.2010.12.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartoyo, A., Cadusch, P. J., Liley, D. T. J., and Hicks, D. G. (2019). Parameter estimation and identifiability in a neural population model for electro-cortical activity. PLoS Comput. Biol. 15:e1006694. doi: 10.1371/journal.pcbi.1006694

PubMed Abstract | CrossRef Full Text | Google Scholar

Helling, R. M., Koppert, M. M. J., Visser, G. H., and Kalitzin, S. (2015). Gap junctions as common cause of high-frequency oscillations and epileptic seizures in a computational cascade of neuronal mass and compartmental modeling. Int. J. Neural Syst. 25:1550021. doi: 10.1142/S0129065715500215

PubMed Abstract | CrossRef Full Text | Google Scholar

Hocepied, G., Legros, B., Van Bogaert, P., Grenez, F., and Nonclercq, A. (2013). Early detection of epileptic seizures based on parameter identification of neural mass model. Comput. Biol. Med. 43, 1773–1782. doi: 10.1016/j.compbiomed.2013.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Iavarone, E., Yi, J., Shi, Y., Zandt, B., Oreilly, C., Van Geit, W., et al. (2019). Experimentally-constrained biophysical models of tonic and burst firing modes in thalamocortical neurons. PLoS Comput. Biol. 15:e1006753. doi: 10.1371/journal.pcbi.1006753

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

Jones, E. G. (2007). The Thalamus, 2nd Edn. Cambridge: Cambridge University Press.

Google Scholar

Kellaway, P. (2000). The electroencephalographic features of benign centrotemporal (rolandic) epilepsy of childhood. Epilepsia 41, 1053–1056. doi: 10.1111/j.1528-1157.2000.tb00296.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Knox, A. T., Glauser, T., Tenney, J., Lytton, W. W., and Holland, K. (2018). Modeling pathogenesis and treatment response in childhood absence epilepsy. Epilepsia 59, 135–145. doi: 10.1111/epi.13962

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramer, M. A., Stoyell, S. M., Chinappen, D., Ostrowski, L. M., Spencer, E. R., Morgan, A. K., et al. (2021). Focal sleep spindle deficits reveal focal thalamocortical dysfunction and predict cognitive deficits in childhood epilepsy with centrotemporal spikes. J. Neurosci. 41, 1816–1829. doi: 10.1523/JNEUROSCI.2009-20.2020

CrossRef Full Text | Google Scholar

Lal, D., Reinthaler, E. M., Schubert, J., Muhle, H., Riesch, E., Kluger, G., et al. (2014). Depdc5 mutations in genetic focal epilepsies of childhood. Ann. Neurol. 75, 788–792. doi: 10.1002/ana.24127

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y., Hwang, S. K., and Kwon, S. (2017). The clinical spectrum of benign epilepsy with centro-temporal spikes: a challenge in categorization and predictability. J. Epilepsy Res. 7, 1–6. doi: 10.14581/jer.17001

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemke, J. R., Lal, D., Reinthaler, E. M., Steiner, I., Nothnagel, M., Alber, M., et al. (2013). Mutations in grin2a cause idiopathic focal epilepsy with rolandic spikes. Nat. Genet. 45, 1067–1072. doi: 10.1038/ng.2728

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., and Cleland, T. A. (2017). A coupled-oscillator model of olfactory bulb gamma oscillations. PLoS Comput. Biol. 13:e1005760. doi: 10.1371/journal.pcbi.1005760

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Ji, G., Yu, Y., Yu, Y., Ding, M., Tang, Y., et al. (2017). Epileptic discharge related functional connectivity within and between networks in benign epilepsy with centrotemporal spikes. Int. J. Neural Syst. 27:1750018. doi: 10.1142/S0129065717500186

PubMed Abstract | CrossRef Full Text | Google Scholar

Marten, F., Rodrigues, S., Benjamin, O., Richardson, M. P., and Terry, J. R. (2009). Onset of polyspike complexes in a mean-field model of human electroencephalography and its application to absence epilepsy. Philos. Trans. R. Soc. A 367, 1145–1161. doi: 10.1098/rsta.2008.0255

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirandola, L., Cantalupo, G., Vaudano, A. E., Avanzini, P., Ruggieri, A., Pisani, F., et al. (2013). Centrotemporal spikes during nrem sleep: the promoting action of thalamus revealed by simultaneous EEG and fMRI coregistration. Epilepsy Behav. Case Rep. 1, 106–109. doi: 10.1016/j.ebcr.2013.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Moussaoui, S., Carteret, C., Brie, D., and Mohammaddjafari, A. (2006). Bayesian analysis of spectral mixture data using markov chain monte carlo methods. Chemometr. Intell. Lab. Syst. 81, 137–148. doi: 10.1016/j.chemolab.2005.11.004

CrossRef Full Text | Google Scholar

Ostrowski, L. M., Song, D. Y., Thorn, E. L., Ross, E. E., Stoyell, S. M., Chinappen, D. M., et al. (2019). Dysmature superficial white matter microstructure in developmental focal epilepsy. Brain Commun. 1:fcz002. doi: 10.1093/braincomms/fcz002

PubMed Abstract | CrossRef Full Text | Google Scholar

Paoletti, P. (2011). Molecular basis of nmda receptor functional diversity. Eur. J. Neurosci. 33, 1351–1365. doi: 10.1111/j.1460-9568.2011.07628.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavlou, E., Gkampeta, A., Evangeliou, A., and Athanasiadoupiperopoulou, F. (2012). Benign epilepsy with centro-temporal spikes (BECTS): relationship between unilateral or bilateral localization of interictal stereotyped focal spikes on EEG and the effectiveness of anti-epileptic medication. Hippokratia 16, 221–224.

PubMed Abstract | Google Scholar

Poolos, N. P. (2004). The Yin and Yang of the h-channel and its role in epilepsy. Epilepsy Curr. 4, 3–6. doi: 10.1111/j.1535-7597.2004.04101.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Reif, P. S., Tsai, M., Helbig, I., Rosenow, F., and Klein, K. M. (2017). Precision medicine in genetic epilepsies: break of dawn? Expert Rev. Neurotherapeut. 17, 381–392. doi: 10.1080/14737175.2017.1253476

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritter-Makinson, S., Clemente-Perez, A., Higashikubo, B., Cho, F. S., Holden, S. S., Bennett, E., et al. (2019). Augmented reticular thalamic bursting and seizures in Scn1a-dravet syndrome. Cell Rep. 26, 54–64. doi: 10.1016/j.celrep.2018.12.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosler, 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

PubMed Abstract | CrossRef Full Text | Google Scholar

Ross, E. E., Stoyell, S. M., Kramer, M. A., Berg, A. T., and Chu, C. J. (2020). The natural history of seizures and neuropsychiatric symptoms in childhood epilepsy with centrotemporal spikes (CECTS). Epilepsy Behav. 103:106437. doi: 10.1016/j.yebeh.2019.07.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Sánchez Fernández, I., Chapman, K. E., Peters, J. M., Harini, C., Rotenberg, A., and Loddenkemper, T. (2013). Continuous spikes and waves during sleep: electroclinical presentation and suggestions for management. Epilepsy Res. Treat. 2013:583531. doi: 10.1155/2013/583531

PubMed Abstract | CrossRef Full Text

Scheuer, M. L., Bagic, A., and Wilson, S. B. (2017). Spike detection: inter-reader agreement and a statistical turing test on a large data set. Clin. Neurophysiol. 128, 243–250. doi: 10.1016/j.clinph.2016.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Steynross, D. A., Steynross, M. L., Sleigh, J. W., Wilson, M. T., Gillies, I. P., and Wright, J. J. (2005). The sleep cycle modelled as a cortical phase transition. J. Biol. Phys. 31, 547–569. doi: 10.1007/s10867-005-1285-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Suffczynski, P., Kalitzin, S., and Silva, F. H. L. D. (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

Thomson, D. J., and Haley, C. L. (2014). Spacing and shape of random peaks in non-parametric spectrum estimates. Proc. R. Soc. A Math. Phys. Eng. Sci. 470, 20140101–20140101. doi: 10.1098/rspa.2014.0101

PubMed Abstract | CrossRef Full Text | Google Scholar

Thorn, E. L., Ostrowski, L. M., Chinappen, D. M., Jing, J., Westover, M. B., Stufflebeam, S. M., et al. (2020). Persistent abnormalities in rolandic thalamocortical white matter circuits in childhood epilepsy with centrotemporal spikes. Epilepsia. 61, 2500–2508. doi: 10.1111/epi.16681

PubMed Abstract | CrossRef Full Text | Google Scholar

Traub, R. D., Wong, R. K., Miles, R., and Michelson, H. (1991). A model of a ca3 hippocampal pyramidal neuron incorporating voltage-clamp data on intrinsic conductances. J. Neurophysiol. 66, 635–650. doi: 10.1152/jn.1991.66.2.635

PubMed Abstract | CrossRef Full Text | Google Scholar

Wamsley, E. J., Tucker, M. A., Shinn, A. K., Ono, K. E., Mckinley, S. K., Ely, A. V., et al. (2012). Reduced sleep spindles and spindle coherence in schizophrenia: mechanisms of impaired memory consolidation? Biol. Psychiatry 71, 154–161. doi: 10.1016/j.biopsych.2011.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Weigenand, A., Costa, M. S., Ngo, H. 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

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Wimmer, R. D., Astori, S., Bond, C. T., Rovó, Z., Chatton, J.-Y., Adelman, J. P., et al. (2012). Sustaining sleep spindles through enhanced sk2-channel activity consolidates sleep and elevates arousal threshold. J. Neurosci. 32, 13917–13928. doi: 10.1523/JNEUROSCI.2313-12.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, W., and Zhou, D. (2017). Progress in unraveling the genetic etiology of rolandic epilepsy. Seizure 47, 99–104. doi: 10.1016/j.seizure.2017.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., and Luo, J. (2018). Mutations of n-methyl-d-aspartate receptor subunits in epilepsy. Neurosci. Bull. 34, 549–565. doi: 10.1007/s12264-017-0191-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zador, A. M., Koch, C., and Brown, T. H. (1990). Biophysical model of a hebbian synapse. Proc. Natl. Acad. Sci. U.S.A. 87, 6718–6722. doi: 10.1073/pnas.87.17.6718

CrossRef Full Text | Google Scholar

Zhao, X., and Robinson, P. A. (2015). Generalized seizures in a neural field model with bursting dynamics. J. Comput. Neurosci. 39, 197–216. doi: 10.1007/s10827-015-0571-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Żygierewicz, J., Suffczyński, P., and Blinowska, K. (2001). A model of sleep spindles generation. Neurocomputing 38, 1619–1625. doi: 10.1016/S0925-2312(01)00516-1

CrossRef Full Text | Google Scholar