Our 475-cell thalamocortical network consists of single-compartment neuron models of the following types: layer 5 regular spiking (RS), layer 5 intrinsically bursting (IB), layer 6 non-tufted regular spiking (NRS) pyramidals and a deep low-threshold spiking (LTS) interneuron in the cortex, and thalamic reticular cells (RE) and thalamocortical cells (TC) in the thalamus. The membrane potential of each cell type is described by the following equations:

where V_i is the membrane potential for i = RS, IB, NRS, LTS, TC, RE, C_m1 = 0.9 and C_m2 = 1.0 is the specific membrane capacitance given in μF/cm². I_hold is the input current required to hold and maintain the voltage of the neuron (also used in this context for setting the resting membrane potential). A consistent set of units were maintained such that voltage is given in mV, ionic currents in mA/cm², membrane conductance densities in S/cm², and time in msec. Note that while input currents are given in nA, they are converted to a current density by dividing by the cell compartment's area.

All cortical cells have the same set of ionic currents with different kinetics in some instances. These consist of both transient and persistent sodium currents, potassium currents of the delayed rectifier, transient, slowly inactivating, and Ca²⁺-dependent types, in addition to both low- and high-threshold Ca²⁺ currents. The models for these cells were adapted from our previous work (Ahmed, 2019), which modified single-neuron models from Traub et al. (2005). The original model by Traub et al. featured multi-compartment neurons in a large-scale network, successfully replicating various thalamocortical phenomena. In our work, we simplified these complex neuronal models to single-compartment representations of the soma while preserving essential dynamics. This reduction required careful recalibration of membrane conductance densities to maintain physiologically realistic behavior. In contrast, the thalamic cells were based on the model by Destexhe et al. (1998), which describes thalamic neurons with fewer currents but still exhibits bursting behavior due to the presence of the low-threshold calcium current, with slower kinetics in the RE cells. All intrinsic currents follow the same formalism, described by the following equation:

where ḡ_ion represents the maximal conductance, and m and h represent the voltage gating of the ion channels. All synaptic currents are described by the following equations:

where ḡ_syn denotes the maximal synaptic conductance, and AMPA/GABA_A- and GABA_B- mediated synapses are described appropriately by the synaptic conductance function, s(t) which denotes the fraction of open synaptic channels at time t. The time of a presynaptic spike arrival is denoted by t₀, while α and β denote the forward (binding) and backward (unbinding) rate constants describing the binding of neurotransmitter, respectively. The duration of the neurotransmitter-mediated pulse is given by C_dur, and C_max refers to the maximal neurotransmitter concentration value. The details of the individual currents, synaptic dynamics, and model parameters are provided in Supplementary Sections 1, 2.

The network consists of the following cell populations: 25 layer 5 RS pyramidals, 75 layer 5 IB pyramidals, 75 layer 6 non-tufted RS pyramidals, 100 deep LTS interneurons, 100 thalamocortical cells, and 100 thalamic reticular cells. At baseline, a higher proportion of bursting-type pyramidal cells compared to regular-spiking cells in layer 5 was used to both introduce asymmetry in this layer, as well as model the frontal cortex. All cell layers were of equal size, except for layer 6, which had a reduced population size relative to layer 5 to reflect its typically lower density of pyramidal cells (DeFelipe et al., 2002). Each layer of cells is arranged in one dimension as shown in Figure 1 with synaptic connectivity varying between cell types as described. Cortical input into the thalamus is only from layer 6 of the cortex with excitatory synapses formed with ascending thalamic axons. All excitatory connections in the network are mediated by AMPA receptors, and inhibitory connections are mediated by either GABA_A only or a combination of GABA_A and GABA_B receptors.

Due to the unequal number of cells in some layers, the number of presynaptic neurons each postsynaptic neuron connects to (nPrePost) was layer-dependent, as given in Table 1. In each network layer, presynaptic cells were indexed relative to the indexing of postsynaptic cells. Within the thalamus and the cortex, each neuron connects to postsynaptic neurons within a radius of 11 indices, while between the thalamus and cortex, each presynaptic neuron connects with neurons within a radius of 21 postsynaptic indices. This cell count affects the weight of individual synapses which is defined by ḡ_syn/nPrePost where ḡ_syn is the maximal synaptic conductance for a particular synapse, as given in Table 2.

After allowing the network to reach steady state, a 1 nA stimulus was applied for 100 ms to a group of five layer 6 NRS neurons, inducing additional oscillations. Network behavior and states were characterized using these triggered oscillations. The model fitting parameters, given in Tables 2–4 were selected based on the method described in Supplementary Section 3.1, and were intended to be interpreted qualitatively. The same set of fitted parameters was used to model both the default and diseased states. In the default mode, the network exhibits spontaneous oscillations driven by spontaneous discharges in TC neurons, while in the diseased mode, the model exhibits synchronized SWDs due to alterations in the strength of cortical inhibition.

The model was implemented in NetPyNE, a Python package developed by Dura-Bernal et al. that allows for high-level specification of network connectivity of biological neuronal networks and simulations using the NEURON simulator (Dura-Bernal et al., 2019).

Mutations in genes encoding GABA_A receptor subunits, including the GABRG2 gene, have been implicated in CAE (Kananura et al., 2002; Kang and Macdonald, 2004; Ito et al., 2005; Marini et al., 2003). Functional studies of mutant receptors, particularly those containing the γ2(R82Q) and γ2(R43Q) subunits, have demonstrated reduced receptor surface expression in cortical pyramidal neurons and decreased GABA_A receptor current amplitudes (Kang and Macdonald, 2004; Macdonald et al., 2010). These alterations are consistent with a reduction in neuronal feedforward inhibition within cortical circuits (Bianchi et al., 2002; Currie et al., 2017). In our model, the effects of the GABRG2 mutation were implemented specifically in cortical pyramidal neurons across both layers by modifying the parameter corresponding to the maximum conductance of the GABA_A current (ḡ_max), as outlined below:

The effect of ALLO was modeled at the level of synapses, particularly GABA_A receptor-mediated synapses, by adjusting the ḡ_max, α, β, and C_dur parameters, as highlighted in Equations 9–10, 13.

Modifications to the highlighted parameters were informed by fitting our synapse model to experimental data, using the relative change between control GABA_A receptor activity and activity after ALLO application. It is important to note that the available experimental data vary considerably in methodology across studies, including differences in the species of central nervous system neurons used, and the concentrations of both GABA (0.005 to 0.1 mM) and ALLO (30 to 1,000 nM). While the physiological concentration of GABA released into the synaptic cleft ranges from 1 to 10 mM, the experimental data used to model the effects of ALLO are based on GABA concentrations that are at least an order of magnitude lower (Brickley and Mody, 2012; Overstreet et al., 2002; Tretter and Moss, 2008; Mozrzymas et al., 2003). Similarly, physiological ALLO levels in rat cortex typically range from 1 to 20 nM, yet experimental protocols generally apply concentrations that are one to two orders of magnitude higher (Pinna et al., 2000; Kelley et al., 2011). The lack of significant differences in response amplitude (i.e., relative change in ḡ_max) following the application of ALLO in studies using low mM concentrations of GABA could perhaps be explained by these variations. Given these discrepancies, our modeling of these effects is meant to be interpreted qualitatively. The details of the curve fits can be found in Supplementary Section 3.2.

While the C_dur parameter corresponding to Control synapses (C_durControl) for each GABA_A receptor-mediated synapse was set to 0.3 ms, the values of ḡ_max, α, and β corresponding to Control synapses (i.e., ḡ_maxControl, α_Control, β_Control) were taken from the values given in Tables 2–4. For post-ALLO synapses, each parameter was modified based on the relative changes between control GABA_A receptor activity and activity following the application of ALLO, as follows:

We specifically chose to model the effects of ALLO using synaptic and kinetic parameters fitted at physiological ALLO concentrations. A comparison of these parameters across physiological and supraphysiological ALLO levels is provided in Supplementary Figure S2.

To model the regional differences in neuronal heterogeneity observed experimentally, as discussed previously, we developed two distinct model configurations. Given that layer 5 IB neurons are more abundant in the frontal cortex and RS neurons in the parietal cortex, we modeled the frontal cortex using IB neurons while the parietal cortex was modeled using RS neurons. We chose this binary framework because it provides a manageable level of model complexity while still allowing us to probe specific mechanistic questions. Accordingly, we varied their proportions in each model such that the 95-5 (nIB:nRS) configuration is representative of a cortex most influenced by the frontal cortex, while the 5-95 (nIB:nRS) configuration reflects a cortical component that is most influenced by the parietal cortex. Given the differences in model architecture, the number of presynaptic neurons each postsynaptic neuron connects to (nPrePost) differs between the two model configurations, as presented in Supplementary Tables S3, S4.

Additionally, for each model configuration, we modeled enhanced connectivity in the frontal cortex by modifying the strength of the following synapses, as illustrated in Figure 2: IBIB, IBNRS, and NRSIB. Specifically, we modified the maximal synaptic conductances as follows: ḡ_IBIB was increased by a factor 7, ḡ_IBNRS by a factor of 5, and ḡ_NRSIB by a factor of 3. These specific adjustments were chosen for demonstration purposes, but similar results were observed across a physiologically reasonable range of synaptic strengths, as shown in Supplementary Figure S4. A detailed description of the method used to determine these parameter adjustments is provided in Supplementary Figure S4.

The goal of our model fitting was to replicate behavior that is qualitatively consistent with that of prior established models. Specifically, the model needed to exhibit two critical baseline states: a healthy state characterized by spindle oscillations with a network frequency of 7–10 Hz, and a diseased state marked by SWDs with a network frequency of 2.5–5 Hz. The 75-25 (nIB:nRS) model was simulated using GABA_A synapses that were all of the Control type, as described in Table 2. At 100% baseline cortical GABA_A conductance, the network produces spindle oscillations with a peak network frequency of 7 Hz, as shown in Figures 3A, C. In this healthy state, cell populations exhibit organized, alternating bursting patterns—particularly within the TC population—with moderately synchronous activity across other populations, as shown in Figure 3B. This activity is reflected in the LFP trace as a low-amplitude rhythmic oscillation. The network could be transformed to a diseased state by means of cortical disinhibition. By reducing GABA_A-receptor mediated inhibition in the cortex to 10% baseline conductance, the network produces a dominant 5 Hz SWD pattern as shown in Figures 3D, F. In this diseased state, neural firing becomes more synchronous across all populations, with distinctive spike-wave complexes visible in the LFP, as shown in Figure 3E. The reduction in cortical GABA_A conductance was determined by examining how changes in conductance influence network frequency and the relative power in the SWD frequency range, as shown in Figure 4.

Next, we modeled the effect of ALLO by altering parameters associated with the GABA_A current described in Section 2.4. The network was initialized to exhibit SWDs, as in Figures 3D–F, with cortical GABA_A conductance reduced to 10% baseline. The effect of ALLO was implemented uniformly across all GABA_A synapses in the model. Our results show that ALLO effectively suppresses pathological SWDs by altering the temporal firing patterns of neural populations throughout the thalamocortical circuit. As shown in Figure 5, there is a shift in the dominant network frequency to 7 Hz, which falls within the frequency band associated with normal spindle oscillations. This transition is also characterized by alternating bursts of synchrony across populations which is reflected in the low-amplitude filtered LFP trace. The shift from pathological to physiological dynamics is consistent with our previous findings and supports the therapeutic potential of positive GABA_A modulators in absence epilepsy.

Using the previously described frontal cortex composition (95-5 nIB:nRS) and parietal cortex composition (5-95 nIB:nRS) in Section 2.5, we next examined whether regional neuronal heterogeneity affects the network's propensity to generate the two critical network states of interest: healthy spindles and pathological SWDs. Additionally, we investigated whether the effect of ALLO differs across networks with different cortical compositions, and importantly, whether one network composition is more prone to non-resolution than the other.

We first investigated the excitability profiles of both network configurations by examining network frequency responses to varying levels of GABA_A-receptor mediated inhibition in the cortex, as shown in Figure 6. Both models display a transition from SWD to spindle activity as cortical inhibition increases. However, the thresholds at which this transition occurs differs slightly between the two, with the 95-5 (nIB:nRS) model exhibiting more intermediate oscillatory states compared to the 5-95 (nIB:nRS) model. The latter surprisingly maintains SWD activity across a broader inhibition range and requires stronger inhibition to shift to spindle-generating states. To further explore the influence of cortical composition, we simulated both networks under the same conditions as in Section 3.1—100% and 10% baseline cortical GABA_A conductance—to compare their abilities to sustain either physiological spindles or pathological SWDs, and to determine whether baseline dynamics differ across the two models. As shown in Figure 7, both network configurations exhibit qualitatively similar spectral characteristics, across the default mode and under cortical disinhibition. The spectral power distributions show similar frequency bands of activity (both exhibiting peak network frequencies of 7.5 Hz and 5 Hz under each condition), with only subtle differences in power distribution, likely due to variations in oscillatory patterns, as reflected in the LFP traces. This is also characterized by prominent power spectral density peaks within the 7–10 Hz (Figures 8A1, C1) and 2.5–5 Hz ranges (Figures 8A2, C2). At 100% baseline cortical GABA_A conductance, both models exhibit 7.5 Hz peak network frequencies with similar relative power in the spindles frequency range (Figure 8B1). Similarly, at 10% baseline conductance, both models exhibit network frequencies peaking at 5 Hz with similar relative power in the SWDs frequency range (Figure 8B2). These results suggest that despite substantial variation in cortical neuronal composition between models, the fundamental network dynamics remain conserved. This suggests that our fitted model (75-25 nIB:nRS configuration) is robust to changes in network architecture in producing baseline network states.

Following the approach in Section 2.4, we modeled the effect of ALLO in all GABA_A synapses, in each model configuration by first initializing to the diseased state. After implementing post-ALLO synapses in each compromised network, we observe a restorative effect, with each network transitioning back to spindle oscillations similar to the control state, as shown in Figures 8A3, C3. This transition is characterized by peak network frequencies of 7 Hz and 7.5 Hz for the 5-95 and 95-5 (nIB:nRS) models respectively, with comparable relative power in the spindles frequency range (Figure 8B3). Additionally, these qualitative restorative effects were preserved across the range of physiological and supraphysiological ALLO concentrations tested, with all levels showing ALLO-mediated resolution of SWDs, as shown in Supplementary Figure S3.

While the differences between the two model configurations is small, there is a more obvious difference between the control network states with cortical GABA_A at 100% and network states under the effect of ALLO. Specifically, the relative power in the spindles frequency range is higher in the control condition as compared to the post-ALLO state (Figures 8B1, B3). This is mainly due to the more sharply defined spectral density peak curves in the control condition, resulting in greater area under the curve within the spindles frequency range. In the post-ALLO state, the spectral profiles appear more spread out, resulting in relatively less concentrated power contribution falling within the spindles frequency range, despite restoring the dominant network frequency. Overall, these results suggest that, despite differences in the composition of cortical cell types (with one network having a higher proportion of IB cells), both networks respond similarly to cortical GABA_A modulation and the effect of ALLO, with neither showing significantly greater vulnerability to the non-resolution of SWDs.

In this part of the study, we focused on modeling the connectivity differences observed in treatment non-responders prior to treatment in newly diagnosed CAE patients, particularly the increased connectivity in the frontal cortex. By simulating a 50-50 (nIB:nRS) model configuration, we aimed to investigate how, in a network with a balanced composition of cortical cell types in layer 5, increased frontocortical connectivity might influence the modulation of SWDs—particularly in the non-resolution of SWDs (i.e., a lack of effect from ALLO on SWDs). Specifically, we modeled enhanced frontal cortical connectivity by modifying the following synaptic conductance parameters with a multiplicative factor to increase synaptic strength: (i) ḡ_IBIB·F_IBIB; (ii) ḡ_IBNRS·F_IBNRS; (iii) ḡ_NRSIB·F_NRSIB. The factors used and the rationale for their selection are detailed in Section 2.5. Additionally, we explored how these enhanced synapses affect the other two model configurations [5-95 and 95-5 (nIB:nRS)] to assess whether, in addition to connectivity differences, the composition of cell types within the network influences its responsiveness to the effects of ALLO.

Building on the approach in Section 3.1, we simulated the 50-50 (nIB:nRS) model with enhanced synaptic strengths, using (F_IBIB, F_IBNRS, F_NRSIB) = (7, 5, 3). We performed the simulation under the two conditions of cortical inhibition (100% and 10% baseline conductance) to investigate the impact of these changes on network activity. As shown in Figure 9, we observed distinct changes in spectral power across different frequencies, depending on the level of cortical inhibition and the presence of ALLO. Specifically, under full cortical inhibition, the network exhibits a maximal power at 5.5 Hz (Figure 9A), marking the first instance where the peak network frequency did not coincide with the frequency exhibiting the highest power. However, the relative power within the spindles range (0.282 V²) was greater than that within the range associated with SWDs (0.195 V²). When partial disinhibition was applied, the network shifted its maximal power to 4 Hz (Figure 9B). After introducing post-ALLO synapses, the network showed a lack of a restorative effect from ALLO, with the maximal frequency reducing to 2.5 Hz (Figure 9C). The continued SWD-like behavior is further reflected in the LFP traces of network activity post-ALLO application, especially when compared to network activity under cortical disinhibition, as shown in Figures 9B, C (in blue).

Next, we examined the impact of enhanced synaptic activity, with (F_IBIB, F_IBNRS, F_NRSIB) = (7, 5, 3), on the other two model configurations, 5-95 and 95-5 (nIB:nRS). As illustrated in Figures 10A1, B1, C1, both models showed similar network frequencies [6.5 Hz for 5-95 and 7 Hz for 95-5 (nIB:nRS)] and spindle-range relative power (0.339 V² and 0.376 V², respectively) under the condition of 100% baseline cortical GABA_A conductance. Reducing cortical GABA_A conductance to 10% resulted in comparable network behavior across models, with peak power at 4 Hz and 3 Hz, respectively, as shown in Figures 10A2, B2, C2. Both of these findings align with the results in Section 3.2, particularly the shift from spindle oscillations to SWDs under reduced cortical inhibition. Notably, the 95-5 (nIB:nRS) model consistently exhibited higher relative power across both frequency ranges, regardless of frontocortical connectivity alterations, as illustrated in Figures 11A, C. Without enhanced synaptic activity, network frequencies were identical between models and changed similarly across levels of cortical inhibition, as shown in Figures 11B, D (orange bars). However, when frontocortical connectivity was enhanced, the 95-5 (nIB:nRS) model demonstrated higher frequency under full inhibition and a lower frequency under partial disinhibition, as illustrated in Figures 11B, D (purple bars).

Introducing post-ALLO synapses in networks initialized to a pathological state, however, leads to interesting differences between the models. Surprisingly, the effect of ALLO on the 5-95 model mirrors the results from Section 3.2, where ALLO exerts a restorative effect, shifting the network's maximal power into the spindles range with a peak network frequency of 7.5 Hz, as shown in Figure 10A3. This shift occurs despite there being enhanced frontal cortical connectivity. In contrast, the effect of ALLO on the 95-5 (nIB:nRS) model resembles that seen in the 50-50 (nIB:nRS) model, with maximal power in the SWD-associated frequency range and a peak network frequency of 4.5 Hz.