To evaluate the contribution of GABA transporters to thalamic network activity in the context of GABA_B receptor-mediated inhibition, we used a standard, acute rat brain slice model in which electrical oscillations are evoked by extracellular stimulation of the synaptic inputs to the reticular thalamic nucleus in the presence of the GABA_A receptor blocker bicuculline (Huguenard and Prince, 1994; Jacobsen et al., 2001; Kleiman-Weiner et al., 2009). We evoked oscillations at intervals producing no rundown (once per minute, Jacobsen et al., 2001) and monitored neuronal activity with extracellular multiunit electrodes placed within the ventrobasal complex of the thalamus. By detecting evoked bursts, we found that oscillations last between 2 and 13 s at baseline (Figure 1A). The autocorrelogram of binned spike times revealed pronounced secondary peaks at multiples of approximately 500 ms (Figure 1A). After recording evoked oscillations for 20 min under baseline conditions, we then applied one of four experimental solutions to the perfusate (Figure 1B). Experimental solutions consisted of: (1) a control solution identical to the baseline solution, (2) 4 µM NO-711, a specific GAT1 blocker (Sitte et al., 2002), (3) 100 µM SNAP-5114, a specific GAT3 blocker (Borden et al., 1994), or (4) a combination of 4 µM NO-711 and 100 µM SNAP-5114. These blocker concentrations achieve full GAT blockade (Beenhakker and Huguenard, 2010). Experimental solutions were applied for 40 min, a duration sufficient for solution wash in and response stabilization (Beenhakker and Huguenard, 2010), and then washed out over another 20 min.

Figure 1

Download asset Open asset

Figure 1—source data 1

Oscillation measures in response to different GABA transporter blockade conditions.

Oscillation durations for each slice that was perfused with control, NO-711 (GAT1 blockade), SNAP-5114 (GAT3 blockade) or combined NO-711+ SNAP-5114 (dual blockade) (Figure 1D). Oscillation periods for each slice that was perfused with control, NO-711 or SNAP-5114 (Figure 1E). Oscillation durations for each slice that was perfused with combined NO-711+ SNAP-5114 for 60 min (Figure 1F). Phase 1 is baseline and phase 2 is drug perfusion.

https://cdn.elifesciences.org/articles/59548/elife-59548-fig1-data1-v2.zip

Download elife-59548-fig1-data1-v2.zip

When individually applied, either GAT1 or GAT3 blockade prolonged oscillations (Figure 1C), consistent with the absence seizure exacerbation seen with a clinically used GAT1 blocker, tiagabine (Ettinger et al., 1999; Knake et al., 1999; Vinton et al., 2005). We measured the duration of each evoked oscillation, then computed the average duration of the last five stable oscillations in baseline and experimental solutions (Figure 1D). Relative to baseline, individual GAT1 or GAT3 blockade increased oscillation duration by 36% (n = 8 slices from five animals, p=0.0027; here and in all results, percentages refer to relative change from control conditions) and 99%, (n = 7 slices from five animals, p=0.0076), respectively. We also evaluated the effects of individual GAT1 or GAT3 blockade on the period of evoked oscillations (Figure 1E). Relative to baseline, blocking either GAT1 or GAT3 increased the oscillation period by 13% (n = 8 slices from five animals, p=0.0021) and 32% (n = 7 slices from five animals, p=0.0050), respectively. Collectively, the effects of GAT blockade on oscillation properties generally agreed with the previously reported actions of GAT blockers on isolated GABA_B receptor-mediated IPSCs. That is, the 1.4-fold and 2-fold increase in oscillation duration corresponds roughly to the reported 1.5-fold and 2.2-fold increase in GABA_B IPSC amplitude produced by GAT1 or GAT3 blockade, respectively (Beenhakker and Huguenard, 2010). Additionally, GAT3 blockade significantly prolonged oscillation period, while the effect for GAT1 blockade was modest, consistent with reported effects on isolated GABA_B IPSC decay (Beenhakker and Huguenard, 2010).

Surprisingly, the effects of NO-711 and SNAP-5114 co-perfusion on evoked oscillations were not additive. Rather than prolonging evoked oscillations, dual GAT1+GAT3 blockade ultimately eliminated oscillations. Following a brief prolongation during the early phases of drug perfusion (see Figure 1C), dual GAT1+GAT3 blockade eventually decreased oscillation duration by 48% (Figure 1D, n = 9 slices from four animals, p=0.026). As the effects of dual blockade on oscillation duration did not appear to reach a steady state by 40 min, we extended the drug application to 60 min in a subset of experiments. For those slices, dual GAT1+GAT3 blockade invariably abolished oscillations (Figure 1F, n = 5 slices from three animals, p=0.0062).

In summary, the observed effects of individual GAT1 or GAT3 blockade on oscillation duration and period generally reflect the actions these individual blockers have on GABA_B receptor-mediated IPSCs isolated from thalamocortical neurons. GAT blockade-dependent increases in IPSC amplitude were associated with increased strength of oscillation, as measured by duration. However, the effects of dual GAT1+GAT3 blockade on oscillation duration did not reflect the additive effects of combined blockade on GABA_B IPSCs (Beenhakker and Huguenard, 2010). To better understand the discrepancy between GAT regulation of IPSCs and GAT regulation of thalamic oscillations, we next examined how IPSC amplitude and kinetics regulate the activity of single thalamocortical neurons.

We investigated the effects of GAT-modulated, GABA_B receptor-mediated currents on thalamocortical neuron rebound bursting, as this property is likely critical for the initiation of each successive cycle of an evoked oscillation (von Krosigk et al., 1993; Huguenard and Prince, 1994; Warren et al., 1994). Experimentally evoked GABA_B receptor-mediated IPSCs isolated in acute thalamic slices vary considerably in amplitude (Beenhakker and Huguenard, 2010), likely reflecting differences in synaptic activation of reticular thalamic neurons by the electrical stimulus. We therefore utilized an alternative approach to systematically examine the effects of GAT blockade on the firing properties of thalamocortical neurons: dynamic clamp (Sharp et al., 1993; Ulrich and Huguenard, 1996). We used GABA_B receptor-mediated IPSC waveforms isolated under voltage clamp during each pharmacological condition (control, GAT1 blockade, GAT3 blockade, dual GAT1+GAT3 blockade) as conductance waveform commands applied to single thalamocortical neurons (Figure 2A). We refer to these dynamic clamp-mediated conductance waveforms as dynIPSCs. We applied the dynIPSC corresponding to each pharmacological condition (dynControl, dynGAT1-Block, dynGAT3-Block, dynDual-Block; Figure 2B) to each recorded thalamocortical neuron.

Figure 2 with 1 supplement see all

Download asset Open asset

Figure 2—source data 1

LTS and burst features in response to IPSCs recorded with dynamic clamp (dynIPSCs).

Analyzed LTS and burst features for each recorded IPSC response. Averaged LTS and burst features for each neuron, for dynIPSCs scaled by 200% (Figure 2D, Figure 2—figure supplement 1B). Averaged LTS and burst features for each neuron, for dynIPSCs across all conductance amplitude scales (Figure 2E, Figure 2—figure supplement 1C).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig2-data1-v2.zip

Download elife-59548-fig2-data1-v2.zip

Since neurons likely receive variable numbers of inhibitory inputs, we scaled the conductance amplitudes for each pharmacological dynIPSC by 25%, 50%, 100%, 200%, 400% and 800%, yielding 24 possible dynIPSC waveforms (i.e. four pharmacological conditions x six amplitude scales). Additionally, we delivered each of the 24 waveforms at three approximate holding potentials: −60 mV, −65 mV or −70 mV. Five non-consecutive repetitions were performed for each dynIPSC waveform and holding potential condition. We found that differences across pharmacological dynIPSCs were most significant when the amplitudes were scaled by 200%. Throughout our analyses, we compare across dynIPSCs using this amplitude scale unless otherwise specified.

Figure 2C shows example responses to dynIPSCs delivered with dynamic clamp. Post-inhibitory, low-threshold calcium spikes (LTS) frequently followed each dynIPSC, with sodium-dependent action potentials often crowning each LTS. Herein, LTS refers to the slow, broad (~50 ms) event following inhibition. Burst, on the other hand, specifically refers to the collection of action potentials crowning the LTS. We quantified several properties of each post-inhibitory LTS and burst in response to each dynIPSC, including the probability of occurrence and the latency from dynIPSC onset (Figure 2C). We also computed LTS features such as peak voltage value, maximum rising slope and the number of spikes per LTS, averaged across trials for each neuron (Figure 2—figure supplement 1A).

We first examined LTS and burst probability distributions over all recorded neurons following delivery of dynIPSCs (LTS: Figure 2—figure supplement 1B–C, burst: Figure 2D–E). Relative to dynControl responses, LTS and burst probabilities were higher following either dynGAT1-Block (n = 47 cells, LTS: +26%, p=0.0080, burst: +63%, p=0.0018) or dynGAT3-Block (n = 47 cells, LTS: +39%, p=0.0015, burst: +106%, p=1.7×10⁻⁷), but lower following dynDual-Block (n = 47 cells, LTS: −88%, p=4.8×10⁻⁶, burst: −82%, p=0.030). We observed the same trend across pharmacological dynIPSCs for all other conductance scales (LTS: Figure 2—figure supplement 1C, burst: Figure 2E). Not surprisingly, increasing the conductance scale increased LTS and burst probability following either dynControl, dynGAT1-Block or dynGAT3-Block. However, both probabilities were consistently very low, below 6%, across all conductance scales following dynDual-Block. These changes in thalamocortical neuron rebound burst probability parallel the prolonged oscillation duration following individual GAT1 or GAT3 blockade, and the decrease in oscillation duration following dual GAT1+GAT3 blockade (Figure 1D).

Next, we examined distributions of average LTS and burst latencies over neurons responsive to dynIPSCs (LTS: Figure 2—figure supplement 1B–C, burst: Figure 2D–E). We restricted this analysis to dynGAT1- and dynGAT3-Block because dynDual-Block did not reliably evoke LTSs. Relative to dynControl responses, average LTS latency was not significantly different following dynGAT1-Block (n = 32 cells, p=0.97), while average burst latency was modestly prolonged (+4.6%, n = 21 cells, p=0.034). In contrast, dynGAT3-Block IPSCs significantly prolonged both LTS (+53%, n = 32 cells, p=3.7×10⁻⁹) and burst latency (+58%, n = 21 cells, p=1.1×10⁻⁹). We observed the same trend across pharmacological conditions for all other conductance amplitude scales (Figure 2E). As the inter-burst interval (latency from last burst) separates each cycle of seizure-like oscillations, and is dominated by inhibition of TC cells (Bal et al., 1995), the above results are consistent with the increased oscillation period following either individual GAT1 or GAT3 blockade, but a more pronounced effect for the latter (Figure 1E). Table 1 provides a summary of all LTS and burst features (see also Figure 2—figure supplement 1B–C).

In summary, the bidirectional differences in thalamocortical neuron rebound burst properties in response to different GABA_B activation waveforms (i.e. burst enhancement with either GAT1 or GAT3 blockade, but burst elimination during dual blockade) were correlated with the bidirectional differences in thalamic oscillation duration and period during the corresponding pharmacological manipulations. That is, by ultimately regulating thalamocortical neuron bursting, GATs appear to powerfully control thalamic network oscillations through differential activation of GABA_B receptors.

We next sought to determine the essential components of the thalamocortical neuron that contributes to the differential dynIPSC responses observed during our dynamic clamp experiments. We also sought to better understand the underlying channel dynamics contributing to the differential responses. Toward these ends, we established a conductance-based, multi-compartment, single neuron model for each of the 36 experimentally recorded thalamocortical neurons for which we had stable responses across all acquired conductance amplitude scales (Figure 3A).

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 3—source data 1

Parameters and errors for optimized model thalamocortical neurons.

Final error values between simulated and recorded data for each of 33 model thalamocortical neurons (Figure 3E). Optimized parameter values for each model thalamocortical neuron (Figure 3F).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig3-data1-v2.zip

Download elife-59548-fig3-data1-v2.zip

Our preliminary modeling results using existing TC cell models (Destexhe et al., 1998; Amarillo et al., 2014) failed to recapitulate two key features of GABA_B receptor-mediated post-inhibitory rebound LTSs that are likely critical in determining network level responses. Notably, the average LTS latencies of the model responses were routinely much earlier (400 ~ 1000 ms) than the biological ones (400 ~ 4000 ms). In addition, the model LTS and burst responses were graded in amplitude as a function of inhibitory strength, in contrast to the more characteristic all-or-none responses of recorded neurons. Therefore, we developed a gradient descent fitting approach to obtain suitable multicompartment models compatible with the data. As prior computational and experimental work demonstrates the importance of higher T channel densities in dendritic versus somatic compartments (Destexhe et al., 1998; Munsch et al., 1997; Williams and Stuart, 2000; Zhou et al., 1997), we modeled each thalamocortical neuron by a cylindrical somatic compartment and two cylindrical dendritic compartments in series (Figure 3C).

To reduce the number of fitted parameters, some simplifying assumptions were made: the somatic length and diameter were equivalent, the two dendritic compartments had equal dimensions, and passive leak channels were inserted into all three compartments at uniform densities. Four voltage-independent (passive) parameters were allowed to vary across model neurons: the somatic diameter (${diam}_{soma}$), the dendritic diameter (${diam}_{dend}$), the dendritic length (${length}_{dend}$) and the passive leak conductance density ($g_{pas}$). As prior work has identified ionic currents that contribute to the resting membrane potential of thalamocortical neurons (Amarillo et al., 2014), we inserted the following five voltage-dependent channels in all three compartments: the T-type calcium channel ($I_T$), the hyperpolarization-activated nonspecific cationic channel ($I_h$), the A-type transient potassium channel ($I_A$), the inward-rectifying potassium channel ($I_{Kir}$) and the persistent sodium channel ($I_{NaP}$). The densities of voltage-dependent channels were allowed to vary across compartments, resulting in a total of 15 (5 currents x 3 compartments) voltage-dependent (active) parameters that were allowed to vary across model neurons.

To provide an initial estimate of the geometric parameters that corresponded to each recorded thalamocortical neuron, we applied the short pulse methodology described by Johnston and Wu (1994, Chapter 4). During dynamic clamp experiments, we applied a short current pulse at the beginning of each recorded sweep. The average current pulse response for each neuron was well-fitted by a double exponential function (Figure 3B). From the coefficients and time constants of the two exponential components, we inferred the following four parameters for a ball-and-stick model (Rall, 1962): input conductance, electrotonic length, dendritic-to-somatic conductance ratio and the membrane time constant. These Rall model values were then converted into initial passive parameter seed values (${diam}_{Soma}$, ${diam}_{Dend}$, ${length}_{dend}$, $g_{pas}$) of each three-compartment model neuron (see Materials and methods). Notably, the average electrotonic length across all recorded neurons was 0.71 ± 0.08 space constants (range: 0.1-1.5 space constants), mitigating voltage attenuation concerns that might arise from applying currents through the soma with dynamic clamp.

Single thalamocortical neuron responses recorded during dynamic clamp experiments served to optimize passive and active parameters of each three-compartment, model thalamocortical neuron. GABA_B receptors were placed in the somatic compartment of each model neuron, and activation waveforms identical to the conductance waveforms used in dynamic clamp (i.e. dynIPSCs) were applied. We refer to these simulated GABA_B receptor activation waveforms as simIPSCs, corresponding to each pharmacological condition (simControl, simGAT1-Block, simGAT3-Block, simDual-Block). For each model neuron, responses to simIPSCs were iteratively compared to recorded dynIPSC responses. We evaluated the goodness-of-fit for each iteration with a total error defined by a weighted combination of component errors (see Materials and methods). Examples of resultant geometry and voltage response fits are shown in Figure 3C and D, respectively.

Each model neuron was trained using a set of 12 recorded traces, each selected from a different dynIPSC waveform, but evaluated against all 60–180 recorded traces for the neuron and ranked by the total error (Figure 3E). By removing neurons with a total error greater than two standard deviations above the mean, we designated the top 31 neurons as the set of well-fitted model neurons. Voltage response fits for selected well-fitted and excluded neurons are shown in Figure 3—figure supplement 1. All well-fitted neurons were characterized by high T channel densities in the distal dendrite and high A-type potassium channel densities in the soma (Figure 3F). Considerable variability among model neurons was observed in the densities of other ionic channels, likely contributing to the heterogeneity in LTS and burst measures among recorded neurons in response to each GABA_B dynIPSC waveform (Figure 2D).

We compared the distribution of LTS probabilities and features over the 31 well-fitted model neurons to the corresponding distributions derived from their biological counterparts recorded in dynamic clamp. In general, there was high agreement between the model simulations and dynamic clamp recordings (Figure 4A and B). Relative to sim/dynControl responses, LTS probability was increased following sim/dynGAT1-Block (n = 31, simulation: +30%, p=0.0086, dynamic clamp: +31%, p=0.038) or sim/dynGAT3-Block (simulation: +44%, p=7.4×10⁻⁴, dynamic clamp: +45%, p=0.011) and decreased following sim/dynDual-Block (simulation: −75%, p=9.8×10⁻⁷, dynamic clamp: −93%, p=5.7×10⁻⁴). Relative to sim/dynControl responses, average LTS latency was not different following sim/dynGAT1-Block (simulation: n = 27, p=0.85, dynamic clamp: n = 21, p=0.99) but was increased following sim/dynGAT3-Block (simulation: +42%, n = 27, p=7.3×10⁻⁸, dynamic clamp: +52%, n = 21, p=2.4×10⁻⁶). Differences in average number of spikes per LTS, average LTS peak value and average LTS maximum slope across dynIPSC waveforms were not sufficiently captured by model neurons (Figure 4—figure supplement 1). The same trends across pharmacological conditions were observed for all other conductance amplitude scales, showing high agreement between model and recorded neurons for LTS probability and latency, but not for other LTS features (Figure 4C and D, Figure 4—figure supplement 1).

Figure 4 with 1 supplement see all

Download asset Open asset

Figure 4—source data 1

Simulated and recorded LTS and burst features across well-fitted neurons.

Analyzed LTS and burst features for each simulated IPSC response with fast sodium-potassium mechanism (HH2.mod) inserted. Averaged LTS and burst features for each recorded and model neuron, for dynIPSCs scaled by 200% (Figure 4A–B, Figure 4—figure supplement 1A–B). Averaged LTS and burst features for each recorded and model neuron, for dynIPSCs across all conductance amplitude scales (Figure 4C–D, Figure 4—figure supplement 1C–D).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig4-data1-v2.zip

Download elife-59548-fig4-data1-v2.zip

In summary, many single neuron models were established that sufficiently recapitulated the probability and timing of post-inhibitory rebound bursts in response to 12 different physiological GABA_B-receptor IPSC waveforms. A commonality among well-fitted model neurons is that T channel densities were high in the dendrites and A-type potassium channel densities were high in the soma, while there was heterogeneity in other channel densities.

We next sought to understand the underlying mechanisms contributing to the different burst responses following different GAT-modulated, GABA_B-mediated IPSCs. The IPSC-evoked, post-inhibitory rebound LTS has been shown to be caused by T-type calcium currents (Llinás and Jahnsen, 1982; Kim et al., 2001, Figure 5—figure supplement 1A–C). While much is known regarding the gating properties of T channels, little is known how these properties behave in response to physiological inhibition. Although it has been proposed from artificial voltage ramp studies that rebound bursting is sensitive to the slope of voltage depolarization (Gutierrez et al., 2001), the underlying T channel dynamics that confer such voltage sensitivity remain unknown.

Low-threshold spikes are produced when T channel open probability discrepancy passes a threshold

To understand how IPSCs shape T channel dynamics, we first compared examples of LTS-producing responses evoked by simGAT1- and simGAT3-Block waveforms with examples of LTS-lacking responses evoked by simControl and simDual-Block waveforms in a well-fitted model neuron (Figure 5A). We tracked the activation ($m_T$) and inactivation ($h_T$) gating variables of the T channel, as a function of time (Figure 5B). By convention, $m_T=1$ when all channels are activated, and $h_T=0$ when all channels are inactivated (Hodgkin and Huxley, 1952). Also by convention, T channel open probability is given by ${m_T}^2{h}_T$ (Huguenard and McCormick, 1992). We also distinguished between steady-state values $\left[m_{T,\infty }\right(V)$ and $h_{T,\infty }\left(V\right)]$ that depend only on voltage, from instantaneous values $\left[m_T\right(V,t)$ and $h_T(V,t)$] that depend on both voltage and time. At every time point, instantaneous values attempt to reach steady-state values exponentially through voltage-dependent time constants $\left[{\tau }_m\right(V)$ and ${\tau }_h\left(V\right)$].

Figure 5 with 1 supplement see all

Download asset Open asset

Figure 5—source data 1

Open probability discrepancy measures for all simulated traces.

Analyzed open probability discrepancy measures and LTS features for each simulated IPSC response with no fast sodium-potassium mechanism inserted. These are averaged for each model thalamocortical neuron, for traces with and without an LTS (Figure 5C, Figure 5—figure supplement 1).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig5-data1-v2.zip

Download elife-59548-fig5-data1-v2.zip

One notable feature of the T channel is that the activation time constant ${\tau }_m$ is rapid, indicating that the instantaneous activation variable $m_T$ quickly achieves its voltage-dependent steady-state value $m_{T,\infty }$. Thus, throughout the course of the relatively slow membrane voltage changes associated with the GABA_B-receptor-mediated IPSC, $m_T$ was nearly identical to $m_{T,\infty }$ [Figure 5B(i–ii)]. By contrast, the T channel inactivation time constant ${\tau }_h$ is slow, with a value about 10-fold higher than ${\tau }_m$ (Coulter et al., 1989). Thus, during the IPSC, the instantaneous T channel inactivation variable $h_T$ rarely achieved its voltage-dependent steady-state value $h_{T,\infty }$ [Figure 5B(iii) and (iv)].

At each time point, the capacity of $h_T$ to achieve $h_{T,\infty }$ can be quantified by calculating the difference between the two values. If the discrepancy between the two values is low (i.e. $h_T-h_{T,\infty }\approx 0$), then $h_T$ has successfully achieved its steady-state value at that time point. A high discrepancy, on the other hand, reflects an $h_T$ that has incompletely reached steady state. Importantly, the magnitude of discrepancy between $h_T$ and $h_{T,\infty }$ affected the T channel open probability ${m_T}^2{h}_T$, but only when the activation variable $m_T$ was high. We therefore quantified the difference between instantaneous versus steady-state open probability (${m_T}^2{h}_T-{m_{T,\infty }}^2{h}_{T,\infty }$), a parameter we simply refer to as T channel open probability discrepancy [Figure 5B(vi)]. When all 31 well-fitted model neurons were considered, the maximum T channel open probability discrepancy ${\displaystyle \left[\underset{t}{max}\left(m_T^2{h}_T-m_{T,\mathrm{\infty }}^2{h}_{T,\mathrm{\infty }}\right)\right]}$ was on average 2.9 orders of magnitude higher for LTS-producing responses than for LTS-lacking responses (Figure 5C, n = 31 cells, p = 2.4 x 10⁻²⁹). In fact, when all traces were considered, a threshold open probability discrepancy of 10⁻² separated LTS-producing responses from LTS-lacking responses (not shown, herein referred to as the observed LTS threshold).

Two criteria drive open probability discrepancy past LTS threshold

As high T channel open probability discrepancy predicted an LTS response, we sought to determine criteria that allowed this measure to pass LTS threshold (Figure 5D). For all simIPSC conditions, when voltage depolarized during the decay phase of the IPSC, open probability discrepancy increased as ${\displaystyle h_T}$ (which was high from the hyperpolarization) slowly approached low $h_{T,\infty }$ (Figure 5D(ii),(iii): label a); discrepancy was largely defined by the slow kinetics of $h_T$ because $m_T$ matches $m_{T,\infty }$ at all time points. However, near the end of the IPSC, the open probability discrepancy curve for LTS-producing responses reached an inflection point (point of zero concavity, Figure 5D(ii),(iv): circles) that progressed towards positive concavity (concave upward, Figure 5D(ii),(iv): label b) to eventually pass LTS threshold [blue and red curves]. In contrast, open probability discrepancy curves in LTS-lacking responses either failed to reach an inflection point [Figure 5D(ii),(iv): black curve, label c], or reached an inflection point but veered back towards negative concavity [concave downward, Figure 5D(ii),(iv): purple curve, label d]. In either case, the response failed to reach LTS threshold. We define the time point at which the open probability discrepancy curve reaches zero (or maximal negative) concavity as the decision point (circles in Figure 5D).

A comparison between the simGAT3-Block (i.e. LTS-producing, red) and simDual-Block (i.e. LTS-lacking, purple) responses showed that at the decision point, the slope of open probability discrepancy, rather than its value, predicted whether positive or negative concavity was ultimately achieved [Figure 5D(iv), circles]. A rapidly changing open probability discrepancy (i.e., steep slope) causes a greater acceleration (i.e. positive concavity) in this measure. In fact, for all open probability curves reaching zero concavity, the slope of the open probability discrepancy curve at the decision point was always higher for LTS-producing responses than for LTS-lacking responses. We show the progression of trajectories aligned to the decision points for two conditions (Video 1, simGAT3-Block and simDual-Block) and for all conditions (Video 2). When all 31 well-fitted model neurons were considered, there was a significant difference between LTS-producing and LTS-lacking responses for either the maximum open probability discrepancy concavity (Figure 5—figure supplement 1D), the discrepancy slope at the decision point (Figure 5—figure supplement 1E) or the voltage slope at the decision point (Figure 5—figure supplement 1F). We conclude that an LTS is produced only if two criteria are satisfied: (1) the open probability discrepancy curve reaches an inflection point (e.g. Figure 5D(iv): blue, red and purple curves), and (2) the slope of the open probability discrepancy curve at that inflection point is high (e.g. Figure 5D(iv): blue and red curves).

Video 1

Download asset

Video 2

Download asset

Open probability discrepancy depends on T channel inactivation kinetics

To test the contribution of T channel inactivation kinetics to LTS production, we bidirectionally altered the T channel inactivation time constant (${\tau }_h$) for the same simIPSC response simulations as in Figure 5. When ${\tau }_h$ was halved in simGAT1- and simGAT3-Block simulations, T channel open probability discrepancy remained less than LTS threshold and LTS responses normally observed during GAT1 and GAT3 blockade were abolished (Figure 6A with Figure 5A–B). In contrast, doubling ${\tau }_h$ in simDual-Block simulations drove T channel open probability discrepancies past LTS threshold and, consequently, normally absent LTSs appeared (Figure 6B with Figure 5A–B). Note that simControl waveforms continued to result in LTS failure as voltage hyperpolarization was weak and T channel recovery was low.

Figure 6

Download asset Open asset

Open probability discrepancy depends on IPSC kinetics

To test the contribution of inhibition kinetics to T channel open probability discrepancy, we systematically varied the time constant of IPSCs while fixing the amplitude. Prolonging an LTS-producing simGAT3-Block IPSC decreased T channel open probability discrepancy and abolished LTS responses as time constants increased above four-fold (Figure 7A1). Conversely, shortening a non-LTS-producing simDual-Block IPSC increased T channel open probability discrepancy and produced LTS responses as time constants decreased by 20% (Figure 7A2). We show the progression of trajectories aligned to the decision points for two time constants (Video 3, last LTS success and first LTS failure) and for all time constants (Video 4). As changing the kinetics of inhibition also changes the total amount of inhibition delivered to a cell, we also changed inhibition kinetics while keeping charge (i.e. the area under the curve) constant. Nonetheless, we continued to observe a decrease in T channel open probability discrepancy and eventual LTS failure as simIPSC kinetics slowed (Figure 7—figure supplement 1). Therefore, inhibition kinetics appear to be important for LTS production through its influence on T channel open probability discrepancy.

Figure 7 with 1 supplement see all

Download asset Open asset

Video 3

Download asset

Video 4

Download asset

Since hyperpolarization promotes T channel recovery (Coulter et al., 1989), it remains possible that sufficiently strong hyperpolarization – regardless of waveform – will produce an LTS. To test this possibility, we varied the inhibition amplitude using either the simGAT3-Block (fast kinetics) or simDual-Block (slow kinetics) waveform. As we increased the amplitude of the simGAT3-Block waveform while fixing the rise and decay time constants, LTSs emerged as T channel open probability discrepancy increased (Figure 7B1). In contrast, as the amplitude of the simDual-Block waveform increased while fixing the rise and decay time constants, T channel open probability discrepancy nonetheless remained low and robust LTS responses never emerged (Figure 7B2). Thus, fast inhibition kinetics are important for driving T channel open probability discrepancy, and slow kinetics provides an explanation for the consistently low LTS and burst probability across all conductance amplitude scales following dyn/simDual Block (Figures 2E, 4C and D).

In summary, LTS production following physiological inhibition is largely controlled by the dynamics of T channel open probability discrepancy. We show that LTS production appears to depend not only on the amplitude of inhibition, but also the temporal envelope of inhibition. Large inhibition amplitude is required for sufficient T channel recovery, whereas fast inhibition decay is required for driving the T channel open probability discrepancy beyond an inflection point and creating a brief time window with sufficiently high T channel open probability for LTS production.

We next explored whether the interplay between GABA_B-mediated inhibition and T type calcium channel dynamics in thalamocortical neurons contributes to the observed changes in network-level oscillations following GAT blockade. We first examined the effects of GABA_B receptor-mediated inhibition in a simplified two-cell network configuration. In each two-cell network, we connected a single compartment, GABAergic reticular thalamic model neuron (Klein et al., 2018) to one of the 31 well-fitted model thalamocortical neurons (Figure 8A). To generate action potentials, Hodgkin-Huxley type sodium and potassium channels were inserted into the somatic compartment of each model neuron (Williams and Stuart, 2000). The reticular thalamic neuron was connected to the thalamocortical neuron only via a GABA_B receptor-mediated inhibitory synapse (GABA_A receptors were blocked during experimentally evoked oscillations, see Figure 1). Consistent with previous intra- thalamic models (Destexhe et al., 1996), the thalamocortical neuron provided AMPA receptor-mediated excitation to reticular thalamic neurons. By applying a brief stimulating current to the reticular thalamic neuron and varying the GABA_B receptor activation parameters (simIPSCs), GABA_B conductance waveforms comparable to those used by dynamic clamp were evoked in each thalamocortical neuron (Figure 8B). To simulate variable resting membrane potentials of thalamocortical neurons, we varied the thalamocortical neuron leak reversal potential between −73 and −60 mV (14 leak reversal potentials). To generate trial-to-trial variability (five trials per leak reversal potential), we randomized the leak conductance of each neuron to within 10% of the original value. Of all 31 possible two-cell networks, 24 had a quiescent, pre-stimulation baseline over this range of leak reversal potentials. In those networks with quiescent baseline, oscillations persisted only when the GABA_B-mediated inhibition promoted T channel open probability discrepancy in thalamocortical neurons to produce rebound bursting (Figure 8C), consistent with our single neuron models. An oscillation probability was computed for each two-cell network over the 14 × 5 = 70 trials. In addition, an oscillatory period and an oscillatory index based on the autocorrelation function of pooled spikes was computed for each successfully evoked oscillation (see Materials and methods).

Figure 8 with 1 supplement see all

Download asset Open asset

Figure 8—source data 1

Oscillation measures for two-cell model thalamic networks using different GABA_B receptor activation waveforms (simIPSCs).

Oscillation measures for each network simulation. Averaged oscillation measures for each two-cell network (Figure 8D).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig8-data1-v2.zip

Download elife-59548-fig8-data1-v2.zip

We examined the distributions of oscillation probability, average oscillation period and average oscillatory index over the 24 different two-cell networks (Figure 8D). Relative to using simControl parameters, oscillation probability increased when using simGAT1-Block parameters (+66%, n = 24 networks, p=0.016) or simGAT3-Block parameters (+93%, p=0.048). These results are consistent with the experimental observation that individual GAT1 or GAT3 blockade prolonged oscillations (Figure 1D). Relative to using simControl parameters, average oscillation period increased when using simGAT3-Block parameters (+39%, n = 10, p=4.2×10⁻⁴), also consistent with the experimental observation that individual GAT3 blockade increased oscillation periods (Figure 1E). In contrast, when simDual-Block parameters were applied in the network, oscillations did not arise, consistent with the experimental observation that dual GAT1+GAT3 blockade inevitably abolished oscillations (Figure 1D and F). Additionally, across all two-cell networks, we found that oscillations only appeared when T channel open probability discrepancy passed the observed LTS threshold of 10⁻² identified earlier (Figure 8—figure supplement 1B, Figure 5C). Indeed, T channel open probability discrepancy correlated well with oscillation probability (Figure 8E) and underscores the importance of the interplay between GABA_B-mediated inhibition and T channel gating in regulating oscillations. Although the two-cell networks recapitulated several effects of GAT blockade on oscillations, the oscillations generated by each two-cell network were nonetheless extremely stereotyped and regular, resulting in unrealistically high oscillatory indices (Figure 8—figure supplement 1A). Moreover, network responses across all two-cell networks were highly heterogeneous, with many of them failing to oscillate at the conductance amplitude scale used.

We sought to determine whether larger, more complex model networks could more realistically simulate experimental oscillations and recapitulate the bidirectional effects of GAT blockade by varying simIPSCs. We scaled up the network to include one circular layer of 100 reticular thalamic (RT) neurons and one circular layer of 100 thalamocortical (TC) neurons (Figure 9A). RT-TC inhibitory connections and TC-RT excitatory connections were both convergent and divergent. To assess the importance of the geometric and conductance heterogeneity we observed in the single cell models (Figure 3F), we established two sets of model thalamic networks: (1) 24 TC-homogeneous networks with TC parameters taken one at a time from each of the 24 model neurons used in the two-cell networks and (2) 24 TC-heterogeneous networks with TC parameters taken from all of the 24 model neurons, randomly ordered. All networks had a quiescent, pre-stimulation baseline when the thalamocortical neuron leak reversal potential was varied between −73 and −60 mV (14 leak reversal potentials). To generate trial-to-trial variability (five trials per leak reversal potential), we randomized the leak conductance for each of the 200 neurons to within 10% of the original value. For both TC-homogeneous and TC-heterogeneous networks, oscillations emerged and spread in response to some but not all GABA_B receptor activation parameters (Figure 9B and C). An oscillation probability was computed for each 200 cell network over the 14 × 5 = 70 trials. In addition, an oscillatory period, an oscillatory index and a half activation time was computed for each successfully-evoked oscillation (see Materials and methods).

Figure 9

Download asset Open asset

Figure 9—source data 1

Oscillation measures for 200 cell model thalamic networks using different GABA_B receptor activation waveforms (simIPSCs).

List of candidate model thalamocortical neurons. Oscillation measures for each network simulation, for both TC-homogeneous and TC-heterogeneous networks. Averaged oscillation measures for each TC-homogeneous 200 cell network (Figure 9D). Averaged oscillation measures for each TC-heterogeneous 200 cell network (Figure 9E).

https://cdn.elifesciences.org/articles/59548/elife-59548-fig9-data1-v2.zip

Download elife-59548-fig9-data1-v2.zip

We examined the distributions of oscillation measures over the set of TC-homogeneous networks (Figure 9D) and the set of TC-heterogeneous networks (Figure 9E). The values of oscillation periods and oscillatory indices for TC-heterogenous networks were similar to values extracted from experimental recordings (oscillation period: Figure 1E, oscillatory index: not shown). In response to the same simIPSC conditions, we observed highly varied (often bi-modal) oscillation responses across TC-homogeneous networks, which reflects the highly varied LTS responses across individual model TC neurons (Figure 4A). In contrast, oscillation measures were less variable across the different TC-heterogeneous networks. Therefore, cell heterogeneity averages out the LTS response variability, provides more robust network responses and amplifies the differences across simIPSC conditions. Indeed, for the set of TC-heterogeneous networks, relative to using simControl parameters, oscillation probability increased when using either simGAT1-Block (+8.2%, n = 24 networks, p=0.049) or simGAT3-Block parameters (+8.5%, p=0.013), but decreased when using simDual-Block parameters (−82%, p=0.0010). Consistent with this result, the percent of active TC cells increased when using either simGAT1-Block (+135%, p=0.0044) or simGAT3-Block parameters (+193%, p=1.6×10⁻⁵), but decreased when using simDual-Block parameters (−79%, p=0.037). Finally, average oscillation period increased when using either simGAT3-Block (+23%, p=1.3×10⁻⁷) or simDual-Block parameters (+50%, p=3.8×10⁻⁹). These results are consistent with experimental findings that individual GAT1 or GAT3 blockade increased oscillation durations (Figure 1D) and oscillation periods (Figure 1E), whereas dual GAT1+GAT3 blockade eliminated oscillations (Figure 1F) in acute thalamic slices.

In summary, a population of thalamic network models was established that sufficiently recapitulated the bidirectional effects of individual versus dual GAT blockade on thalamic oscillations by merely altering the kinetics of GABA_B-receptor inhibition. The same interplay between GABA_B-receptor-mediated inhibition and T channel open probability that governs thalamocortical neuron rebound bursting appears to regulate network-level oscillations. Furthermore, we found that including cell heterogeneity in network simulations provides more robust oscillations and more realistically recapitulates experimental oscillation periods by averaging out LTS response heterogeneity. Using heterogeneous networks in future studies would thus facilitate comparison across experimental conditions and provide more meaningful predictions.

Thalamocortical neuron rebound bursting has long been implicated in spike-wave discharges (SWDs) observed in generalized seizures, and T channels mediate thalamocortical neuron rebound bursting (Kim et al., 2001; Porcello et al., 2003). The expression of the CaV3.1 T channel subtype by thalamocortical neurons correlates with SWD expression in both animal models (Kim et al., 2001; Broicher et al., 2008; Ernst et al., 2009) and human patients (Singh et al., 2007). Blocking T channels reduces both thalamocortical neuron bursting and oscillations in thalamic slice models (Huguenard and Prince, 1994). However, the importance of thalamocortical neuron bursting in generalized seizures remains unresolved. One recent study found reduced thalamocortical neuron firing during SWDs and a lack of seizure reduction with weak, local T channel blockade in the ventrobasal nucleus (McCafferty et al., 2018; however, strong blockade diminished SWDs), while a second recent study found that increasing/decreasing thalamocortical neuron bursting increases/decreases SWDs in both epileptic mice and rats (Sorokin et al., 2017). One possibility accounting for these discrepant findings is that the population of active thalamocortical neurons during SWDs is sparse (Huguenard, 2019). Interestingly, our TC-heterogenous network models – models that are likely more representative of biological networks – were largely characterized by robust, yet sparse, oscillations (in contrast, TC-homogenous network models produced oscillations that were widespread, Figure 9B with Figure 9C).

Our study found a strong relationship between thalamocortical neuron bursting and epileptiform thalamic oscillations. First, the pharmacological conditions in which oscillation durations were increased (individual GAT1 or GAT3 blockade) or decreased (dual GAT1+GAT3 blockade) relative to baseline were the same conditions in which thalamocortical neuron rebound burst probability increased (individual GAT1- or GAT3-Block GABA_B dynIPSC waveform) or decreased (Dual-Block GABA_B dynIPSC waveform), relative to control dynIPSCs. Second, the pharmacological conditions in which oscillation periods were increased (individual GAT1 or GAT3 blockade) relative to baseline were also the same conditions in which thalamocortical neuron rebound burst latencies increased (individual GAT1- or GAT3-Block dynIPSC waveform) relative to control dynIPSCs. Finally, in our two-cell model thalamic network, each successive oscillation cycle was initiated by a thalamocortical neuron rebound burst (Figure 8C), similar to what has been reported in experiments when a reticular thalamic neuron and a thalamocortical neuron were simultaneously recorded during a thalamic oscillation (Bal et al., 1995).

Prior experimental and computational work has suggested that a shift from GABA_A receptor-mediated to GABA_B receptor-mediated inhibition at the RT-TC synapse transforms oscillations in acute thalamic slices from a 10 Hz, sparse, spindle-like activity to a 3 Hz, hyper-synchronized, seizure-like state (von Krosigk et al., 1993; Destexhe et al., 1996; Destexhe, 1998; Blumenfeld and McCormick, 2000). Notably, the shift in oscillation frequency is consistent with differences in IPSC decay constants (GABA_A: <100 ms; GABA_B: about 300 ms) recorded in thalamocortical neurons (Huguenard and Prince, 1994). Multiple animal model studies support the hypothesis that generalized spike-wave seizures rely on robust GABA_B receptor-mediated inhibition. Specifically, systemic injection of GABA_B receptor agonists increases SWDs (Liu et al., 1992; Bortolato et al., 2010), while injection of GABA_B receptor antagonists reduces or even abolishes SWDs (Liu et al., 1992; Vergnes et al., 1997). Notably, however, other studies were not able to record rhythmic GABA_B IPSCs during SWDs in vivo (Charpier et al., 1999) or find significant SWD changes with GABA_B receptor modulation (Staak and Pape, 2001), and have highlighted a particular role for extrasynaptic GABA_A receptors that produce a non-dynamic, tonic current in thalamocortical neurons (Cope et al., 2009). Presumably, tonic GABA_A currents hyperpolarize the resting membrane potential, promoting T channel recovery and increasing rebound bursting (Cope et al., 2005). Thus, the disparate conclusions regarding the importance of GABA_A receptor- versus GABA_B-receptor-mediated inhibition appear to nonetheless converge on similar conclusions regarding the importance of T channel recovery, a process that involves membrane potential hyperpolarization. However, we show here that hyperpolarization must be dynamic for T channels to produce an LTS. Under physiological conditions, it seems reasonable to expect multiple, convergent inhibitory mechanisms that promote T channel recovery.

Our study establishes a clear link between GABA_B receptor-mediated inhibition, thalamocortical neuron rebound bursting and epileptiform thalamic oscillations. First, we found a 1.4-fold or a twofold increase in oscillation duration in acute thalamic slices after perfusing with the GAT1 blocker NO-711 or the GAT3 blocker SNAP-5114 (Figure 1D), closely corresponding with the 1.5-fold or 2.2-fold increase in GABA_B IPSC amplitude recorded under the same conditions, respectively (Beenhakker and Huguenard, 2010). Second, in both dynamic clamp recordings and model neuron simulations, the GAT1-Block and GAT3-Block GABA_B IPSC waveforms produced higher thalamocortical neuron rebound LTS or burst probability relative to the Control waveform, with the GAT3-Block waveform increasing LTS or burst probability more (Figure 2D and Figure 4A–B). Third, we found that increasing conductances for each of the Control, GAT1-Block and GAT3-Block waveforms led to an increase in LTS or burst probability, for both dynamic clamp recordings and model neuron simulations (Figure 2E and Figure 4C–D). Finally, a comparison of the Control (black) versus GAT1-Block (blue) IPSC responses in Figure 5 shows that the lack of LTS in the former correlates with a decreased level of TC hyperpolarization, a lack of T channel de-inactivation (both $h_T$ and $h_{\infty ,T}$ are low) and a deficiency in T current production (no spike in $I_T$), agreeing with prior studies that show how the initiation of a low-threshold rebound spike depends on the sufficient removal of T channel inactivation through membrane potential hyperpolarization (Llinás and Jahnsen, 1982; Coulter et al., 1989).

More interestingly, we found that not only is the overall amount of inhibition important, how such inhibition distributes over time is equally important. For instance, the dynDual-Block waveform has an area under the curve (i.e. ≈ charge) about twice that of the dynGAT3-Block waveform (Figure 2B). Nevertheless, both dynamic clamp recordings and model neuron simulations showed that the dynDual-Block waveform decreased rebound burst probability relative to the control, whereas the dynGAT3-Block waveform increased rebound burst probability (Figure 2D and Figure 4). In fact, the dynDual-Block waveform largely failed to produce rebound bursts even when the conductance amplitude was scaled so high that the burst probability was close to one in all other conditions, that is following dynControl, dynGAT1-Block and dynGAT3-Block waveforms scaled at 400%+ (Figure 2E and Figure 4). Dual GAT blockade also abolished oscillations, in stark contrast to the robust prolongation of oscillations observed during GAT3 blockade only (Figure 1D and F). Therefore, even high levels of synaptic inhibition, if decayed too slowly, can abolish both thalamocortical neuron rebound bursts and epileptiform thalamic oscillations, a conclusion recapitulated by our model thalamic networks (Figure 8D and 9E).

Our model neurons allowed us to identify a novel measure of T channel gating that provides an explanation for how GAT-modulated, GABA_B IPSCs bidirectionally regulates LTS production. We found that when the discrepancy between high instantaneous T channel open probabilities (${m_T}^2{h}_T$) and low steady-state open probabilities (${m_{T,\infty }}^2{h}_{T,\infty }$) passed a threshold of 10⁻², thalamocortical neurons produced an LTS (Figure 5C) and two-cell thalamic networks produced oscillations (Figure 8—figure supplement 1B). To reach LTS threshold, we identified two criteria that must occur in the open probability discrepancy curve: (1) reach an inflection point and (2) the inflection point must have a steep slope. We compared three qualitatively different IPSC response trajectories in Figure 5 and Video 2. We observed that simControl (black) did not produce an LTS response simply because inhibition was insufficient for T channel recovery ($h_T$ was always below 0.2, Figure 5B). Under these circumstances, the open probability discrepancy curve never reached zero concavity [failed criterion (1)]. Although T channels were sufficiently recovered ($h_T$ reached above 0.6) by strong hyperpolarization associated with simGAT1-, simGAT3-, and simDual-Block, only the former two waveforms produced an LTS response. Rapid depolarization during the decay of simGAT1- and simGAT3-Block (blue and red) caused the open probability discrepancy curve to not only reach an inflection point [satisfied criterion (1)], but also had a steeper slope and progressed towards positive concavity [satisfied criterion (2)], ultimately driving the open probability past LTS threshold. In contrast, the simDual-Block waveform (purple) produced a prolonged inhibition. As a result of the slower depolarization during IPSC decay, the open probability discrepancy curve did reach an inflection point [satisfied criterion (1)], but had a shallower slope at that inflection point [failed criterion (2)], causing the T channel open probability discrepancy to veer away from LTS threshold.

We conclude that a combination of high inhibition amplitude and fast inhibition kinetics is necessary for LTS production. A high inhibition amplitude causes sufficient T channel recovery, which appears to be necessary for the open probability discrepancy curve to reach an inflection point [criterion (1) above]. Fast inhibition kinetics causes rapid depolarization during IPSC decay, which appears to be necessary for the slope of the open probability discrepancy curve to be high at the inflection point [criterion (2) above], eventually driving the open probability discrepancy past LTS threshold. This open probability discrepancy depends on a slow T channel inactivation time constant, as doubling the T channel inactivation time constant produces LTSs in response to the simDual-Block waveform (Figure 6B), a waveform that normally results in LTS failure. In order to confirm that inhibition amplitude and kinetics influence LTS production by controlling T channel opening, we varied the GABA_B simIPSC waveform systematically five different ways (Figure 7A–B, Figure 7—figure supplement 1). In each case, an LTS was only produced when there was both sufficient T channel recovery and an increase in T channel open probability discrepancy. Collectively, the results of these manipulations (Figure 7B, Video 3 and Video 4) support the hypothesis that dual GAT1+GAT3 blockade eliminates thalamic oscillations by promoting sustained, GABA_B receptor-mediated inhibition that, in turn, promotes a convergence of instantaneous and steady-state T channel open probabilities, a condition that results in LTS failure.

Prior experimental work has shown that T channels are critical for the generation of post-inhibitory low-threshold rebound spikes in thalamocortical neurons (Kim et al., 2001; Porcello et al., 2003). A computational study by Destexhe et al., 1998 showed that membrane voltage trajectories of low-threshold rebound spikes in thalamocortical neurons can be recapitulated using a three-compartment model, but not a single-compartment model. A recent study by Amarillo et al., 2014 used pharmacological approaches to identify seven intrinsic channels ($I_T$, $I_h$, $I_A$, $I_{Kir}$, $I_{NaP}$, $I_{NaLeak}$, $I_{KLeak}$) contributing to the resting membrane potential of thalamocortical neurons. The same study created a single-compartment computational model neuron that produced thalamocortical neuron rebound bursting and showed that a balance between the five voltage-dependent channels shapes the low-threshold rebound spike.

Our model thalamocortical neurons extend the studies by Destexhe et al., 1998 and Amarillo et al., 2014. That is, we incorporated the same intrinsic channels described by Amarillo et al., 2014, but allowed channel density to vary across three compartments [note that the two different leak channels in Amarillo et al. were simplified by combining them into a single passive leak conductance ($g_{pas}$) with a reversal potential ($E_{pas}$)]. Furthermore, we established a different model neuron for each set of single thalamocortical neuron recordings (Figure 3), in an effort to capture the heterogeneity in response to physiologically relevant inhibitory inputs such as GABA_B receptor-mediated IPSCs (Figure 2). By creating a total of 31 well-fitted model neurons, we were also able to generate population statistics that were in many ways similar to those of the corresponding population of recorded neurons (Figure 4). Having many distinct model thalamocortical neurons enabled us to introduce TC neuron heterogeneity into thalamic networks, which we found essential for robustly reproducing oscillations and capturing the differences across GAT blockade conditions seen in our slice recordings (Figure 9). Future work may utilize the distinct sets of optimized parameters (Figure 3F) to probe for mechanisms underlying the high IPSC response heterogeneity seen in thalamocortical neurons in response to the same IPSCs (Figure 2D).

Drugs that increase synaptic inhibition in the brain may paradoxically exacerbate generalized seizures. For example, the GAT1 blocker tiagabine increases extracellular GABA concentrations (Beenhakker and Huguenard, 2010) and is an effective anti-epileptic drug used to treat focal epilepsy (Nielsen et al., 1991; Uthman et al., 1998). However, tiagabine can also induce continuous SWDs (i.e. absence status epilepticus) in patients with absence epilepsy (Knake et al., 1999) and non-convulsive status epilepticus in patients with focal epilepsy (Ettinger et al., 1999; Vinton et al., 2005). GAT1 knockout mice, which have increased inhibitory currents, and presumed increased activation of GABA_B receptors (Beenhakker and Huguenard, 2010) in thalamocortical neurons, also exhibit increased incidence of SWDs (Cope et al., 2009). Furthermore, application of GABA_B receptor antagonists increases seizures in rats susceptible to convulsive focal seizures, but suppresses seizures in rats with non-convulsive absence seizures (Vergnes et al., 1997). Together, these observations suggest that, although a seizure is a manifestation of hyperexcitable neuronal firing, increasing inhibition in the form of neuronal hyperpolarization may not always reduce seizures (Beenhakker and Huguenard, 2009).

In fact, as we discovered in this study, inhibition kinetics seems to play an important role in determining whether thalamocortical rebound bursts and seizure-like oscillations occur. The finding that dual GAT1+GAT3 blockade abolished epileptiform oscillations may at first glance suggest that nonspecific GABA transporter blockers may be used to treat generalized seizures. However, this approach would be relatively non-specific and, for example, likely disrupt normal spindle oscillation formation, leading to undesirable side effects such as sedation and coma. Nevertheless, it is possible that temporary alteration of GABA_B inhibition kinetics, either through GABA transporter modulation or other means, may create pockets of IPSC response heterogeneity in the overall thalamocortical neuron population, making it less likely for the thalamic network to develop a hypersynchronous state (Pita-Almenar et al., 2014).

Sprague-Dawley rats [Charles River Laboratories (RRID:RGD_734476), Wilmington, MA] of postnatal day 11 to 17 (P11-P17), of either sex, were used in oscillation experiments, which were performed in accordance with protocols approved by the Institutional Animal Care and Use Committee at the University of Virginia (Charlottesville, VA). Rats were deeply anesthetized with pentobarbital, then transcardially perfused with ice-cold protective recovery solution containing the following (in mM): 92 NMDG, 25 glucose, 5 Na-ascorbate, 3 Na-pyruvate, 2.5 KCl, two thiourea, 20 HEPES, 26 NaHCO₃, 1.25 NaH₂PO₄, 0.5 CaCl₂, 10 MgSO₄, titrated to a pH of 7.3–7.4 with HCl (Ting et al., 2014). Horizontal slices (400 μm) containing the thalamus were cut in ice-cold protective recovery solution using a vibratome (VT1200, Leica Biosystems, Wetzlar, Germany). Slices were trimmed to remove the hippocampus and the hypothalamus, and then transferred to protective recovery solution maintained at 32–34°C for 12 min. Brain slices were kept in room temperature ACSF consisting of the following (in mM): 126 NaCl, 2.5 KCl, 10 glucose, 26 NaHCO₃, 1.25 NaH₂PO₄, 2 CaCl₂, and 1 MgSO₄. All solutions were equilibrated with 95% O₂/5% CO₂.

Slices were placed in a humidified, oxygenated interface recording chamber and perfused with oxygenated ACSF (2 mL/min) at 32–34°C. 10 μM bicuculline was added to the ACSF (bicuculline-ACSF) to block GABA_A receptors. Oscillations were evoked by a square voltage pulse (10 V, 0.5 ms duration) delivered once every 60 s through two parallel tungsten electrodes (50–100 kΩ, FHC) 50–100 μm apart and placed in either the internal capsule or the reticular thalamus, which stimulated traversing corticothalamic and thalamocortical axons. Extracellular potentials were recorded in a differential manner with two tungsten electrodes (50–100 kΩ, FHC) by placing one in the somatosensory ventrobasal nuclei of the thalamus close to the stimulating electrode and one far away from the stimulating electrode. One experiment was performed per slice. Multi-unit recordings were amplified 10,000 times with a P511 AC amplifier (Grass), digitized at 10 kHz with Digidata 1440A, band-pass filtered between 100 Hz and 3 kHz, and acquired using Clampex 10.7 software (Molecular Devices, San Jose, CA).

After at least 20 min in baseline bicuculline-ACSF, slices were perfused with one of 4 possible solutions: (1) bicuculline-ACSF, (2) bicuculline-ACSF plus 4 μM NO-711 to block GAT1 transport (Sitte et al., 2002), (3) bicuculline-ACSF plus 100 μM SNAP-5114 to block GAT3 transport (Borden et al., 1994) or (4) bicuculline-ACSF plus a combination of 4 μM NO-711 and 100 μM SNAP-5114 to simultaneously block both GAT1 and GAT3. After 40 min, the perfusion solution was switched back to bicuculline-ACSF for at least 20 min.

Male Sprague-Dawley rats of postnatal day 11 to 15 (P11-P15) were used in dynamic clamp experiments, which were performed in accordance with protocols approved by the Administrative Panel on Laboratory Animal Care at Stanford University (Palo Alto, CA). Rats were deeply anesthetized with pentobarbital, then brains were rapidly extracted and placed in ice-cold protective recovery solution containing the following (in mM): 34 sucrose, 2.5 KCl, 11 glucose, 26 NaHCO₃, 1.25 NaH₂PO₄, 0.5 CaCl₂, and 10 MgSO₄, titrated to a pH of 7.4 with HCl. Horizontal slices (300 μm) containing the thalamus were cut in ice-cold protective recovery solution using a vibratome [VT1200 (RRID:SCR_018453), Leica Biosystems]. Slices were transferred to artificial cerebrospinal fluid (ACSF) maintained at 32°C for 45–60 min, then gradually brought to room temperature. The ACSF contained the following (in mM): 126 NaCl, 2.5 KCl, 10 glucose, 26 NaHCO₃, 1.25 NaH₂PO₄, 2 CaCl₂, and 1 MgSO₄. All solutions were equilibrated with 95% O₂/5% CO₂.

Slices were placed in a submerged recording chamber and perfused with oxygenated ACSF (2 mL/min) at 32–34°C. The chamber contained nylon netting which suspended the slice 1–2 mm from the chamber floor and enhanced slice perfusion. Slices were visualized with Dodt-contrast optics (Luigs and Newmann, Ratingen, Germany) on an Axioskop microscope (Zeiss, Pleasanton, CA). Recordings were obtained with a MultiClamp 700A patch amplifier (RRID:SCR_011323, Molecular Devices), digitized with Digidata 1322A and acquired using Clampex software. Borosilicate glass pipettes (1.5–3 MΩ) pulled on a P-87 micropipette puller (Sutter Instruments, Novato, CA) were filled with an internal solution containing (in mM): 100 potassium gluconate, 13 KCl, 10 EGTA, 10 HEPES, 9 MgCl₂, 2 Na₂-ATP, 0.5 Na-GTP, 0.07 CaCl₂ (pH 7.4).

Dynamic clamp experiments were conducted using a computer running the RealTime Application Interface for Linux (RTAI, www.rtai.org) sampling the intracellular potential at 50 kHz. Custom-written software modified from previous work (Sohal et al., 2006) used the sampled membrane potential and the pre-specified GABA_B-mediated inhibitory conductance to, at every timestep, update the inhibitory current that was injected into the thalamocortical neuron at 50 kHz. All measured voltages were corrected for a junction potential of −10 mV.

Inhibitory conductance waveforms (dynControl, dynGAT1-Block, dynGAT3-Block, dynDual-Block, Figure 2B) for use in dynamic clamp experiments were created based on previous recordings of GABA_B-mediated IPSCs (Beenhakker and Huguenard, 2010). In the previous study, voltage-clamped thalamocortical neurons were recorded during electrical stimulation of presynaptic reticular thalamic neurons while the thalamic slices were bathed in ionotropic glutamate receptor and GABA_A antagonists to isolate GABA_B-mediated IPSCs. After acquiring baseline data, one of three drugs was added to the perfusing ACSF: 4 μM NO-711 to block GAT1, 100 μM SNAP-5114 to block GAT3, or both 4 μM NO-711 and 100 μM SNAP-5114 to block both GAT1 and GAT3 (Borden et al., 1994; Sitte et al., 2002). After at least 10 min of drug perfusion and stabilization of drug effect, a series of GABA_B-mediated IPSCs was acquired. In the present study, we converted the previously obtained averaged IPSCs for each drug condition into a conductance waveform fitted with the equation:

where $A$ is the amplitude coefficient, ${\tau }_{rise}$ is the rise time constant, ${\tau }_{fallfast}$ and ${\tau }_{fallslow}$ are the fast and slow decay constants, respectively, and $w$ is the weighting factor for the two decay terms (Otis et al., 1993). The percent change in each parameter from baseline to drug condition was computed for each experiment. Baseline-normalized parameter values were compared between conditions, and for parameters with significant population differences in values (p<0.05, Wilcoxon signed rank test), the mean change in parameter value was found and implemented for that drug condition. For parameters that did not change significantly in the drug condition compared to baseline, the parameter value was kept identical to that of the baseline (dynControl) condition for the templates (Figure 2B). Parameter values for the four experiment conditions (dynControl, dynGAT1-Block, dynGAT3-Block, dynDual-Block) are listed in Table 2. A value of -115 mV was used for the reversal potential associated with the GIRK conductance following post-synaptic GABA_B receptor activation.

Since the number of GABA_B receptors activated in a physiological setting was not determined, all GABA_B IPSC waveforms (dynIPSCs) had amplitude scaled between 25–800%. A total of 5–15 repetitions were performed for each dynIPSC waveform, with a holding current adjusted at the beginning of each recording so the holding membrane potential is in the range of −73 to −60 mV. For each sweep, a brief current pulse (−50 pA, 10 ms) was applied at around 100 ms to assess electrode resistance and passive membrane properties, then the dynIPSC was applied at 1000 ms.

MATLAB R2018a (RRID:SCR_001622, MathWorks, Natick, MA) was used for all data analysis. Single action potential spikes were detected from raw multi-unit activity as follows (Sohal et al., 2003). The raw signal was first bandpass-filtered between 100 and 1000 Hz. Slopes between consecutive sample points were computed from the filtered signal. For each sweep $i$, the baseline slope noise ${\alpha }_i$ was computed from the root-mean-square average of the slope vector over the baseline region before stimulation start, and the maximum slope ${\beta }_i$ was defined as the maximum slope value at least 25 ms after stimulation start (to account for the stimulus artifact). For each slice, if $\alpha ={\sum }_i{\alpha }_i$ is the baseline slope noise averaged over all sweeps and $\beta ={\sum }_i{\beta }_i$ is the maximum slope averaged over all sweeps, then the slice-dependent signal-to-noise ratio was given by $r=1+0.1(\frac{\beta }{\alpha }-1)$. The resulting signal-to-noise ratios fell between 2 and 5. A sweep-dependent slope threshold was then defined by ${\theta }_i=r{\alpha }_i$. Finally, single spikes were defined as all local maxima of the slope vector at least 25 ms after stimulation start with values exceeding a slope threshold.

To compute the oscillation duration for each sweep, spikes were first binned by 10 ms intervals to yield a spike histogram. Evoked bursts were detected by joining consecutive bins with a minimum spike rate of 100 Hz, using a minimum burst length of 60 ms, a maximum delay after stimulation start of 2000 ms and a maximum inter-burst interval of 2000 ms (Figure 1A). The oscillation duration was defined as the time difference between the end of the last evoked burst and stimulation start.

To compute the oscillation period for each sweep, an autocorrelation function (ACF) was first computed from the binned spikes, then moving-average-filtered using a 100 ms window. Peaks (local maxima) were detected from the filtered ACF using a minimum peak prominence that is 0.02 of the amplitude of the primary peak (the first value of the ACF). Peaks with lags greater than the oscillation duration were ignored. The oscillation period was computed by searching for the lag value $\delta$ that minimizes the distance function $f\left(\delta \right)={\sum }_j|p_j-q_j(\delta \left)\right|$, where $p_j$ is the lag value of peak $j$ and $q_j\left(\delta \right)$ is the closest multiple of $\delta$ to $p_j$. The search was initialized with the estimate ${\delta }_0=p_2-p_1$ and confined to the bounds $\left[\frac{2}{3}{\delta }_0,\frac{3}{2}{\delta }_0\right]$.

An oscillatory index (Sohal et al., 2006) was also computed from the filtered autocorrelation function (fACF). For each sweep, the non-primary peak with largest fACF value was designated as the secondary peak. The oscillatory index was then defined by $OI=(A_{peak2}-A_{trough})/(A_{peak1}-A_{trough})$, where $A_{peak1}$ is the fACF value of the primary peak, $A_{peak2}$ is the fACF value of the secondary peak and $A_{trough}$ is the minimum fACF value between the primary peak and the secondary peak.

The average oscillation duration or period at the end of each phase (baseline or drug) was computed by choosing the five values from the last 10 sweeps of each phase that were within 40% of the average of the group of values. This approach was used to minimize bias caused by abnormally-shortened evoked oscillations due to the presence of spontaneous oscillations.

MATLAB R2018a was used for all data analyses. Responses to dynIPSCs were analyzed as follows: All voltage traces were manually examined and noisy recordings were excluded. The peak of the current trace (IPSC peak) within the first 300 ms of dynIPSC start was first detected. The most likely candidate for a calcium-dependent low-threshold spike (LTS) was then detected from the raw voltage trace in between the time of IPSC peak and 7000 ms after dynIPSC start.

To detect the LTS candidate, the voltage trace was first median-filtered with a time window of 30 ms to remove action potentials (Chung et al., 2002), then moving-average-filtered with a time window of 30 ms. First and second derivatives of the doubly-filtered voltage traces were computed by differences between consecutive sample points, and the first derivative vector was moving-average-filtered with a time window of 30 ms before taking the second derivative. The local maximum of the doubly-filtered voltage trace with the most negative second derivative was chosen as the LTS candidate. A histogram of LTS candidate second derivatives for traces was fitted to a sum of 3 Gaussian distributions and the minimum of the probability density function between the first two peaks was used as a second derivative threshold.

Features were computed for each LTS candidate. The LTS peak value was the absolute voltage value for the peak. The LTS latency was defined by the time difference between the LTS peak time and dynIPSC start. A time region that was bounded by the first local minimum of the doubly-filtered voltage trace on either side of the peak was used to compute the LTS maximum slope and detect action potentials. For accuracy, the slope value was computed from a different doubly filtered trace with a smaller moving-average-filter time window corresponding to roughly a 3 mV-change in amplitude for that trace. Action potentials spikes were detected by a relative amplitude threshold 10 mV above the LTS amplitude for that trace.

An LTS candidate was then assigned as an LTS if the following criteria were all satisfied: (1) Peak prominence must be greater than the standard deviation of the filtered voltage values before IPSC start; (2) Peak second derivative must be more negative than the threshold (−2.3 V/s²) described above; (3) If there were action potentials riding on the LTS, the LTS peak time must occur after the time of the first action potential. Detection results were manually examined by blinded experts and the algorithm’s decision for LTS determination was overturned only if three of the four polled electrophysiology experts agreed. The LTS or burst probability was defined as the proportion of traces producing an LTS or a burst (an LTS with at least one riding action potential) across all 5–15 repetitions (with varying holding potentials) for a particular GABA_B IPSC waveform.

All computational simulations were performed using NEURON (RRID:SCR_005393, Carnevale and Hines, 2009) version 7.5 with a temperature of 33°C. Each model thalamocortical neuron had one cylindrical somatic compartment (Soma) and two cylindrical dendritic compartments (Dend1 and Dend2) in series (Figure 3C). The somatic length and diameter were set to be equivalent with the parameter ${diam}_{soma}$. The two dendritic compartments were set to have equal diameters (${diam}_{dend}$), with each having lengths equal to half of the parameter $L_{dend}$. Values for specific membrane capacitance and axial resistivity were equivalent to those reported by Destexhe et al., 1998. Passive leak channels were inserted in all three compartments at equivalent densities ($g_{pas}$) with a fixed reversal potential of -70 mV. The 4 voltage-independent (passive) parameters described above were allowed to vary across neurons.

The following mechanisms were inserted in all three compartments: the T-type calcium current ($I_T$) and the submembranal calcium extrusion mechanism (${Ca}_{decay}$) was adapted from Destexhe et al., 1998; the hyperpolarization-activated nonspecific cationic current ($I_h$), the A-type transient potassium current ($I_A$), the inward-rectifying potassium current ($I_{Kir}$) and the persistent sodium current ($I_{NaP}$) were adapted from Amarillo et al., 2014. The following parameters were allowed to vary across compartments and across neurons: the maximum permeability ${\stackrel{-}{p}}_T$ (in cm/s) of $I_T$, which was described by the Goldman–Hodgkin–Katz flux equation (Huguenard and McCormick, 1992), and the maximum conductance densities ${\stackrel{-}{g}}_h$, ${\stackrel{-}{g}}_A$, ${\stackrel{-}{g}}_{Kir}$, ${\stackrel{-}{g}}_{NaP}$ (in S/cm²) of other currents described by Ohm’s law. All parameters that were not varied during optimization were identical for all model neurons and taken from literature values. Based on calculations from solutions used in experiments, the potassium reversal potential was -100 mV, the sodium reversal potential was 88 mV, the calcium concentration outside neurons was 2 mM and the initial calcium concentration inside neurons was 240 nM.

All intrinsic current mechanisms were described by Hodgkin-Huxley type equations with voltage and/or time-dependent activation ($m$) and inactivation ($h$) gating variables and have been described in detail in the corresponding sources (Amarillo et al., 2014; Destexhe et al., 1998). In particular, the open probability of $I_T$ was described by ${m\left(V,t\right)}^{2}h\left(V,t\right)$ (Huguenard and Prince, 1992). At every point in time, the activation variable $m$ converges exponentially to its steady-state value $m_{\infty }\left(V\right)$ with a time constant ${\tau }_m\left(V\right)$ and the inactivation variable $h$ converges exponentially to its steady-state value $h_{\infty }\left(V\right)$ with a time constant ${\tau }_h\left(V\right)$. The equations for $m_{\infty }\left(V\right)$, ${\tau }_m\left(V\right)$, $h_{\infty }\left(V\right)$ and ${\tau }_h\left(V\right)$ were fitted from electrophysiological recordings and described by Huguenard and McCormick, 1992.

For simulations involving action potentials, a fast sodium and potassium mechanism adapted from Sohal and Huguenard, 2003 was inserted in the somatic compartment and made identical across neurons. Based on the comparison of the number of spikes per LTS between model neurons and corresponding recorded neurons (Figure 4—figure supplement 1A–B), the threshold parameter $V_{Traub}$ was changed to −65 mV.

A custom GABA_B receptor mechanism was inserted in the somatic compartment that produced IPSCs with the equation form given by Equation 1. The parameters used were identical to that used for dynamic clamp, with kinetics varying for each pharmacological condition as in Table 2. Although parameter measurements were performed at a same temperature (33°C) as simulations, a Q₁₀ of 2.1 (Otis et al., 1993) was included in the model.

All simulations performed for each model neuron were matched to traces recorded using dynamic clamp for the corresponding recorded neuron as follows: To match the holding potentials, we first performed a test simulation to voltage clamp the model neuron (with initial potential −70 mV) for 2 s to bring it to a quasi-steady state. The resultant steady-state current was used as holding for the current clamp simulation. To simulate dynamic clamp of the recorded neuron, after 2 s to allow state variables in the model to stabilize, a brief current pulse (−50 pA, 10 ms) was applied at around 2.1 s, then an simIPSC identical to the dynIPSC applied in dynamic clamp was applied at 3 s. The analysis of LTS and burst features for simulated simIPSC responses is identical to that for recorded dynIPSC responses.

Initial values for passive parameters were different for each model neuron and were estimated from the current pulse responses from the corresponding recorded neuron by a strategy adopted from Johnston and Wu, 1994. A portion of raw current pulse responses had a systematic voltage shift at the beginning and end of the pulse, consistent with an unbalanced bridge. These shifts were detected by a slope threshold determined through a histogram of all initial slopes and then corrected by shifting the entire portion of the response during the brief current pulse, by the calculated amount. The corrected current pulse responses for a particular neuron were then fitted to the following first order response equation with two exponential components:

where $C_0$, $C_1$ are the amplitudes (in mV) of the two components, ${\tau }_0$, ${\tau }_1$ are the time constants (in ms) of the two components, $t_p=10\mathrm{}\mathrm{m}\mathrm{s}$ is the width of the current pulse, and $B\left(x\right)$ is a Boolean function defined by $B\left(x\right)=1$ if $x=\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}$ and $B\left(x\right)=0$ if $x=\mathrm{f}\mathrm{a}\mathrm{l}\mathrm{s}\mathrm{e}$. Initial values for the curve fit were $C_{0,i}=∆V$, $C_{1,i}=∆V$, ${\tau }_{0,i}=10\mathrm{m}\mathrm{s}$, ${\tau }_{1,i}=1\mathrm{m}\mathrm{s}$ where $∆V$ is the mean voltage change for the neuron after each 10 ms stimulus.

Next, given the fitted coefficients, we estimated the defining parameters for a ball-and-stick model. The input resistance was computed by

where $I_p=50\mathrm{}\mathrm{p}\mathrm{A}$ is the amplitude of the current pulse. The membrane time constant ${\tau }_m$ was set to equal that of the slow component ${\tau }_m={\tau }_0$. The length constant was computed by solving for $L$ in equation 4.5.57 of Johnston and Wu:

where ${\alpha }_1={\left(\left({\tau }_1/{\tau }_0\right)-1\right)}^{1/2}$, starting with the initial guess of $L_i=\pi /{\alpha }_1$. The dendritic-to-somatic conductance ratio $\rho$ was then computed by equation 4.5.58 of Johnston and Wu:

Finally, we converted the defining ball-and-stick parameters to initial estimates for the passive parameters in our 3-compartment NEURON model that are optimized. The input resistances of the soma and the dendrite were computed by

where $R_s$ is the series resistance. The specific membrane resistivity $R_m$ was computed by equation 4.3.3 of Johnston and Wu:

where a fixed value of 0.88 µF/cm² (Destexhe et al., 1998) for the specific membrane capacitance $C_m$. Then the somatic diameter ${diam}_{Soma}$ was computed by equation 4.3.8 of Johnston and Wu:

the dendritic diameter ${diam}_{Dend}$ was computed by equation 4.5.47 of Johnston and Wu:

where a fixed value of 173 Ω·cm (Destexhe et al., 1998) for the axial resistivity $R_a$, the dendrite length $L_{Dend}$ was computed by equation 4.4.15 of Johnston and Wu:

and the passive conductance $g_{pas}$ was estimated by

where $A=\pi \left({{diam}_{Soma}}^2+{diam}_{Dend}{L}_{Dend}\right)$ is the total surface area of the neuron.

Initial values for active parameters were identical for all model neurons and were taken from literature values (Amarillo et al., 2014; Destexhe, 1998) as follows: ${\stackrel{-}{p}}_{T,Soma}={\stackrel{-}{p}}_{T,Dend1}={\stackrel{-}{p}}_{T,Dend2}=2.0\mathrm{x}\mathrm{}{10}^{-4}\mathrm{c}\mathrm{m}/\mathrm{s}$, ${\stackrel{-}{g}}_{h,Soma}={\stackrel{-}{g}}_{h,Dend1}={\stackrel{-}{g}}_{h,Dend2}=2.2\mathrm{x}\mathrm{}{10}^{-5}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$, ${\stackrel{-}{g}}_{Kir,Soma}={\stackrel{-}{g}}_{Kir,Dend1}={\stackrel{-}{g}}_{Kir,Dend2}=2.0\mathrm{x}\mathrm{}{10}^{-5}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$, ${\stackrel{-}{g}}_{A,Soma}={\stackrel{-}{g}}_{A,Dend1}={\stackrel{-}{g}}_{A,Dend2}=5.5\mathrm{x}\mathrm{}{10}^{-3}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$ and ${\stackrel{-}{g}}_{NaP,Soma}={\stackrel{-}{g}}_{NaP,Dend1}={\stackrel{-}{g}}_{NaP,Dend2}=5.5\mathrm{x}\mathrm{}{10}^{-6}\mathrm{S}/{\mathrm{c}\mathrm{m}}^2$.

The voltage responses of a model neuron to a set of 12 different GABA_B IPSC waveforms (four pharmacological conditions at three different conductance amplitude scales) was compared with the responses of the corresponding recorded neuron. To emphasize the LTS response, model neurons excluded fast sodium and potassium channels and the experimental responses were median-filtered with a time window of 30 ms to remove action potentials. As multiple objective functions were of interest (Druckmann et al., 2007), a total error $E_{total}$ was computed by a weighted average of component errors (Figure 3E):

where $E_m$ is the average error across all traces for whether the presence or absence of LTS matches (if an LTS is missed an error of 18 is assigned; if an LTS is falsely produced an error of 6 is assigned), $E_{sw}$ is the average root-mean-square error of voltage values, $E_a$ is the average LTS peak value error across all LTS-matching traces, $E_t$ is the average LTS latency error across all LTS-matching traces, $E_{sl}$ is the average LTS maximum slope error across all LTS-matching traces, and the $w_i$’s denote the corresponding weights. Based on the maximum LTS latency observed in recorded neurons, the fitting region was restricted to between 0 and 1.8 s after the start of the IPSC. Traces were also allowed to have different sweep weights, in which case the average across traces were computed using root-mean-square.

A modified version of the Nelder-Mead simplex algorithm (Lagarias et al., 1998) was used to update parameter values within defined parameter bounds, with the objective of minimizing $E_{total}$. For each simplex, the maximum number of iterations was 2000, the maximum number of error evaluations was 4000, the relative error tolerance was 0.1, the relative parameter change tolerance was 0.1, the coefficient of reflection was $\rho =1$, the coefficient of expansion was $\chi =2$, the coefficient of contraction was $\gamma =0.75$, the coefficient of shrinkage was $\sigma =0.8$. The initial simplex was set up so that each vertex deviates in a parameter to be varied by $\delta =2/3$ of the parameter’s total range. To transform values to an unconstrained parameter space so that the Nelder-Mead method could be applied, each bounded parameter value was first linearly transformed to the region [−1, 1] than transformed to $(-\infty ,\infty )$ using the inverse tangent function.

A total of 21 different iterations were run to optimize each model neuron, with the weights among different error types and weights among different traces varying from iteration to iteration. For some poorly fit neurons, the initial parameters for that neuron were replaced by the best-fit neuron’s parameters from the previous iteration. One trace per GABA_B IPSC waveform was used during model optimization, but all of 60–180 traces were used for the final model evaluation (Figure 3E).

Each 2 cell network includes one single-compartment model reticular thalamic neuron, described in previous work (Klein et al., 2018), and one of the well-fitted, three-compartment model thalamocortical neurons described earlier. The reticular thalamic neuron has an AMPA receptor that is synaptically activated by the thalamocortical neuron, and the thalamocortical neuron has a GABA_B receptor in the somatic compartment that is synaptically activated by the reticular thalamic neuron. Synaptic currents were evoked with 100% probability and a delay of 1 ms whenever the presynaptic neuron reached a voltage threshold of −30 mV.

AMPA receptors were adapted from Sohal et al., 2000. Reflecting more recent physiological measurements (Deleuze and Huguenard, 2006), AMPA currents were adjusted to bring rise time to 0.5 ms, decay time to 5.6 ms and maximal conductance of 7 nS per synapse. The reversal potential for AMPA currents was maintained at 0 mV.

GABA_B receptors were as described before for single model TC neurons. In order for overlapping IPSCs to be summed linearly, we expanded Equation 1 into 18 terms, each having its own rise and decay exponential time constants. To be consistent with dynamic clamp experiments, a reversal potential of −115 mV was used in simulations, although the networks behave similarly if a reversal potential of −100 mV was used instead (data not shown). Relative to the values used in dynamic clamp, the conductance amplitudes were scaled by 1/12 to account for temporal summation from an RT burst, and the synaptically-evoked IPSCs was verified to be comparable to recorded values (Figure 8B).

Each 200 cell network assumed a bilayer architecture, with one 100 cell circular layer of single-compartment model reticular thalamic neurons and one 100 cell circular layer of three-compartment model TC neurons. The identity of the TC neuron was chosen from each of the 31 well-fitted model TC neurons. Each reticular thalamic neuron inhibited the nine nearest TC cells, whereas each TC neuron excited the five nearest reticular thalamic neurons (Sohal and Huguenard, 2003). The passive leak conductance of each model reticular thalamic neuron was randomly selected from a uniform distribution between 45 and 55 µS/cm² (Sohal and Huguenard, 2003). The passive leak conductance of each model TC neuron was randomly selected from a uniform distribution between −10 and 10% of the optimized value.

Network simulations were performed with leak reversal potentials and initial membrane potentials set to a value between −73 mV and −60 mV, at 1 mV increments, so that a total of 14 repetitions were applied for each model network. After a delay of 3 s to allow for state variable stabilization, either the reticular thalamic neuron in the two-cell network or each of the center 20 reticular thalamic neurons in the 200 cell network was injected with a square current pulse (0.2 nA, 40 ms). Simulations were performed with a 0.1 ms integration time step and continued for a total of 30 s. Pooling all action potential spikes after stimulation end, we detected bursts and computed an oscillation period and an oscillatory index using the same algorithm as that for multiunit recordings, except for a bin width of 100 ms for spike histograms and a minimum peak prominence of 0.5 relative to the largest secondary peak for the filtered autocorrelation function. Oscillation probability was defined as the proportion of simulations (with varying holding potentials) that induced oscillations with at least three bursts. Percent of active TC cells was defined as the percentage of model TC neurons in the network that produced at least one spike. Half activation latency was defined as the time it took for half of the final percentage of active cells to be activated. Mean oscillation period, mean oscillatory index and mean half activation latency were computed by restricting to trials with a successfully evoked oscillation.

MATLAB R2019b was used for all statistical analysis. Since the number of available data points for LTS and burst features was significantly different for the Dual-Block condition, we performed a paired t-test or a signed-rank test between the Control condition and the Dual Block condition. We used either repeated-measures ANOVA or the Friedman’s test for comparison across the Control, GAT1-Block and GAT3-Block conditions.

Paired comparisons were applied if not otherwise specified. Normality of the differences to the within-subject mean was assess for all groups using a combination of the Lilliefors test, the Anderson-Darling test and the Jarque-Bera test. Normality was satisfied when the geometric mean of three p values was at least 0.05 for all groups. When normality was satisfied, the paired-sample t-test was used when there are two groups and repeated-measures ANOVA (with multiple comparison) was used when there are more than two groups. When normality was not satisfied, the Wilcoxon signed-rank test was used when there are two groups and Friedman’s test (with multiple comparison) was used when there are more than two groups. Tests were two-tailed with a significance level of 0.05. Error bars reflect 95% confidence intervals. All violin plots used a bandwidth that was 10% of the maximum data range.

Bicuculline methiodide was purchased from Sigma-Aldrich (St. Louis, MO). The GAT1 blocker NO-711 [(1,2,5,6-tetrahydro-1-[2-[[(diphenylmethylene)amino]oxy]ethyl]-3-pyridinecarboxylic acid hydrochloride)] and GAT3 blocker SNAP-5114 [1-[2-[tris(4-methoxyphenyl)methoxy]ethyl]-(S)-3-piperidinecarboxylic acid] were purchased from Tocris Bioscience (Minneapolis, MN).

Data analysis, data visualization, statistical tests, NEURON model setup and NEURON model optimization were performed using custom MATLAB code. All code for reproducing results (including NEURON. hoc,. tem and. mod files used) are available online at https://github.com/luadam4c/m3ha_published/. (Lu, 2020; copy archived at https://github.com/elifesciences-publications/m3ha_published).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank the University of Virginia Medical Scientist Training Program, the Whitfield-Randolph Scholarship and the NIH grants R01-NS099586 and R01-NS034774 for funding support. We thank Shin-Shin Nien for help with figure organization. We thank Dr. Peter Klein, Dr. Lise Harbom and Katie Salvati for help with manual LTS scoring. We thank Dr. Farzan Nadim for providing guidance on computational modeling.

Animal experimentation: Oscillation experiments were performed in accordance with Protocol #3892 approved by the Institutional Animal Care and Use Committee at the University of Virginia. Dynamic clamp experiments were performed in accordance with protocols approved by the Administrative Panel on Laboratory Animal Care at Stanford University. Rats were deeply anesthetized with pentobarbital before transcardial perfusion, and every effort was made to minimize suffering.

© 2020, Lu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.