Abstract
Absence seizures are characterized by brief interruptions of conscious experience accompanied by oscillations of activity synchronized across many brain areas. Although the dynamics of the thalamocortical circuits are traditionally thought to underlie absence seizures, converging experimental evidence supports the key involvement of the basal ganglia (BG). In this theoretical work, we argue that the BG are essential for the maintenance of absence seizures. To this end, we combine analytical calculations with numerical simulations to investigate a computational model of the BG-thalamo-cortical network. We demonstrate that abnormally strong striatal feedforward inhibition can promote synchronous oscillatory activity that persists in the network over several tens of seconds as observed during seizures. We show that these maintained oscillations result from an interplay between the negative feedback through the cortico-subthalamo-nigral pathway and the striatal feedforward inhibition. The negative feedback promotes epileptic oscillations whereas the striatal feedforward inhibition suppresses the positive feedback provided by the cortico-striato-nigral pathway. Our theory is consistent with experimental evidence regarding the influence of BG on seizures (e.g., with the fact that a pharmacological blockade of the subthalamo-nigral pathway suppresses seizures). It also accounts for the observed strong suppression of the striatal output during seizures. Our theory predicts that well-timed transient excitatory inputs to the cortex advance the termination of absence seizures. In contrast with the thalamocortical theory, it also predicts that reducing the synaptic transmission along the cortico-subthalamo-nigral pathway while keeping constant the average firing rate of substantia nigra pars reticulata reduces the incidence of seizures.
SIGNIFICANCE STATEMENT Absence seizures are characterized by brief interruptions of consciousness accompanied by abnormal brain oscillations persisting tens of seconds. Thalamocortical circuits are traditionally thought to underlie absence seizures. However, recent experiments have highlighted the key role of the basal ganglia (BG). This work argues for a novel theory according to which the BG drive the oscillatory patterns of activity occurring during the seizures. It demonstrates that abnormally strong striatal feedforward inhibition promotes synchronous oscillatory activity in the BG-thalamo-cortical network and relate this property to the observed strong suppression of the striatal output during seizures. The theory is compatible with virtually all known experimental results, and it predicts that well-timed transient excitatory inputs to the cortex advance the termination of absence seizures.
Introduction
Absence epileptic seizures are characterized by brief interruptions of conscious experience accompanied by abnormal brain oscillatory activity (2.5–4 Hz) (Crunelli and Leresche, 2002). These oscillations are synchronous over many brain areas. They can be recorded in electroencephalograms (EEGs) as spike-and-wave discharges (SWDs) (Gibbs et al., 1935). In animal models of absence epilepsy, SWDs are abolished by lesions or inactivations of the cortex or the thalamus (Avoli and Gloor, 1982; Vergnes and Marescaux, 1992; Meeren et al., 2005, 2009; Polack et al., 2009). These experimental results led to the thalamocortical theory that the interactions between the thalamus and the cortex generate absence seizures (Prince and Farrell, 1969; Destexhe, 1998, 2007; Avoli, 2012). Although there is converging evidence that absence seizures are triggered from a specific cortical focus (Meeren et al., 2002, 2005; Polack et al., 2007, 2009), to the best of our knowledge, there is no proof that the thalamus and cortex are sufficient to maintain seizures (Danober et al., 1998; Depaulis et al., 2016). Polack et al. (2009), for instance, recently showed that inactivating the thalamic region reciprocally connected to the initiating cortical area did not suppress SWD. Thus, it remains unclear which network maintains seizures that can last tens of seconds.
Electrophysiological recordings in animal models have revealed synchronous oscillatory activity in the basal ganglia (BG) (Depaulis et al., 1994; Paz et al., 2009). Importantly, lesions or inactivations in the BG suppress seizures (Depaulis et al., 1994; Paz et al., 2007). Furthermore, fMRI studies have shown that the striatum is deactivated during seizures in both human patients (Moeller et al., 2008) and in the well-established rodent model of absence epilepsy, the Genetic Absence Epilepsy Rat from Strasbourg (GAERS) (David et al., 2008). This is consistent with the suppression of the firing activity of striatal medium spiny neurons (MSNs) during seizures in GAERS, probably due to a strong feedforward inhibition by striatal fast-spiking interneurons (FSIs) (Slaght et al., 2004; Paz et al., 2009). However, the precise role of the BG during SWD and, in particular, the implications of MSN firing suppression remain unclear.
This paper investigates the hypothesis that the maintenance of absence seizures over several tens of seconds emerges from the dynamics of the BG-thalamo-cortical network. Specifically, we hypothesize that absence seizures are an electroclinical symptom associated with bistability between two states of the dynamics of this network: a normal state in which the network activity is nonoscillatory and asynchronous, and an abnormal state in which the activity is oscillatory and spatially synchronized. We explore this hypothesis theoretically in computational models of the BG-thalamo-cortical network at different levels of physiological detail. This leads us to argue that, although absence seizures are initiated in the cortex, they are maintained by the recurrent dynamics of the BG-thalamo-cortical loop in a way that depends crucially on striatal FSI. We conclude with some predictions that can be tested experimentally.
Materials and Methods
Our model of the GAERS BG-thalamo-cortical network is made up of neuronal populations representing the somatosensory cortex (Ctx), the ventromedial (VM) nucleus of the thalamus (Th), and the BG (see Fig. 1). The BG nuclei considered here are the dorsal striatum, the subthalamic nucleus (STN), the substantia nigra pars reticulata (SNr), and the globus pallidus pars externa (GPe). For simplicity, each of these nuclei was modeled as a single population of neurons (Leblois et al., 2006), except for the striatum for which we include two populations, one of which representing the MSN that project to the SNr and GPe, and the other representing the FSI that project to the MSN. We included the VM because it is a BG-receiving “nonspecific” thalamic nucleus projecting to layer I of almost all the cortical areas, including the somatosensory cortex (Herkenham, 1979), and it displays bursting activity during SWD (Paz et al., 2007). The cortex is also described as a single population of neurons, thus neglecting the complex corticocortical circuitry. In other words, we assumed that the global input to the cortex is simply forwarded to the global cortical output, and we ignored any dynamics that are faster or more local than the dynamics of the SWD. Thus, the transmission delay from the thalamus to cortical neurons in our model (Table 1) should be seen as the sum of the delay from the thalamocortical neurons to the thalamus-receiving cortical neurons and the delay from these cortical neurons to the cortico-subthalamic and cortico-striatal neurons.
Anatomical studies (Alexander et al., 1986; Kelly and Strick, 2004) indicate that the BG together with the thalamus and cortex constitute feedback loops (i.e., that the output from the cortex propagates back to the same portion of the cortex through the BG and thalamus). Taking into account the polarities of the synaptic connections that go through the cortex, the STN and the SNr, the hyperdirect pathway (cortex-STN-SNr) (Nambu et al., 2000) inhibits the thalamus and therefore effectively provides negative (dis-facilitatory) feedback to the cortex (hyperdirect feedback loop; see Fig. 1, blue). By contrast, the direct pathway, which includes the cortex, MSN population, and the SNr, provides a disinhibition of the thalamus (Chevalier and Deniau, 1990) and thus positive (disinhibition) feedback to the cortex (direct feedback loop; red). The third pathway, which comprises the cortex, the FSI, the MSN, and the SNr, involves three serial inhibitory connections and therefore provides negative feedback to the cortex (loop through FSI; green). The key idea in our theory is that absence seizures persist primarily as a result of the competition between these three loops. In particular, loops that include the GPe, such as the one involving the indirect pathway (cortex-striatum-GPe-STN-SNr) or the STN-GPe-STN loop do not play a significant role. In our model, each major component is described by only one population of neurons, except for the striatum, which comprises two populations, corresponding to the MSN and the FSI. Thus, the complex cortical and thalamic internal circuitry (e.g., the different cortical layers and the reciprocal interactions between nucleus reticularis thalami [nRT] and thalamocortical neurons) are not included. This model is sufficient to investigate the basic role of the feedback loops in this circuit in the emergence of bistable dynamics.
We use two levels of description for the neuronal dynamics: namely, minimal rate-based dynamics (hereafter, rate model) in which the activity of each population is described by a single rate variable (Wilson and Cowan, 1972; Dayan and Abbott, 2001); and an integrate-and-fire model (hereafter, spiking model) describing the dynamics of the membrane potential of each neuron in each population (Lapicque, 1907; Dayan and Abbott, 2001). The parameters listed in Tables 1 and 2 are used as the reference set.
The rate model
In this model, the activity of population p = MSN, FSI, STN, SNr, Th, Ctx, is described by a variable, mp, the dynamics of which are governed by the following: Here, τp is the time constant of the dynamics, Gp(Ip(t)) the input–output relationship of population p and Ip(t) the total input current in that population, at time t: where ∑q indicates summation over presynaptic populations, Jpq is a synaptic weight, Δq is a delay, and hp a constant background input to population p.
We model the synaptic inputs as current entries. Specifically, we treat glutamatergic and GABAergic postsynaptic responses in a simplified manner as positive (excitatory) and negative (inhibitory) input currents, respectively. For simplicity, we take a threshold-linear input–output relationship, Gp(Ip) = [Ip]+ where [x]+ = x if x > 0 otherwise 0, for all populations except FSI. For the latter, we take an expanding nonlinear function as follows: where c is a constant controlling the level of nonlinearity (see Fig. 5B). Our main result does not depend on the value of c or on the exact shape of GFSI. However, the mathematical proof is beyond the scope of this paper. The quantitative dependency of the bistability on a similar nonlinearity parameter η of the spiking model is described in Results. Variables and constants mp(t), Ip(t) and hp are all scaled to the unit of the firing rate Gp(Ip(t)) and the synaptic weight Jpq is a dimensionless parameter.
This model has many parameters. We constrained our exploration of the model properties to a combination of parameters leading to a physiologically reasonable firing rate in each population when the network is in a normal state (Table 1). These firing rates should be compared with the level of activities reported in Slaght et al. (2004) and Polack et al. (2007) for the cortex, Slaght et al. (2004) for the MSN, Mallet et al. (2005) for the FSI, Paz et al. (2005) for the STN, and Deransart et al. (2003) for the SNr.
The spiking model
Our spiking model consists of the six or seven aforementioned populations (with or without GPe, respectively). In each population, single cells are modeled as integrate-and-fire neurons. The subthreshold dynamics of the membrane potential, Vp,i, of a neuron i in population p are as follows:
where, τp is the membrane time constant. We use τp = 10ms for all populations, except MSN and FSI. For MSN, we take τMSN = 10 ms whereas for FSI, τFSI = 5ms, to reflect slow dynamics of MSN due to the powerful potassium inwardly rectifying current (Nisenbaum and Wilson, 1995) and the fast dynamics of FSI whose effect on MSN is effectively as fast as the direct cortical excitation (Mallet et al., 2006; Pidoux et al., 2011). Neurons receive Gaussian noisy input from outside the network, ξp,i(t), such that: 〈ξp,i(t)〉 = 0 and 〈ξp,i(t) ξp,i(t′)〉 =
The quantity Ip,i(t) represents the total synaptic input into the neuron from its interactions with other neurons in the BG-thalamo-cortical network modeled as a current entry. It is given by the following: where ∑k indicates the summation over all spikes of all neurons in population q presynaptic to neuron i, tk is the timing of kth spike, Kq is the average number of neurons in population q projecting to a neuron in a postsynaptic population, τd and τr are the synaptic decay and rise times, which were assumed to be identical for all the interactions, hp is a background current from outside the network. We also assumed that neurons in p and q are randomly connected with probability Nq/Kq where Nq is the number of neurons in population q. Jpq denotes the strength of the postsynaptic potential in a neuron from population p induced by an action potential of a neuron in population q.
When the membrane potential Vp,i reaches the threshold θ, it is reset to Vr, that is: For constant input, μ = Ip(t), the firing rate of the neuron can be calculated analytically (Ricciardi, 1977) yielding the following:
Analysis of the rate model
In our model, the existence of oscillatory and nonoscillatory dynamics representing the normal and seizure state, respectively, depends on the gain of the feedback loops through the BG. The gain of a connection from an upstream to a downstream population (e.g., cortex-STN) characterizes the extent to which the activity of the downstream population (e.g., STN) changes with the activity of the upstream population (e.g., cortex) when keeping the activities of all other populations in the network fixed. The value of this gain can be calculated as the ratio of the change induced in the downstream population to the change imposed in the upstream population. To define the gain of a closed loop (e.g., cortex-STN-SNr-thalamus-cortex), consider an “unfolded” open path by assuming that the initial and final populations are different (e.g., cortex1-STN-SNr-thalamus-cortex2). The gain of the loop is defined as the ratio of the change in firing rate of the final population (e.g., cortex2) to the change in the firing rate of the initial population (e.g., cortex1), whereas inputs from the populations outside of this pathway remain constant. Overall, the gain of a closed loop describes how a perturbation applied to the activity of a given population propagates back to it through the loop. In our model, this gain is the product of all synaptic weights and slopes of the input–output relationships of the populations (see Eq. 3 below).
We analytically investigated the stability of fixed point solutions to the rate dynamics. This depends on the feedback gains of the closed loop along the hyperdirect, the direct loop via MSN, and the direct loop via FSI, defined as follows: where J̃pq is the effective synaptic interaction strength defined by the following: and νp is the firing rate of population p in the fixed point (normal state) given in Table 1. J̃pq = Jpq for p ≠ FSI because Gp is threshold linear in this case.
To characterize the stability of a fixed point solution, the evolution of a small perturbation ϵ around this solution needs to be calculated. For this purpose, we substituted mp(t) = ϵpexp(λt) + νp in Equation 1 yielding the following: Here ϵ is the six dimensional vector of the perturbation with element ϵp. The elements of the 6 × 6 stability matrix M are as follows: The growth rates of the perturbations λ are then the solutions of the equation P(λ) = 0 where P(λ) = det M. The fixed point solution is stable if all the real part of all the solutions λ, the characteristic exponents of the fixed point, are negative. Therefore, a necessary condition for instability onset when some parameter changes is for one of the characteristic exponents to change sign. The parameter at this onset can be found by imposing the characteristic exponent to be zero (λ = iω), that is, by solving the following: where ω is real.
For instantaneous striatal feedforward inhibition, i.e., τFSI = 0 and ΔFSI = 0, the function P(iω) depends linearly on A and B + C. Equations 4 and 5 can then be solved for A and B + C as a function of ω, yielding the stability boundary of the fixed point as a curve parameterized by ω in the A − (B + C). On this curve, a Hopf bifurcation occurs. The unstable mode oscillates with a frequency ω determined by Equations 4 and 5. For given (A − C) off the Hopf bifurcation boundary, we approximate the frequency ω by the nearest point, ω*, on the boundary, namely: where Aω, Bω, Cω are the solutions to Equations 4 and 5 with a given ω. This estimate is in good agreement with the numerical simulations as shown in Figure 5D.
In Figure 3A, we use B + C instead of B as the y-axis of the bifurcation diagram to align the Hopf bifurcation boundaries with different combinations of B and C. This is possible because the Hopf bifurcation depends solely on the value of A and B + C, as explained above. Without the shift by C, the bistable region without FSI (see Fig. 3A, black striped region; A − (B + C)) would be below the bistable region with FSI (see Fig. 5A, black region; A − B).
Simulations
In all our simulations of the rate model, differential equations were integrated using the Euler scheme (Dayan and Abbott, 2001). Simulations of the spiking model were performed using the second-order Runge-Kutta scheme for stochastic differential equations (Honeycutt, 1992; Hansel et al., 1998). The integration time step was 0.1 ms for all simulations. The parameters used in the simulations are those given in Tables 1 and 2 unless specified otherwise.
To determine the domain of parameters in which the network exhibits bistability between a fixed point and oscillations, we proceeded as follows. We simulated the network with JSTN Ctx = 0 and other parameters as in Table 1 or Table 2. In this case, the dynamics converge to a fixed point because the hyperdirect loop is not effective. We then gradually increased |JSTN Ctx| until the network settled into an oscillatory state. Then we gradually reduced |JSTN Ctx| back to 0. The network dynamics are bistable for some value of |JSTN Ctx| if for that value the network is at a fixed point when |JSTN Ctx| increases but exhibits oscillations when |JSTN Ctx| decreases (hysteresis). To draw the bifurcation diagrams (see Figs. 3A, 6), we ran these simulations with different values of JMSN Ctx.
To get a reliable estimate of the bistability range, |JSTN Ctx| must be varied sufficiently slowly so that at each step the network is always in the stationary state. To determine whether the network state was oscillatory for a given value of the parameters, we computed an “oscillation index.” For the rate model, this index is the amplitude of the input to the cortical population. If it is larger than some threshold, O, the state is considered oscillatory. The smaller the threshold, the better the characterization of the network state, but the longer it takes. We determined that O = 5 was a good tradeoff between precision and simulation time. For the spiking model, we first computed the maximum autocorrelation of the population averaged cortical input over a time duration larger than twice ΔA, the sum of all the delays along the hyperdirect feedback loop. This value was normalized to the autocorrelation at the zero time delay so that it was always between 0 and 1. The threshold O = 0.35 is used to detect oscillations.
The spiking model dynamics are noisy, and a simple threshold crossing criterion to detect oscillations may not be precise enough to estimate the extent of the bistable region. To improve this estimate, we used the knowledge of the structure of the bifurcation diagram provided by the study of the rate model. For each simulation at a fixed value of JMSN Ctx, we applied three postprocessing algorithms: (1) We looked for the largest asynchronous state region. (2) We considered the locations detected as oscillatory or bistable to the left of the largest asynchronous state region in the bifurcation diagram as indeed misdetections. We thus classified them as being in the asynchronous region. (3) We considered the right-most transition line from the nonoscillatory to the oscillatory regimen as the boundary of the oscillatory region.
Results
The BG-thalamo-cortical network model
The computational model investigated in this work represents the cortex-basal ganglia thalamocortical circuit in the GAERS. It consists of several populations of neurons representing the pyramidal neurons in the somatosensory cortex, the thalamocortical neurons in the ventromedial thalamic nucleus, striatal MSN, striatal FSI, the subthalamic nucleus, the substantia nigra pars reticulate, and optionally the globus pallidus pars externa (Fig. 1). By reducing the thalamus and the cortex to single populations, we aimed to emphasize the predominant role of feedback loops through the BG. Details of the model and their justifications are presented in Materials and Methods. We first consider a “rate version” of our model (Wilson and Cowan, 1972; Dayan and Abbott, 2001) in which the activity in each population is represented by one global dynamical “rate variable.” In the second version of the model, each population consists of a large number of leaky integrate-and-fire spiking neurons (Lapicque, 1907; Dayan and Abbott, 2001).
Our theory focuses on how persistent oscillations arise from the competition between the hyperdirect (cortex-STN-SNr) (Nambu et al., 2000) and the direct (Cortex-striatum-SNr) pathways (Leblois et al., 2006). Therefore we first investigate in what ways the dynamics depend on the gain of these three loops (for definitions, see Materials and Methods) without including the indirect pathway in the network. This reduces the number of parameters in the model and simplifies the analysis. We then verify that the inclusion of the GPe does not qualitatively affect our results.
Bistability of BG-thalamo-cortical network dynamics
Our key hypothesis is that absence seizures are maintained for lengthy durations because the BG-thalamo-cortical network dynamics are bistable. We therefore begin by demonstrating that our spiking network model can indeed exhibit bistability and that this depends on the strength of the feedforward inhibition in the striatum.
An example of bistability in the simulated dynamics of our spiking model is depicted in Figure 2A. At the beginning of the simulation, the network is in an asynchronous state in which neurons in all the populations fire irregularly and asynchronously. As a result, the population average firing rate of the cortical neurons is essentially constant (∼4 spikes/s; not to be confused with the internal oscillation frequency of SWD) until a transient excitation (amplitude: 2 mV; duration: 8 ms) is applied to the cortex at t = 1 s (Fig. 2A, top). This transient stimulation induces a rapid change in the activity pattern of the network: it becomes synchronous and oscillates at a frequency of 10 Hz (oscillatory state). The network remains in this state until another transient excitatory stimulus (amplitude: 2 mV; duration: 50 ms) is applied to the cortex at t = 8 s. As a result of this stimulus, the network switches back to the asynchronous state it was in before the first stimulation. This behavior is typical of a network exhibiting bistable dynamics (Strogatz, 1994). Similarly, a transient excitatory input to the STN (Fig. 2F) or the thalamus (Fig. 2G) can induce transitions from the asynchronous state to the synchronous oscillatory network state and can also stop oscillations by switching the network back to the asynchronous state.
A notable feature of the network in the oscillatory state in Figure 2A (bottom) is the almost complete suppression of MSN firing. This abolition of MSN discharge is due to the rhythmic bursting activity of the striatal FSI (Fig. 2C–E), which inhibits the MSN. This can be verified by blocking the FSI inhibition on a single MSN, which causes this MSN to become strongly active and to fire periodic bursts of action potentials during the oscillations (Fig. 2E).
The strong feedforward inhibition of the MSN by the FSI plays a key role in the bistability of the dynamics in our BG-thalamo-cortical network model. Figure 2B shows that, if the strength of this inhibition is reduced by 80%, transient stimulations do not trigger a transition to an oscillatory state even if their amplitude is much stronger than the one used in Figure 2A. In such conditions, the network relaxes back very fast into the asynchronous state after the stimulus. This was observed for all the stimulus parameters (amplitude: 1–20 mV; duration: 1–100 ms) we tested, indicating that this reduction in the strength of striatal feedforward inhibition impedes the bistability of the network dynamics.
Unbalanced feedback loops in the BG-thalamo-cortical network sustain absence seizures
To understand the mechanism underlying the bistability in our BG-thalamo-cortical network model and the conditions in which it takes place, we first analyze how changing the interaction strength between the various populations affects the dynamics. To this end, it is convenient to consider the rate-based formulation of the model. The relative simplicity of this description allows us to determine its bifurcation diagram by combining analytical calculations with extensive numerical simulations.
The bifurcation diagram plotted in Figure 3A summarizes the ways in which the dynamics of the network depend on the strength of the feedback loops for the parameters given in Table 1. Figure 3B–D presents examples of these dynamics. The four regions in the bifurcation diagram correspond to four qualitatively different dynamical regimens: (1) If the direct and hyperdirect feedback gains (B + C and A) appropriately “balance,” the network always settles at a fixed point in which the activity of the various populations is constant in time (Fig. 3B). Moreover, following transient perturbations, the network eventually converges back to the stable fixed point, possibly with damped oscillations of activity. This is thus a monostable fixed point regimen (FP; Fig. 3A, middle region in white). In terms of our theory, it corresponds to the normal state. (2) When the hyperdirect feedback gain, A, is strong relative to the direct feedback gain, B, the network stationary dynamics are solely oscillatory, and the oscillations cannot be suppressed by transient perturbations (monostable oscillations, OSCs; Fig. 3A,C). These oscillations are sustained by the delayed negative feedback of the hyperdirect loop and therefore must be distinguished from transient oscillations, such as nonpathological beta oscillations related to movement (Leventhal et al., 2012; Feingold et al., 2015). (3) In between these two regimens, there is a region where the dynamics have two coexisting stable stationary states, namely, a fixed point and an oscillatory state (bistability, BST, Fig. 3A, gray). The network can switch between these two states if an appropriate transient excitatory input is applied to the cortex (Fig. 3D). We posit that this regimen underlies absence epilepsy. The fixed point corresponds to the normal state, whereas synchronized oscillations in the whole network correspond to the seizures. The fact that the oscillatory state is stable when experiencing sufficiently small perturbations corresponds to the persistence of the absence seizure oscillations (Fig. 3D). (4) If the positive feedback is overwhelming, the dynamics are unstable (U in Fig. 3A). In this case, the firing rate of the cortical neurons diverges. This nonphysiological regimen is not discussed here.
Strong striatal feedforward inhibition promotes absence epilepsy
The BG-thalamo-cortical model network can have bistable dynamics even if the feedforward inhibition of the MSN by the FSI is blocked (Fig. 3A; C = 0). However, in this case, the corresponding region in the bifurcation diagram is small (Fig. 3A, black striped region). A key result of our study is that increasing this feedforward inhibition results in an expansion of the bistable region (Fig. 3A, gray region), and thus promotes the abnormal epileptic regimen. In other words, a network in the monostable fixed point regimen (i.e., the healthy condition; Fig. 4A) can become bistable (i.e., epileptic condition; Fig. 4B) if the striatal feedforward inhibition is sufficiently strong. This is because the boundary between regions U and FP and the boundary between regions FP and OSC depend solely on the pair (A, B + C) (the ratio of B to C is irrelevant) whereas the boundary between regions BST and FP depends directly on C. The points in Figure 3A corresponding to Figure 4, A and B, are in the same location in the bifurcation diagram (i.e., they correspond to the same values of A and B + C whereas the value of B and C are different).
When the striatal feedforward inhibition is blocked, the MSN population is highly active during oscillations in the bistable regimen (Fig. 3A, striped region). Increasing striatal feedforward inhibition not only enlarges the bistable region but also reduces MSN activity during the oscillations. In the example depicted in Figure 4C, we chose the cortical excitation to the MSN to have the smallest value compatible with bistability (location marked with + in the bifurcation diagram Fig. 3A). Even in this setting, the activity of the MSN during the oscillations is larger than at the fixed point (by a factor of 1.2). More generally, if the cortical excitation is strong relative to the striatal feedforward inhibition (Fig. 3A, top right area of the gray region), MSNs are highly active during oscillations (Fig. 4D). Conversely, if the striatal feedforward inhibition is strong compared with the excitation coming directly from the cortex (Fig. 4B), the MSN activity can be suppressed during oscillations. However, this also requires sufficiently fast kinetics of this inhibition. If they are too slow (Fig. 4E; τFSI = 5ms; in other panels in Fig. 4; τFSI = 0), the MSN activity cannot be suppressed during the oscillations because the inhibition affects the MSN too late to cancel the excitation.
The mechanisms underlying bistability and the suppression of MSN activity
For the sake of analytical tractability, all the populations in our rate model have a threshold linear input–output relationship, that is, G(I) = I if I > 0 otherwise 0, except for FSI for which the input–output relationship exhibits an expanding nonlinearity (Eq. 2). We now argue that the latter nonlinearity has an important influence on the size of the region in the parameter space where the network dynamics are bistable and on the MSN activity in this regimen.
The expansion of the bistable region upon the increase of the striatal feedforward inhibition (Fig. 3A) requires nonlinearity in the FSI input–output relationship. Indeed, if this relationship is linear, the bifurcation diagram plotted as a function of the gains A and B + C (as in Fig. 3A) does not depend on the strength of the feedforward inhibition of the MSN (Fig. 5A). This is because, as far as the net input to the MSN is concerned, the network is virtually equivalent to one with no striatal feedforward inhibition.
The net input, IMSN, received by the MSN is the sum of the excitation coming directly from the cortical population and the effectively inhibitory disynaptic input coming from the cortex through the FSI. Because of the nonlinearity in the input–output relationship of the FSI (Fig. 5B), this net input can be expressed as a function of the cortical activity, mCtx (Fig. 5C). The combination of the nonlinearity in the input–output relationship of FSI and the feedforward inhibition makes the net input to MSN a nonmonotonic function of the cortical activity (Fig. 5C, solid line). When the cortex is nearly inactive (mCtx ∼ 0 Hz) MSNs are not active. On the other hand, when the cortical activity is higher than the normal state (mCtx ∼ 10 Hz), MSNs are also not active because the striatal feedforward inhibition overcomes the direct excitation from the cortex. Upon blockade of the striatal feedforward inhibition, the MSN would be highly active in this case (Fig. 5C, dashed line).
The fact that the bistable region becomes larger with the increase in striatal feedforward inhibition only requires the input–output relationship, GFSI (x), of FSI to exhibit an expanding nonlinearity (positive second-order derivative, GFSI″(x) > 0 for x > 0). It is thus a general feature of our network model. This can be shown mathematically if one assumes that the synaptic time constant is very small compared with the delay in the feedback loops τp ≪ Δ. Because of its technicality, the description of this proof is beyond the scope of this article. Of course, this enlargement depends quantitatively on the strength of the nonlinearity of GFSI(x) around the fixed point (i.e., how large the second-order derivative, GFSI″(x), is).
Because during oscillations MSNs are suppressed, and thus the activity cannot propagate along the direct pathway, the timescale of the oscillations is determined by the hyperdirect loop. Indeed, the internal frequency of the oscillations decreases monotonically when the delay along the hyperdirect loop is increased (Fig. 5D). When using the parameters in Table 1, the frequency is 9.6 Hz, a value that is close to the frequency of the oscillations reported in GAERS (Danober et al., 1998; Depaulis et al., 2016) (7–11 Hz). In our theory, the frequency of the oscillations depends on the effective delay along hyperdirect feedback loop (Fig. 5D). In particular, to account for the frequency of the oscillations of absence seizures in humans (∼3 Hz) requires an effective delay along the hyperdirect feedback loop of 150 ms, which is longer than in Table 1 (Fig. 5D). This contrasts with the thalamo-cortical mechanism of SWD in which the frequency of the oscillations depends primarily on the kinetics of intrathalamic inhibition (Pinault et al., 1998; Destexhe, 1999).
Asynchronous firing and synchronous oscillations in the BG-thalamo-cortical spiking network
In this section, we return to the dynamics of our integrate-and-fire model of the BG-thalamo-cortical network to show that the dynamical features we derived in the simple rate model also hold for this spiking neural network model. Building on the analysis of our rate model network, we determined the regions in the parameter space of this spiking network model where: (1) neurons fire in an essentially asynchronous manner, exhibiting only very weak correlations (ASYNC, as in the example shown in Fig. 2B); (2) the network develops oscillatory patterns of activity in which the neurons across the whole system fire with a substantial degree of synchrony (OSC); and (3) the network exhibits bistability between an asynchronous state and a synchronous oscillatory state (BST; Fig. 2A).
The bifurcation diagram of the spiking model is plotted in Figure 6A for the parameters given in Table 2. It is qualitatively similar to the one found in the rate model (Fig. 3A). When the positive feedback is sufficiently strong in the direct loop, the activity of the network is asynchronous (ASYNC). In this state, the population averaged activity is constant over time in each population. This corresponds to the fixed point regimen in the rate model. By contrast, for sufficiently strong negative feedback through the hyperdirect pathway, the asynchronous state becomes destabilized and the network settles in a state in which the neurons display synchronous oscillatory activity. In a broad region in between these two regimens, the dynamics are bistable. Depending on initial conditions, this activity can be asynchronous or oscillatory and synchronous. A transient input to the cortical population can switch the network activity from asynchronous spiking to synchronous oscillations (Fig. 2A).
As in the rate model, the bistable regimen depends crucially on the strength of the striatal feedforward inhibition. The range of parameters giving rise to bistability shrinks when the strength of this inhibition is reduced, and eventually vanishes in the example shown in Figure 6D. The mechanism enabling this strong control over bistability can be understood in terms of the result we have obtained in the rate model, which showed that transmission along the cortex-FSI-MSN pathway needs to be sufficiently nonlinear (Fig. 6E). To demonstrate this, we parametrized the synaptic strengths JFSI Ctx and JMSN FSI according to the following: where JFSI Ctx0 and JMSN FSI0 are fixed (Table 2). Varying parameter η changes the balance in the contributions of the cortex-FSI and FSI-MSN interactions without changing the overall gain of the cortex-FSI-MSN pathway. The bifurcation diagram plotted in Figure 6A corresponds to η = 1. Because the nonlinearity of this pathway stems from the nonlinearity of the input–output relationship of the FSI, it becomes more linear as η decreases. The nonlinearity of the pathway can come from short-term potentiation of cortex-FSI and/or FSI-MSN synapses. Our theory suggests that the source of nonlinearity is irrelevant to the global dynamics, and only the net effect on MSN (Figs. 5C, 6E) matters. Figure 6, B and C, plots the bifurcation diagrams of the network for two values of η. A comparison of these two bifurcation diagrams with the one plotted in Figure 6A indicates that the size of the bistable region reduces with η.
Finally, as was the case in our rate model, the FSI input–output relationship (Fig. 6F) in our spiking network is an expanding nonlinear function of the activity, mCtx, of the cortical population. Because of this expanding nonlinearity, the input to MSN, IMSN, is subthreshold at the minimum (∼0 spikes/s) and maximum (∼10 spikes/s) cortical population firing rate during the oscillations (Fig. 6E). This feature is essential for the suppression of MSN activity during the seizures as we found in the rate model.
Appropriately timed excitatory stimulation of the cortex terminates seizures
Figure 2 depicts an example of bistable dynamics in our spiking network model in which persistent oscillations can be both induced and terminated by a transient excitatory input to the cortical population. The fact that excitation can terminate the oscillations is a nontrivial feature of the mechanism underlying them. To clarify the conditions in which this can occur, we simulated our spiking network model with transient input to all cortical neurons while varying three parameters: duration, amplitude and phase (Fig. 7A). The phase φ is defined with respect to the oscillations in the cortical input and measured in units normalized from 0 to 1. A phase φ = 0.5 corresponds to the minimum of the oscillation in the activity. The combinations of parameters that yield a termination are depicted in Figure 7B. When this input has a short duration and a high amplitude, the network resets to the asynchronous state only if it occurs ∼φ ∼ 0.5 (Fig. 7C). On the other hand, if the transient input has a low amplitude and a sufficiently long duration such that it terminates just before φ ∼ 0, a reset always occurs (Fig. 7D). Upon the transient input, the MSN firing rate increases for a short period of time (Fig. 7C,D, star) regardless of the parameters. At the single-cell level, this can be seen as a rebound-like spiking activity (Figs. 2A, 7C,D, bottom). A similar rebound of activity was observed in GAERS striatal MSN at the end of SWD (Slaght et al., 2004) and was previously suggested as a possible termination mechanism for absence seizures (Paz et al., 2009). In the framework of our theory, maintenance of oscillatory activity during absence seizures results from the decreased efficiency of the direct feedback loop due to the decreased activity of MSN. Thus, activation of the MSN for a sufficiently long period of time can terminate the oscillations. Our theory therefore explains the link between the rebound of excitation in MSN and the end of SWD.
Bistable dynamics in the network model with the GPe included
According to our theory, absence seizures emerge primarily from the competition between the hyperdirect and direct feedback loops modulated by striatal feedforward inhibition. To show that the competition between these loops is sufficient to generate the seizures, the model investigated above only included these loops. In particular, we did not take into account the GPe and the feedback loops in which it is involved, namely, the indirect pathway (cortex-striatum-GPe-STN-SNr) or the STN-GPe-STN loop.
It should be noted that the dynamics of this simplified model exhibit bistability between asynchronous (normal) and oscillatory (pathological) activity patterns in a broad range of parameters. It is therefore to be expected that including the GPe should not qualitatively affect the existence of such a bistable regimen for the model network dynamics, at least if the added connections are not overly strong.
Because of the qualitative nature of our network model, proving that these “not overly strong” connections indeed cover a biologically plausible range is not as straightforward as just plugging in some estimated synaptic efficacies. To this end, we first determined for which values of JSTN GPe the GPe-STN subnetwork does not exhibit oscillations on its own, whereas the other parameters were fixed as in Figure 8A. Figure 8 shows that, for these values, the bifurcation diagram of the network model has the same structure (Fig. 8B, compare withFig. 6A) and similar neuronal dynamics (Fig. 8A,C,D; compare withFig. 2A,C,D) as in the network model considered above, which does not include the GPe. In particular, during seizures, the activity of the MSN exhibits subthreshold oscillations and the FSI display periodic bursting in both models.
Because MSNs are not active during the oscillations, the GPe affects the global dynamics mainly through the positive feedback cortex-STN-GPe-SNr-thalamus-cortex loop. As a result, the region in which the network exhibits bistability is shifted to the right in Figure 8B compared with Figure 6A. This is why, in the example of bistable behavior in Figure 8A,C,D, the values of JSTNCtx and JMSNCtx are different from those in Table 2. We conclude that the bistable dynamics (Fig. 8A,C,D) and the bifurcation structure (Fig. 8B) are qualitatively preserved under inclusion of the GPe.
Discussion
We showed that strong feedforward inhibition in the striatum promotes bistability in the dynamics of the BG-thalamo-cortical network. In the bistable regimen, activation of the cortex can trigger persistent oscillations. According to our theory, seizures in absence epilepsy correspond to such dynamics and cannot be maintained over tens of seconds by the cortex alone or the thalamocortical network without the BG. We thus propose that abnormally strong striatal feedforward inhibition is involved in the electroclinical symptoms of absence epilepsy.
Consistency with previous experimental results
Neurons in all BG nuclei exhibit strong oscillatory activity during SWD (Paz et al., 2009). Furthermore, in GAERS, the VM thalamic neurons, which are inhibited by the SNr and project to almost all cortical areas including the somatosensory cortex (Herkenham, 1979), display bursting activity during SWD (Paz et al., 2007). The nonspecificity of these projections fits with the generalized nature of the SWD. The interruption of SWD by the suppression of the glutamatergic input to the SNr coincides with the suppression of the bursting activity in VM neurons (Paz et al., 2007), thus emphasizing the role of BG-receiving thalamus in SWD. Pharmacological blockade of the subthalamo-nigral pathway prevents seizures (Depaulis et al., 1994; Deransart et al., 1998, 2000; Paz et al., 2007; Kase et al., 2012). Our theory naturally accounts for all these findings by highlighting the central role of BG in absence seizures. Striatal activity is strongly reduced during seizures, as reported in fMRI studies (David et al., 2008; Moeller et al., 2008) and electrophysiological recordings. During seizures, FSI increase their activity whereas MSN exhibit reduced spiking activity and subthreshold membrane oscillations synchronous with SWD (Slaght et al., 2004; Paz et al., 2009). Our theory links these features to the mechanisms of SWD. It also explains the almost complete inactivation of MSN during seizures as an outcome of a supralinear input–output relationship of FSI. This supralinear feedforward inhibition also leads to nonmonotonicity in the MSN response to cortical activation or sensory stimuli, as reported in anesthetized rats (Pidoux et al., 2011).
It has been argued that thalamic T-type channels and hence the GABAB receptors that deinactivate them are critical for SWD (Crunelli and Leresche, 1991; Avanzini et al., 1992, 1993; Hosford et al., 1992; Liu et al., 1992; Marescaux et al., 1992; von Krosigk et al., 1993; Bal et al., 1995; Tsakiridou et al., 1995; Danober et al., 1998; Kim et al., 2001). When VM neurons are in the bursting mode, their firing rate increases abruptly at threshold crossing (Sherman, 2006). This may be thought of in our framework as an increased gain of the input–output relationship of thalamic neurons close to their firing threshold. T-type channels, which allow thalamic bursting may therefore increase the amplitude of the overall gain (|A + B + C|) in our theory, yielding a shift toward the bistable regimen (Fig. 3A). Thus, the reported importance of T-type channels in the thalamus is also in line with our theory. It is also likely that the global gain of the VM population increases with the firing rates of VM neurons and is thus modulated by the strong inhibition from the nRT (Slaght et al., 2002; Pinault, 2004). Various animal models of absence epilepsy exhibit increased excitability in cortical neurons (Pumain et al., 1992; Luhmann et al., 1995; Avanzini et al., 1996; Meeren et al., 2005; D'Antuono et al., 2006; Polack et al., 2007; Avoli, 2012; Leresche et al., 2012). This prompted the idea that the thalamocortical loop maintains SWD (Meeren et al., 2005). However, an increase in the excitability of the cortex corresponds in our theory to an increase in the overall gain (|A + B + C|), hence shifting the network dynamics toward the pathological bistable regimen (Fig. 3A). The network can therefore switch from the nonepileptic to the epileptic regimen by changing the thalamocortical gain while keeping BG physiological properties unchanged. Finally, a transient perturbation in cortical activity can initiate and terminate seizures in our theory. This is consistent with reports that absence seizures in rats are initiated in the somatosensory cortex (Meeren et al., 2002, 2005; Polack et al., 2007, 2009) and with recent works showing that electrical stimulation of the cortex can initiate, shorten, or abort seizures (Berényi et al., 2012; Sadighi et al., 2013; van Heukelum et al., 2016).
Comparison with thalamocortical theory
The thalamus is involved in sleep spindles and was hypothesized to play a central role in SWD (Leresche et al., 2012). During sleep spindles, nRT neurons exhibit bursting mediated by T-type channels (McCormick and Bal, 1997). Manipulations of T-type channels also affect SWD and either pharmacologically suppressing nRT (Liu et al., 1991) or lesioning thalamic relay and reticular nuclei (Buzsáki et al., 1988; Vergnes and Marescaux, 1992; Avanzini et al., 1993; Meeren et al., 2009) decrease or abolish absence seizures. It has therefore been hypothesized that the cortex-nRT-TC network is responsible for sleep spindles and SWD (Gloor, 1978; Wang et al., 1995; Destexhe, 1998; Kostopoulos, 2000; Suffczynski et al., 2004; Marten et al., 2009; Taylor et al., 2013, 2014). In the Destexhe (1998, 1999) models, the nRT-TC network has only one stable state for a given set of parameters. Chen et al. (2014, 2015) investigated the ways in which the BG can modulate seizures in such a scenario by modifying the nigral inputs to the thalamus, bringing thalamo-cortical network inside or outside the regimen in which it oscillates. The models reported by Suffczynski et al. (2004), Marten et al. (2009), Taylor et al. (2013, 2014) rely on a bistability in the dynamics of the cortex-thalamus feedback loop. The manner in which the BG can modulate seizures was not investigated in this scenario. In these models, the nRT-TC network sustains oscillations and the BG modulate them (Depaulis et al., 1990, 1994; Chen et al., 2014, 2015). In contrast, our theory posits that the BG sustain oscillations and the nRT-TC network may modulate them by modifying the gain of BG-thalamo-cortical loop. It should be mentioned that recent studies have challenged the idea that SWD may represent “perverted” sleep spindles (Leresche et al., 2012). Indeed, lesions that leave the rostral part of the nRT intact suppress ipsilateral sleep spindles but increase SWD bilaterally (Meeren et al., 2009). Moreover, no differences have been reported between GAERS and nonepileptic rats regarding the properties of GABAergic synapses and neuronal excitability in nRT or TC nuclei (Danober et al., 1998; Slaght et al., 2002; Paz et al., 2007; Polack et al., 2009).
Limitations of the theory, perspectives, and predictions
The bistability in our model contrasts with the monostable oscillatory state reported by Leblois et al. (2006) to account for beta oscillations in the BG of parkinsonian patients and animals. This difference reflects the fact that SWDs and normal oscillation-free dynamics alternate in absence epilepsy, whereas the dynamics are always abnormal in Parkinson's disease (Hutchison et al., 2004; Gatev et al., 2006; Hammond et al., 2007). Our theory addresses the mechanism underlying the persistence of absence seizures but not the shape of spike-and-wave complexes during seizures. Chipaux et al. (2011) argued that fast recurrent activation of cortical GABAergic interneurons, which negatively control pyramidal neurons initiating SWD shapes these complexes. The rhythmic bursting activity of VM thalamocortical neurons during seizures may also contribute to shaping the sharp-spike pattern of spike-and-wave complexes through their excitatory impact on cortical neurons. VM thalamocortical neurons indeed exhibit brief bursts of action potentials followed by a transient period of hyperpolarization during SWD (Paz et al., 2007). This inhibition may be due to recurrent activation of nRT GABAergic neurons (Slaght et al., 2002; Pinault, 2004). Ictal and interictal durations are widely distributed even within individuals (Suffczynski et al., 2006). In our theory, such broad duration distributions may occur if spontaneous initiation and termination of seizures are induced by spatially coherent random fluctuations in neuronal activity in the BG-thalamocortical network. Our theory predicts that transient stimulation to the cortex can initiate but also terminate absence seizures. For termination, the stimulation has to satisfy a specific phase-duration-amplitude relationship (Fig. 7). In line with this prediction, transcranial (Berényi et al., 2012) or intracranial (Sadighi et al., 2013; van Heukelum et al., 2016) electrical stimulations and direct current injections (Perez Velazquez et al., 2007) can shorten or desynchronize SWD in epileptic rats, although the phase dependency has not been characterized. We have recently reported that electrically stimulating a cortical site (e.g., the somatosensory cortex) with a phase dependency as predicted by our model suppressed SWD in GAERS (n = 2) (Arakaki et al., 2013). These results require confirmation with more data. It should be noted, however, that this prediction is not an unequivocal test of our theory. A similar phase-duration-amplitude relationship may occur in a model in which the thalamocortical feedback loop alone underlies SWD (Suffczynski et al., 2004). A more specific test of our theory would be to suppress MSN activity in normal rats or the GABAergic inhibition within the striatum in GAERS using pharmacology or optogenetics. According to our theory, the former manipulation should promote SWD-like oscillations, whereas the latter should impede them. Another although experimentally more difficult test would consist of reducing the gain of the hyperdirect feedback loop without changing the gain of the direct loop or the mean firing rate in the SNr. This would be feasible by applying a combination of glutamatergic antagonist and agonist in the STN. Indeed, while a glutamatergic antagonist can block cortico-subthalamic transmission, the addition of a competitive agonist at the appropriate concentration would increase the average STN firing rate such that the mean interictal activity would remain unchanged. Our theory predicts that this type of manipulation should shorten the duration of the seizures and eventually suppress them if the cortico-subthalamic gain is sufficiently reduced.
Footnotes
This work was supported in part by Agence Nationale de la Recherche Grant SONGLEARN, 10-RPDoc-0016. T.A. was supported by Ecole des Neurosciences de Paris (Paris School of Neuroscience) doctoral fellowship. This work was performed in the framework of the France-Israel Laboratory of Neuroscience. We thank Ran Darshan for comments on an earlier version of the manuscript.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. David Hansel, Laboratory of Neurophysics, Physiology and Pathology, University Paris Descartes, 45 Rue des Saints-Pères, 75270 Paris, France. david.hansel{at}parisdescartes.fr