1. Computational and Systems Biology
  2. Neuroscience
Download icon

Nonlinearities between inhibition and T-type calcium channel activity bidirectionally regulate thalamic oscillations

  1. Adam C Lu  Is a corresponding author
  2. Christine Kyuyoung Lee
  3. Max Kleiman-Weiner
  4. Brian Truong
  5. Megan Wang
  6. John R Huguenard  Is a corresponding author
  7. Mark P Beenhakker  Is a corresponding author
  1. Department of Pharmacology, University of Virginia, United States
  2. Department of Neurosurgery, Massachusetts General Hospital, United States
  3. Department of Psychology, Harvard University, United States
  4. Princeton Neuroscience Institute, Princeton University, United States
  5. Department of Neurology, Stanford University, United States
Research Article
  • Cited 0
  • Views 980
  • Annotations
Cite this article as: eLife 2020;9:e59548 doi: 10.7554/eLife.59548

Abstract

Absence seizures result from 3 to 5 Hz generalized thalamocortical oscillations that depend on highly regulated inhibitory neurotransmission in the thalamus. Efficient reuptake of the inhibitory neurotransmitter GABA is essential, and reuptake failure worsens human seizures. Here, we show that blocking GABA transporters (GATs) in acute rat brain slices containing key parts of the thalamocortical seizure network modulates epileptiform activity. As expected, we found that blocking either GAT1 or GAT3 prolonged oscillations. However, blocking both GATs unexpectedly suppressed oscillations. Integrating experimental observations into single-neuron and network-level computational models shows how a non-linear dependence of T-type calcium channel gating on GABAB receptor activity regulates network oscillations. Receptor activity that is either too brief or too protracted fails to sufficiently open T-type channels necessary for sustaining oscillations. Only within a narrow range does prolonging GABAB receptor activity promote channel opening and intensify oscillations. These results have implications for therapeutics that modulate inhibition kinetics.

Introduction

Neural circuits rely on a combination of intrinsic cellular properties and synaptic connections to generate large-scale electrical oscillations that drive behavior (Getting, 1989; Marder and Calabrese, 1996; Nusbaum and Beenhakker, 2002; Huguenard and McCormick, 2007). Following membrane hyperpolarization, such as that produced by synaptic inhibition, cortically projecting neurons of the thalamus [i.e. thalamocortical (TC) neurons] produce brief bursts of action potentials (Llinás and Jahnsen, 1982), a cellular property that maintains both sleep-related and seizure-related oscillations (McCormick and Contreras, 2001; Huguenard and McCormick, 2007; Beenhakker and Huguenard, 2009). Several studies have shown that CaV3.1 T-type calcium channels (T channels) sustain post-inhibitory rebound bursts in thalamocortical neurons by producing a relatively prolonged calcium-dependent, low-threshold spike (Kim et al., 2001; Kim et al., 2003; Porcello et al., 2003). These channels require membrane depolarization to open and hyperpolarization to recover (Coulter et al., 1989). Hyperpolarization-dependent recovery involves the removal of T channel inactivation (i.e. de-inactivation). As T channels are largely inactivated at resting membrane potentials, membrane hyperpolarization is necessary for robust rebound bursting (Llinás and Jahnsen, 1982; Coulter et al., 1989). While controlled voltage-clamp experiments have informed our understanding of how neuronal membrane potential dynamics can affect T channel opening (Gutierrez et al., 2001), we still know little regarding channel behavior during physiological forms of synaptic inhibition when the underlying conductance rises and falls nonlinearly.

Reticular thalamic (RT) neurons serve as the main source of inhibitory, GABAergic input to thalamocortical neurons, especially in rodents (Shosaku, 1985; Pinault and Deschênes, 1998). Thalamocortical neurons express synaptic α1β2γ2 GABAA receptors, and two types of extrasynaptic receptors: GABAA4β2δ) and GABAB (Pirker et al., 2000; Kulik et al., 2002; Jia et al., 2005). Studies have shown that modulation of synaptic GABAA receptors on TC neurons has little effect on thalamocortical oscillations (Sohal et al., 2003; Rovó et al., 2014). In contrast, extrasynaptic receptors have been implicated in thalamocortical seizures, both for GABAA (Cope et al., 2009) and GABAB (Liu et al., 1992; Vergnes et al., 1997; Bortolato et al., 2010) receptors. Prior experimental and computational work has demonstrated that a shift from GABAA receptor-mediated to GABAB receptor-mediated inhibition at the RT-TC synapse transforms oscillations 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).

GABA transporters (GATs) powerfully control the activation of GABAB receptors (Beenhakker and Huguenard, 2010). GAT1 and GAT3 represent the primary GABA transporters expressed in the brain and normally recycle GABA from the extrasynaptic space, thereby regulating GABA spillover from the synapse and the activation of extrasynaptic GABAA and GABAB receptors (Cope et al., 2009; Scanziani, 2000). In the thalamus, the more abundant transporter, GAT3, is localized farther away from synapses than GAT1 (De Biasi et al., 1998; Beenhakker and Huguenard, 2010). Consequently, specific GAT1 blockade results only in an increase in the amplitude of the GABAB IPSC, reflecting increases in GABA concentration near the synapse. In contrast, specific GAT3 blockade results in an increase in both amplitude and decay of the GABAB-mediated inhibitory post-synaptic current (GABAB IPSC), as GABA is allowed to diffuse far from the synapse where there is an abundance of GABAB receptors (Kulik et al., 2002; Beenhakker and Huguenard, 2010). On the other hand, dual GAT1 and GAT3 blockade results in a roughly 10-fold increase in the decay time constant of the GABAB IPSC, a supralinear effect relative to the modest increases for single GAT1 (no increase) or GAT3 (1.5-fold) blockade. These findings were replicated in a diffusion-based computational model (Beenhakker and Huguenard, 2010).

In this study, we investigate the consequences of physiologically relevant GABAB receptor-mediated inhibition observed during different combinations of GABA transporter blockade: control, GAT1 blockade, GAT3 blockade and dual GAT1+GAT3 blockade (Beenhakker and Huguenard, 2010). As absence seizures are dependent on GABAB receptor signaling, we hypothesized that GAT blockade would regulate both thalamocortical neuron rebound bursting and network-level thalamic oscillations. We examine the effects of different GABAB receptor activation waveforms on both absence seizure-like thalamic oscillations and single thalamocortical neuron responses. We first use pharmacological manipulations to demonstrate that individual GAT1 or GAT3 blockade prolongs seizure-like oscillations, but that dual GAT1+GAT3 blockade surprisingly abolishes oscillations. Next, we apply physiological GABAB IPSC waveforms corresponding to each pharmacological condition to single thalamocortical neurons with dynamic clamp and demonstrate that individual GAT1 or GAT3 blockade increases rebound burst probability, but that dual GAT1+GAT3 blockade suppresses it. We then build computational model neurons to explore how differential GABAB IPSC kinetics shape rebound bursting in TC neurons. These models show how prolonged inhibition that pushes instantaneous T channel inactivation states toward steady-state values is responsible for rebound burst failure, and thus the nonlinear effects of rebound bursting associated with dual GAT1+GAT3 blockade. Finally, we build computational thalamic network models to demonstrate that the same interplay between GABAB-mediated inhibition and T channels is sufficient to explain the effects of GAT blockade on biological thalamic oscillations. Through these experimental and computational approaches, we identify how GABAB-mediated inhibition across both voltage and time dimensions regulates T channel activity and seizure-like oscillations.

Results

Thalamic oscillations

To evaluate the contribution of GABA transporters to thalamic network activity in the context of GABAB 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 GABAA 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.

Individual GAT1 or GAT3 blockade strengthens thalamic oscillations, but dual GAT1+GAT3 blockade abolishes oscillations.

(A) Slice recording setup and sample analysis. Acute thalamic slices were bathed in bicuculline to block GABAA receptors. A brief voltage stimulus (0.5 ms, 10 V) was applied with a bipolar electrode placed in either the reticular thalamic nucleus or the adjacent internal capsule to evoke epileptiform oscillations recorded extracellularly in the ventrobasal complex. Spikes were detected, binned, and grouped into bursts. The oscillation duration was measured from the spike histogram and the oscillation period was computed from the autocorrelation function of binned spikes. (B) Example evoked epileptiform oscillations at baseline and 40 min after perfusing with (i) control (no drug added), (ii) 4 µM NO-711 (GAT1 blocker), (iii) 100 µM SNAP-5114 (GAT3 blocker) or (iv) dual 4 µM NO-711+100 µM SNAP-5114. (C) Example PSTHs over the entire course of a single control, GAT1-, GAT3-, and dual-block experiment. Oscillations were evoked every 60 s, but only the first 17 s after stimulation are shown. Drugs were perfused for 40 min after 20 min of baseline, followed by 20 min of washout. (D–F) Oscillation measures for all slices. Colored circles denote the mean value, colored lines denote the mean change; and shaded areas denote the 95% confidence intervals for the mean change. (D) Oscillation duration did not change following control perfusion, but increased following NO-711 or SNAP-5114 perfusion. Dual NO-711+SNAP-5114 perfusion reduced oscillation duration (*p<0.05, **p<0.01, paired-sample t-test). Inset shows mean change relative to baseline. (E) Oscillation period did not change following control perfusion and lengthened following NO-711 or SNAP-5114 perfusion (**p<0.01, paired-sample t-test). (F) After 60 min of dual NO-711 + SNAP-5114 perfusion, oscillations were abolished in all slices (**p<0.01, paired-sample t-test).

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

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 GABAB 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 GABAB 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 GABAB 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 GABAB 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 GABAB 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.

Single neuron recordings

We investigated the effects of GAT-modulated, GABAB 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 GABAB 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 GABAB 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
Post-inhibitory, low-threshold rebound spikes and bursts in thalamocortical neurons are bidirectionally modulated by GABAB receptor-mediated conductance waveforms.

(A) Dynamic clamp setup. A thalamocortical neuron was patched in the whole-cell configuration. The applied current was computed from the instantaneous voltage and a command conductance waveform over time to simulate GABAB receptor activation. (B) Command GABAB receptor conductance waveforms (dynIPSCs, amplitudes scaled by 200%) that emulated different GAT blockade conditions based on GABAB IPSCs isolated with voltage clamp (Beenhakker and Huguenard, 2010). (C) Sample voltage responses of two neurons to the four different dynIPSCs shown in (B). (D) Distributions of post-inhibitory rebound burst measures over all 47 recorded neurons across dynIPSCs shown in (B). For burst latency, only responsive neurons are included. Relative to dynControl responses, rebound burst probability increased following either dynGAT1- or dynGAT3-Block, but decreased following dynDual-Block (*p<0.05, **p<0.01, ***p<0.001, Friedman’s test for burst probability, repeated-measures ANOVA and paired-sample t-test for burst latency). (E) Mean burst measures over all 47 recorded neurons, across four different dynIPSC waveforms and six different conductance amplitude scales. Error bars denote 95% confidence intervals.

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

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−7), but lower following dynDual-Block (n = 47 cells, LTS: −88%, p=4.8×10−6, 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−9) and burst latency (+58%, n = 21 cells, p=1.1×10−9). 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).

Table 1
Average change of LTS and burst features relative to dynControl responses.

For LTS and burst probability, Friedman’s test with multiple comparison was used across all four groups. Due to the lack of LTS response to dynDual-Block, the comparison between dynControl and dynDual-Block for all other LTS and burst features was conducted separately using the paired-sample t-test. Across dynControl, dynGAT1-Block and dynGAT3-Block responses, Friedman’s test with multiple comparison was used for LTS latency and repeated-measures ANOVA with multiple comparison was used otherwise.

dynGAT1-BlockdynGAT3-BlockdynDual-Block
Burst Probability+63%, p=0.0018+106%, p=1.7×10−7−82%, p=0.030
Burst Latency+4.6%, p=0.034+58%, p=1.1×10−9no change, p=0.066
LTS Probability+26%, p=0.0080+39%, p=0.0015−88%, p=4.8×10−6
LTS Latencyno change, p=0.97+53%, p=3.7×10−9+254%, p=0.044
Spikes Per LTS+62%, p=1.3×10−6+93%, p=9.3×10−7no change, p=0.095
LTS Peak Value (mV)+2.5 ± 0.5 mV, p=1.4×10−4+3.3 ± 0.8 mV, p=5.2×10−4no change, p=0.33
LTS Maximum Slope+52%, p=2.6×10−9+82%, p=1.4×10−9no change, p=0.068

In summary, the bidirectional differences in thalamocortical neuron rebound burst properties in response to different GABAB 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 GABAB receptors.

Single neuron models

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
Model thalamocortical neurons reproduce GABAB IPSC and rebound responses.

(A) Model optimization workflow. (B) Sample double-exponential curve fits (red) to averaged current pulse responses (blue) for an example neuron. The dashed lines correspond to the curves representing the somatic compartment (green) and dendritic compartment (orange). (C) Ball-and-stick geometries estimated from the two exponential components shown in (B) (Johnston and Wu, 1994) were converted to cylindrical, three compartment models (left), which were then optimized (right). (D) Fits of simulated simIPSC responses to recorded dynIPSC responses, for the same example neuron. (E) The 33 model neurons that underwent optimization were ranked by a weighted average of five different types of errors (see Materials and methods). The 31 highest ranked model neurons were considered well-fitted. The example neuron in (D) had rank 6. (F) Final values of parameters that could vary for the 31 well-fitted model neurons. Note that the T channel density is invariably high in the dendrites and the A-type potassium channel density is invariably high in the soma among all model neurons (red boxes). Geometric parameters are in µm, maximal conductance densities (g-) and conductance densities (g) are in S/cm2 and maximal permeability densities (p-) are in cm/s. The x axis is the model neuron rank in (E). We distinguish among (i) geometric parameters, (ii) voltage-independent passive parameters and (iii) voltage-dependent active parameters.

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

Our preliminary modeling results using existing TC cell models (Destexhe et al., 1998; Amarillo et al., 2014) failed to recapitulate two key features of GABAB 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 (diamsoma), the dendritic diameter (diamdend), the dendritic length (lengthdend) and the passive leak conductance density (gpas). 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 (IT), the hyperpolarization-activated nonspecific cationic channel (Ih), the A-type transient potassium channel (IA), the inward-rectifying potassium channel (IKir) and the persistent sodium channel (INaP). 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 (diamSoma, diamDend, lengthdend, gpas) 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. GABAB 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 GABAB 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 GABAB 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−4, dynamic clamp: +45%, p=0.011) and decreased following sim/dynDual-Block (simulation: −75%, p=9.8×10−7, dynamic clamp: −93%, p=5.7×10−4). 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−8, dynamic clamp: +52%, n = 21, p=2.4×10−6). 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
Well-fitted model and recorded neurons show similar low-threshold rebound spike probabilities and latencies in response to different GABAB IPSC waveforms.

(A) Distributions of post-inhibitory, low-threshold rebound spike measures over the 31 well-fitted model neurons across GABAB IPSC waveforms shown in Figure 2B (*p<0.05, **p<0.01, ***p<0.001, repeated-measures ANOVA for LTS probability, Friedman’s test for LTS latency). (B) Same as (A) but for the corresponding 31 recorded neurons (Friedman’s test). (C) Mean low-threshold rebound spike measures over all 31 model neurons, across four different GABAB IPSC waveforms and three different conductance amplitude scales. Error bars denote 95% confidence intervals. (D) Same as (C) but for the corresponding 31 recorded neurons.

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

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 GABAB-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.

Interplay between GABAB receptors and T-type calcium channels

We next sought to understand the underlying mechanisms contributing to the different burst responses following different GAT-modulated, GABAB-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 (mT) and inactivation (hT) gating variables of the T channel, as a function of time (Figure 5B). By convention, mT=1 when all channels are activated, and hT=0 when all channels are inactivated (Hodgkin and Huxley, 1952). Also by convention, T channel open probability is given by mT2hT (Huguenard and McCormick, 1992). We also distinguished between steady-state values [mT,(V) and hT,(V)] that depend only on voltage, from instantaneous values [mT(V,t) and hT(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 [τm(V) and τh(V)].

Figure 5 with 1 supplement see all
T-type calcium channel open probability discrepancy predicts LTS production following GABAB IPSC waveforms.

(A) The LTS response was correlated with the presence of large T-type calcium currents. (i) Command sim/dynIPSCs as in Figure 2B. (ii) Voltage responses of Neuron 1 of Figure 2C recorded using dynamic clamp. (iii-v) Simulated responses of the corresponding model neuron (same as Figure 3D), including: (iii) somatic voltage, (iv) GABAB receptor current, (v) T current. Currents were summed over all three compartments. (B) A combination of T channel activation (high mT), T channel recovery (high hT) and inactivation lag (hT different from hT,) appeared to be necessary for T channel opening (mT2hT much higher than mT,2hT,). State variables for the distal dendritic T channel (other two compartments are similar), including: (i) instantaneous activation gating variable, (ii) steady-state activation gating variable, (iii) instantaneous inactivation gating variable, (iv) steady-state inactivation gating variable, (v) instantaneous open probability (solid line) versus steady-state open probability (dotted line), (vi) difference of instantaneous versus steady-state open probability (open probability discrepancy). Note that only positive discrepancies are plotted. The green dotted line here [also in (C) and (D) (iii) ] denotes the observed LTS threshold. (C) The maximum open probability discrepancy was higher in LTS-producing responses (***p<0.001, n = 31 cells, paired-sample t-test). (D) Within the first voltage peak following hyperpolarization highlighted in both (i) somatic voltage curves and (ii) dendritic T channel open probability discrepancy curves, LTS-producing responses followed trajectories different from LTS-lacking responses in either (iii) the voltage vs. open probability discrepancy phase plot or (iv) the concavity versus slope of discrepancy phase plot. The slope and concavity are the first and second derivatives of log10(open probability discrepancy), respectively. Sample points plotted are 1 ms apart. Circles denote the decision points where either zero (or maximal negative) concavity is reached in the open probability discrepancy curves. Labels a, b, c, d are explained in the text.

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

One notable feature of the T channel is that the activation time constant τm is rapid, indicating that the instantaneous activation variable mT quickly achieves its voltage-dependent steady-state value mT,. Thus, throughout the course of the relatively slow membrane voltage changes associated with the GABAB-receptor-mediated IPSC, mT was nearly identical to mT, [Figure 5B(i–ii)]. By contrast, the T channel inactivation time constant τh is slow, with a value about 10-fold higher than τm (Coulter et al., 1989). Thus, during the IPSC, the instantaneous T channel inactivation variable hT rarely achieved its voltage-dependent steady-state value hT, [Figure 5B(iii) and (iv)].

At each time point, the capacity of hT to achieve hT, can be quantified by calculating the difference between the two values. If the discrepancy between the two values is low (i.e. hT-hT,0), then hT has successfully achieved its steady-state value at that time point. A high discrepancy, on the other hand, reflects an hT that has incompletely reached steady state. Importantly, the magnitude of discrepancy between hT and hT, affected the T channel open probability mT2hT, but only when the activation variable mT was high. We therefore quantified the difference between instantaneous versus steady-state open probability (mT2hT-mT,2hT,), 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 [maxt(mT2hTmT,2hT,)] 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−29). In fact, when all traces were considered, a threshold open probability discrepancy of 10−2 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 hT (which was high from the hyperpolarization) slowly approached low hT, (Figure 5D(ii),(iii): label a); discrepancy was largely defined by the slow kinetics of hT because mT matches mT, 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
Voltage and T-type calcium channel open probability discrepancy trajectories in response to simGAT3-Block and simDual-Block GABAB IPSCs.

Simulated traces for the example neuron of Figure 3, resampled at 1 ms intervals, in response to either simGAT3-Block (red) or simDual-Block (purple), as shown in Figure 5. Traces over simulation time are shown for (A) somatic voltage, (B) distal dendritic T channel activation gating variable, (C) distal dendritic T channel inactivation gating variable and (D) distal dendritic T channel open probability discrepancy. (E) Phase plot of somatic voltage versus distal dendritic T channel open probability discrepancy. (F) Phase plot of the concavity versus slope of the distal dendritic T channel open probability discrepancy. Note that only positive discrepancy values are plotted. All traces are aligned in Video time to the decision points (circles), defined as the point at which the open probability discrepancy curve either reaching an inflection point (zero concavity) or maximal negative concavity. Note that the traces for instantaneous and steady-state activation gates in (B) overlap on this time scale. In (E-F), only the first voltage peak following hyperpolarization is shown. The time points shown correspond to those with larger marker sizes in (A-D). Following either simGAT3-Block or simDual-Block, the open probability discrepancy curve reaches an inflection point with zero concavity. However, only following simGAT3-Block does the open probability discrepancy curve reach a significantly positive concavity, drive past the LTS threshold (10−2) and produce a rebound low-threshold spike (LTS). Note that the slope of the open probability discrepancy curve at the decision point is higher for simGAT3-Block than for simDual-Block.

Video 2
Voltage and T-type calcium channel open probability discrepancy trajectories in response to all pharmacological GABAB IPSCs.

Same as Video 1 but in response to all simIPSCs shown in Figure 5. (A-F) See descriptions for Video 1. Note that following simControl (black), the maximum concavity for the open probability discrepancy curve is negative, so its trajectory in the open probability discrepancy concavity versus slope phase plot (F) never reaches above the y = 0 line. For other simIPSCs, note that the slope of the open probability discrepancy curve at the decision point is higher for LTS-producing trajectories [simGAT1-Block (blue) and simGAT3-Block (red)] than for an LTS-lacking trajectory [simDual-Block (purple)].

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 (τh) for the same simIPSC response simulations as in Figure 5. When τ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 τ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.

Manipulation of the T-type calcium channel inactivation time constant modulates LTS production.

(A) Simulated responses as in Figure 5 with the T channel inactivation time constant τhT halved. The T channel open probability discrepancy failed to reach threshold and no LTS was produced following any simIPSC waveform. The green dotted line denotes the observed LTS threshold from Figure 5C. In the voltage vs. open probability discrepancy phase plot, only the first voltage peak following hyperpolarization is shown. (B) Simulated responses as in Figure 5 with τhT doubled. The T channel open probability discrepancy increased and LTSs appeared following simDual-Block. Note that there is still no LTS produced following simControl.

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
Manipulation of GABAB IPSC kinetics bidirectionally modulates low-threshold rebound spike production.

(A) (1) LTS responses to the simGAT3-Block waveform (blue) gradually disappeared as the time constant was increased (with amplitude fixed) to that of the simDual-Block waveform (red). The green dotted line denotes the observed LTS threshold from Figure 5C. (2) LTS responses to the sDual-Block waveform (red) gradually appeared as the time constant of the waveform was decreased (with amplitude fixed) to that of the sGAT3-Block waveform (blue). (B) (1) LTS responses gradually appeared as the amplitude of the simGAT3-Block waveform was increased (with time constant fixed). (2) Robust LTS responses never appeared as the amplitude of the simDual-Block waveform is increased (with time constant fixed). This panel provides an explanation for the stark contrast between dynGAT3-Block and dynDual-Block in Figure 2E.

Video 3
Voltage and T-type calcium channel open probability discrepancy trajectories in response to GABAB IPSCs with two different time constants.

Simulated traces for the example neuron of Figure 3, resampled at 1 ms intervals, in response to GABAB IPSCs using the simDual-Block waveform, but with time constants set to 2.0 (yellow) or 2.3 (orange) s. These two time constants correspond to the last LTS success and the first LTS failure, respectively, in Figure 7A2. (A-F) See descriptions for Video 1.

Video 4
Voltage and T-type calcium channel open probability discrepancy trajectories in response to GABAB IPSCs with varying time constants.

Same as Video 3 but in response to all GABAB IPSC time constants shown in Figure 7A2. (A-F) See descriptions for Video 1. There appears to be a threshold for the slope of the open probability discrepancy curve that differentiates between an LTS-producing trajectory and an LTS-lacking trajectory.

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.

Network models

We next explored whether the interplay between GABAB-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 GABAB 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 GABAB receptor-mediated inhibitory synapse (GABAA 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 GABAB receptor activation parameters (simIPSCs), GABAB 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 GABAB-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
GABAB-receptor mediated conductance waveforms modulate oscillations produced by two-cell model thalamic networks.

(A) Schematic of a two-cell model network. A reticular thalamic (RT) neuron projected a GABAB receptor-mediated inhibitory synapse (-) to a thalamocortical (TC) neuron, which reciprocally projected an AMPA receptor-mediated excitatory synapse (+) to the reticular thalamic neuron. (B) Evoked GABAB conductance waveforms in network model TC neurons (solid lines) were similar to GABAB dynIPSC waveforms (dashed lines). (C) Example two-cell network responses under different GABAB receptor conditions. Model TC neuron parameters corresponded to the example neuron of Figure 3. A brief (40 ms, 0.2 nA) current stimulus was applied to the reticular thalamic neuron, evoking an initial burst of 12 spikes. Oscillations were evoked under some but not all GABAB receptor conditions. A total of 24 model TC neurons produced oscillations in response to stimulation. The green dotted line denotes the observed LTS threshold from Figure 5C. (D) Distributions of oscillation measures over all 24, two-cell networks. Oscillation probability was increased when either simGAT1-Block parameters or simGAT3-Block parameters were used, but decreased when simDual-Block parameters were used (*p<0.05, ***p<0.001, Friedman’s test). (E) Proportion of simulations with the TC neuron’s maximum T channel open probability discrepancy between 2 and 3 s after stimulation passing the LTS threshold shown in (C) correlated well with oscillation probability. Each circle denotes a different two-cell network using a specific simIPSC color-coded as in (C). Responses across two-cell networks were heterogeneous.

Figure 8—source data 1

Oscillation measures for two-cell model thalamic networks using different GABAB 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

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−4), 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−2 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 GABAB-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 GABAB 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).

GABAB-receptor-mediated conductance waveforms modulate oscillations produced by 200 cell model thalamic networks.

(A) Schematic of a 200 cell model network. Each reticular thalamic (RT) neuron projected GABAB receptor-mediated inhibitory synapses (-) to nine nearby thalamocortical (TC) neurons. Each TC neuron projected AMPA-receptor-mediated synapses (+) to five nearby reticular thalamic neurons. (B) Sample spike raster plots of a TC-homogenous network, using model TC neuron parameters corresponding to the example neuron of Figure 3. A brief (40 ms, 0.2 nA) current stimulus was applied to each of the center 20 reticular thalamic neurons. Spikes within the stimulation period are red; all other evoked spikes are black. (C) Sample spike raster plots of a TC-heterogeneous network, using model TC neuron parameters corresponding to the 24 model TC neurons used in Figure 8. Relative to TC-homogeneous networks, activity was more localized for TC-heterogeneous networks. (D) Distributions of oscillation measures over all 24 TC-homogeneous 200 cell networks (*p<0.05, **p<0.01, ***p<0.001, repeated-measures ANOVA for oscillation period and half activation time, Friedman’s test otherwise). (E) Distributions of oscillation measures over 24 randomly ordered, TC-heterogeneous 200 cell networks. Oscillation probability, oscillation period and percent of active cells increased when simGAT3-Block parameters were used, but decreased when simDual-Block parameters were used (repeated-measures ANOVA for oscillatory index, Friedman’s test otherwise).

Figure 9—source data 1

Oscillation measures for 200 cell model thalamic networks using different GABAB 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

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−5), 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−7) or simDual-Block parameters (+50%, p=3.8×10−9). 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 GABAB-receptor inhibition. The same interplay between GABAB-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.

Discussion

Here, we show that seizure-like thalamic oscillations were prolonged following individual GAT1 or GAT3 blockade, yet abolished following dual GAT1+GAT3 blockade. We have also shown that relative to control GABAB IPSC waveform responses, thalamocortical neuron rebound burst probability increased following individual GAT1 or GAT3 blockade waveforms, but decreased following dual GAT1+GAT3 blockade waveforms. These rebound burst effects were recapitulated in a large population of computational neurons modeled after biological thalamocortical neurons. When constructed into model thalamic networks, the computational neurons also produced oscillations that responded to different GAT blockade conditions in a manner similar to those recorded in situ. We found that including heterogeneity across thalamocortical neurons in the network significantly increased the robustness of evoked oscillations and accentuated the observed differences across GAT blockade conditions. Finally, we have identified and characterized a link between GABAB-mediated inhibition and T channel gating across both voltage and time dimensions. Specifically, we found that the discrepancy between the instantaneous and the steady-state T channel open probability follows one of two trajectories: (1) failing to reach a threshold, resulting in LTS failure, or (2) driving past the threshold to produce a rebound burst. We discovered that inhibition during single GAT1 or GAT3 blockade waveforms could follow trajectory (2) as the amplitude is increased, but inhibition during the dual GAT1+GAT3 blockade waveform always followed trajectory (1) regardless of amplitude. These observations provide an explanation for the bidirectional effects of GAT blockade on both thalamocortical rebound bursting and seizure-like oscillations.

Role of thalamocortical neuron rebound bursting in generalized seizures

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 GABAB dynIPSC waveform) or decreased (Dual-Block GABAB 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).

Role of the inhibitory temporal envelope in seizures

Prior experimental and computational work has suggested that a shift from GABAA receptor-mediated to GABAB 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 (GABAA: <100 ms; GABAB: 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 GABAB receptor-mediated inhibition. Specifically, systemic injection of GABAB receptor agonists increases SWDs (Liu et al., 1992; Bortolato et al., 2010), while injection of GABAB 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 GABAB IPSCs during SWDs in vivo (Charpier et al., 1999) or find significant SWD changes with GABAB receptor modulation (Staak and Pape, 2001), and have highlighted a particular role for extrasynaptic GABAA receptors that produce a non-dynamic, tonic current in thalamocortical neurons (Cope et al., 2009). Presumably, tonic GABAA 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 GABAA receptor- versus GABAB-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 GABAB 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 GABAB 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 GABAB 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 hT and h,T are low) and a deficiency in T current production (no spike in IT), 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).

T channel open probability discrepancy drives thalamocortical neuron rebound bursting

Our model neurons allowed us to identify a novel measure of T channel gating that provides an explanation for how GAT-modulated, GABAB IPSCs bidirectionally regulates LTS production. We found that when the discrepancy between high instantaneous T channel open probabilities (mT2hT) and low steady-state open probabilities (mT,2hT,) passed a threshold of 10−2, 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 (hT 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 (hT 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 GABAB 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, GABAB 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.

Modeling thalamocortical neuron heterogeneity

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 (IT, Ih, IA, IKir, INaP, INaLeak, IKLeak) 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 (gpas) with a reversal potential (Epas)]. 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 GABAB 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).

Implication on anti-epileptic therapies

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 GABAB receptors (Beenhakker and Huguenard, 2010) in thalamocortical neurons, also exhibit increased incidence of SWDs (Cope et al., 2009). Furthermore, application of GABAB 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 GABAB 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).

Materials and methods

Oscillation recordings in acute thalamic slices

Request a detailed protocol

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 NaHCO3, 1.25 NaH2PO4, 0.5 CaCl2, 10 MgSO4, 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 NaHCO3, 1.25 NaH2PO4, 2 CaCl2, and 1 MgSO4. All solutions were equilibrated with 95% O2/5% CO2.

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 GABAA 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.

Dynamic clamp recordings

Request a detailed protocol

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 NaHCO3, 1.25 NaH2PO4, 0.5 CaCl2, and 10 MgSO4, 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 NaHCO3, 1.25 NaH2PO4, 2 CaCl2, and 1 MgSO4. All solutions were equilibrated with 95% O2/5% CO2.

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 MgCl2, 2 Na2-ATP, 0.5 Na-GTP, 0.07 CaCl2 (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 GABAB-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 GABAB-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 GABAA antagonists to isolate GABAB-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 GABAB-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:

(1) gGABAB=A(1-e-t/τrise)8(we-t/τfallfast+(1-w)e-t/τfallslow)

where A is the amplitude coefficient, τrise is the rise time constant, τfallfast and τ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 GABAB receptor activation.

Table 2
Parameter values of GABAB-mediated inhibitory conductance waveforms used in dynamic clamp experiments.

These values (with the amplitudes scaled by 200%) correspond to the conductance templates shown in Figure 2B. A, amplitude coefficient; τrise, rise time constant; τfallfast, fast decay time constant; τfallslow, slow decay time constant; w, weighting factor for the two decay terms.

Aτriseτfallfastτfallsloww
dynControl16.0052.0090.101073.200.952
dynGAT1-Block24.0052.0090.101073.200.952
dynGAT3-Block8.8838.63273.401022.000.775
dynDual-Block6.3239.8865.802600.000.629

Since the number of GABAB receptors activated in a physiological setting was not determined, all GABAB 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.

Analysis of oscillation recordings

Request a detailed protocol

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 αi was computed from the root-mean-square average of the slope vector over the baseline region before stimulation start, and the maximum slope β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 α=iαi is the baseline slope noise averaged over all sweeps and β=iβi is the maximum slope averaged over all sweeps, then the slice-dependent signal-to-noise ratio was given by r=1+0.1(βα-1). The resulting signal-to-noise ratios fell between 2 and 5. A sweep-dependent slope threshold was then defined by θi=rα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 δ that minimizes the distance function f(δ)=j|pj-qj(δ)|, where pj is the lag value of peak j and qj(δ) is the closest multiple of δ to pj. The search was initialized with the estimate δ0=p2-p1 and confined to the bounds 23δ0,32δ0.

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=(Apeak2-Atrough)/(Apeak1-Atrough), where Apeak1 is the fACF value of the primary peak, Apeak2 is the fACF value of the secondary peak and Atrough 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.

Analysis of dynamic clamp recordings

Request a detailed protocol

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/s2) 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 GABAB IPSC waveform.

Single thalamocortical neuron models

Request a detailed protocol

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 diamsoma. The two dendritic compartments were set to have equal diameters (diamdend), with each having lengths equal to half of the parameter Ldend. 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 (gpas) 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 (IT) and the submembranal calcium extrusion mechanism (Cadecay) was adapted from Destexhe et al., 1998; the hyperpolarization-activated nonspecific cationic current (Ih), the A-type transient potassium current (IA), the inward-rectifying potassium current (IKir) and the persistent sodium current (INaP) were adapted from Amarillo et al., 2014. The following parameters were allowed to vary across compartments and across neurons: the maximum permeability p-T (in cm/s) of IT, which was described by the Goldman–Hodgkin–Katz flux equation (Huguenard and McCormick, 1992), and the maximum conductance densities g-h, g-A, g-Kir, g-NaP (in S/cm2) 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 IT was described by mV,t2hV,t (Huguenard and Prince, 1992). At every point in time, the activation variable m converges exponentially to its steady-state value m(V) with a time constant τm(V) and the inactivation variable h converges exponentially to its steady-state value h(V) with a time constant τh(V). The equations for m(V), τm(V), h(V) and τh(V) 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 VTraub was changed to −65 mV.

A custom GABAB 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 Q10 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.

Parameter initialization

Request a detailed protocol

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:

V(t)=[C0(1et/τ0)+C1(1et/τ1)]B(ttp)+[C0(1etp/τ0)e(ttp)/τ0+C1(1etp/τ1)e(ttp)/τ1]B(t>tp)

where C0, C1 are the amplitudes (in mV) of the two components, τ0, τ1 are the time constants (in ms) of the two components, tp=10 ms is the width of the current pulse, and B(x) is a Boolean function defined by Bx=1 if x=true and Bx=0 if x=false. Initial values for the curve fit were C0,i=V, C1,i=V, τ0,i=10 ms, τ1,i=1 ms 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

RInput=(C0+C1)/Ip,

where Ip=50 pA is the amplitude of the current pulse. The membrane time constant τm was set to equal that of the slow component τm=τ0. The length constant was computed by solving for L in equation 4.5.57 of Johnston and Wu:

|C1/((2C0τ1/τ0)C1)|=cot(α1L)[cot(α1L)1/α1L],

where α1=τ1/τ0-11/2, starting with the initial guess of Li=π/α1. The dendritic-to-somatic conductance ratio ρ was then computed by equation 4.5.58 of Johnston and Wu:

ρ=α1cot(α1L)/coth(L).

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

RMemb=RInputRs,RSoma=(1+ρ)RMemb,RDend=RSoma/ρ,

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

Rm=τm/Cm,

where a fixed value of 0.88 µF/cm2 (Destexhe et al., 1998) for the specific membrane capacitance Cm. Then the somatic diameter diamSoma was computed by equation 4.3.8 of Johnston and Wu:

diamSoma=2Rm/4πRSoma,

the dendritic diameter diamDend was computed by equation 4.5.47 of Johnston and Wu:

diamDend=(2RmRa(cothL)/πRDend)2/3,

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

LDend=LdiamDendRm/4Ra,

and the passive conductance gpas was estimated by

gpas=1/RInputA,

where A=πdiamSoma2+diamDendLDend 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: p-T,Soma=p-T,Dend1=p-T,Dend2=2.0 x 10-4 cm/s, g-h,Soma=g-h,Dend1=g-h,Dend2=2.2 x 10-5 S/cm2, g-Kir,Soma=g-Kir,Dend1=g-Kir,Dend2=2.0 x 10-5 S/cm2, g-A,Soma=g-A,Dend1=g-A,Dend2=5.5 x 10-3 S/cm2 and g-NaP,Soma=g-NaP,Dend1=g-NaP,Dend2=5.5 x 10-6 S/cm2.

Parameter optimization

Request a detailed protocol

The voltage responses of a model neuron to a set of 12 different GABAB 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 Etotal was computed by a weighted average of component errors (Figure 3E):

Etotal=(wmEm+wswEsw+waEa+wtEt+wslEsl)/iwi,

where Em 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), Esw is the average root-mean-square error of voltage values, Ea is the average LTS peak value error across all LTS-matching traces, Et is the average LTS latency error across all LTS-matching traces, Esl is the average LTS maximum slope error across all LTS-matching traces, and the wi’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 Etotal. 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 ρ=1, the coefficient of expansion was χ=2, the coefficient of contraction was γ=0.75, the coefficient of shrinkage was σ=0.8. The initial simplex was set up so that each vertex deviates in a parameter to be varied by δ=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 (-, ) 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 GABAB IPSC waveform was used during model optimization, but all of 60–180 traces were used for the final model evaluation (Figure 3E).

Network models

Request a detailed protocol

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 GABAB 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.

GABAB 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/cm2 (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.

Statistics

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.

Drugs

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).

Code

Request a detailed protocol

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).

Data availability

Source data files have been provided for Figures 1-5, 8-9. Oscillations data, dynamic clamp data are available via Dryad (https://doi.org/10.5061/dryad.4xgxd256f). All code for reproducing results are available online at https://github.com/luadam4c/m3ha_published/ (copy archived at https://github.com/elifesciences-publications/m3ha_published).

The following data sets were generated
    1. Beenhakker MP
    2. Lu AC
    3. Lee CK
    4. Kleiman-Weiner M
    5. Truong B
    6. Wang M
    7. Huguenard JR
    (2020) Dryad Digital Repository
    Nonlinearities between inhibition and T-type calcium channel activity bidirectionally regulate thalamic oscillations.
    https://doi.org/10.5061/dryad.4xgxd256f

References

    1. Borden LA
    2. Dhar TG
    3. Smith KE
    4. Branchek TA
    5. Gluchowski C
    6. Weinshank RL
    (1994)
    Cloning of the human homologue of the GABA transporter GAT-3 and identification of a novel inhibitor with selectivity for this site
    Receptors & Channels 2:207–213.
  1. Book
    1. Johnston D
    2. Wu SM
    (1994)
    Foundations of Cellular Neurophysiology
    MIT press.

Decision letter

  1. Ronald L Calabrese
    Senior Editor; Emory University, United States
  2. Frances K Skinner
    Reviewing Editor; Krembil Research Institute, University Health Network, Canada
  3. Frances K Skinner
    Reviewer; Krembil Research Institute, University Health Network, Canada
  4. Maxim Bazhenov
    Reviewer; University of California, San Diego, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This study effectively combines experimental and modeling work to show and explain the synergistic effects of blocking GABA transporters on thalamocortical oscillations. This is shown to be due to a dependence of T-type calcium channel gating on GABAB receptor activity.

Decision letter after peer review:

Thank you for submitting your article "Nonlinearities between inhibition and T-type calcium channel activity bidirectionally regulate thalamic oscillations" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Frances K Skinner as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Ronald Calabrese as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Maxim Bazhenov (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

As the editors have judged that your manuscript is of interest, but as described below that additional information is required before it is published, we would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). First, because many researchers have temporarily lost access to the labs, we will give authors as much time as they need to submit revised manuscripts. We are also offering, if you choose, to post the manuscript to bioRxiv (if it is not already there) along with this decision letter and a formal designation that the manuscript is "in revision at eLife". Please let us know if you would like to pursue this option. (If your work is more suitable for medRxiv, you will need to post the preprint yourself, as the mechanisms for us to do so are still in development.)

Summary:

This paper addresses the ability of thalamic relay neurons to generate low threshold spikes (and bursts) under conditions when GABA-A-receptors are blocked, and a GABA-B type conductance is applied. The GABA-B conductance has amplitude and kinetics determined by the presence or absence, alone or together, of 2 types of GABA transporter (GAT1 and GAT3). The authors examine effects of GAT1, GAT3 and combined GAT1+GAT3 blocking on thalamocortical bursting activity using combination of experimental techniques and computational modeling (single cell and network). They report when GAT1 or GAT3 is blocked, there was an increase in the thalamocortical oscillations and bursting, however when both GAT 1 and 3 are blocked, the oscillations and bursting were abolished. Based on dynamic-clamp techniques and fitting to the model, the authors identified that interaction of GABA-B with T-channel activation is critical for the differential effect. Specifically, the fast repolarization from a hyperpolarization state led to LTS and eventually to oscillations in case of GAT1/3 blocking. In contrast, blocking both GAT1 and 3 led to prolonged inhibition that prevented LTS. These finding suggests a careful balance between inhibition (GABA-A and B) and calcium currents determines the LTS properties, and bursting in the thalamocortical neurons. Overall, solid work that combines experiment and modeling.

Essential revisions:

While all of the reviewers thought that this was very good work, they all struggled with presentation clarity and a clear summary message about mechanisms involved. The authors are asked to distill and improve the presentation of their work which is very heavy in places.

It was thought that if authors have a clear understanding of principles and mechanisms, they should be able to express it in a clear way in the Abstract/Introduction/Discussion.

For some of the points raised by the reviewers it was decided that it would be at the authors' discretion whether to include them as part of the present revisions (referred to as optional below). The authors could consider them in subsequent work.

1) Experiments: There seem to be washout effect differences with the 2 different blockers. This would seem to potentially play a role on the synergistic abolition? Did you examine oscillations for longer with GAT1 or GAT3 blockers separately as done for the combination? Also, based on the raw recordings shown in Figure 1(iv), the “elimination” could potentially be lower amplitude oscillations? I wondered whether spectral/frequency analyses would be helpful and/or show something additional to PSTHs.

2) Single cell experiments: The dynamic clamp GABAB conductance “approach” would circumvent the experimental variability referred to, it does introduce a different issue in considering the biology – i.e., spatial aspects since one is only looking somatically. So, the comparison can be more “robust” but then the single neuron effect should perhaps be tempered? The authors should comment on how they expect that this would affect their results and interpretation. That is, the GAT3 blockage has an extended temporal effect which is further extended with the combination. I think this should be unpacked a bit more for interpretation purposes.

3) Once one sees Figure 2B, the rest of the aspects of the paper are pretty straightforward and then is mainly about heterogeneous responses moving forward. I am not sure if it is fair to so clearly link the thalamic oscillations directly to the rebound burst aspect (single cell) for these “bidirectional differences”. It then started to get a bit confusing about interpretation of it all (see above point). That is, I think that the authors are assuming an underlying mechanism about TC cells rebound and thalamic oscillations – presumably it is based on previous network modeling studies that they build on later – but the heterogeneous responses and 1 and 2 points above give me a bit of pause. At this stage, they have shown a correlation of network oscillations and single cell responses with GAT blockers.

4) Model: For the examples shown in Figure 3D, could the authors say which cases they are in 3E – presumably well-fit ones? I think that it would be helpful to show the non well-fit cases to see the difference (supplementary figures as in eLife structure could help with this presentation).

5) While the cellular heterogeneity modeling and importance is interesting, and the authors' use of 3 compartment models seems reasonable, the general statements about commonalities of T channel/A channel densities need to be tempered and/or explained a bit more, especially given the limited nature of the models. Figure 3F are where the various parameters are shown but hard to appreciate/understand them in light of the statements. A statement about negative correlation is mentioned as an example. The authors go on to fully explore T-type calcium channels, and bring about the temporal aspects which make sense given kinetic differences noted from the experimental work.

However, I did wonder what one would get if a thorough exploration was done with h-channels to understand the differential effects? Or was this clear already? (this latter point is optional to address in the revised manuscript).

6) The network oscillation modeling and heterogeneous aspects claimed seemed to leave out a lot without mention (see points 1 and 2 above for example). I was also not sure why and whether it was “fair” to add and consider the trial-to-trial variability via the leak reversal potential (given the different TC models)? For example, did it matter which TC model one chose (Figure 3E) to use for the homogeneous networks to compare with heterogeneous? Was some aspect from the 2-cell network used as guidance in the choice for the larger networks? (I may have missed that?)

7) Figure 1C: Why do oscillations seem much shorter "before" drug application in NO and SNAP conditions compare to Control? In fact, in NO condition oscillations seem to be lasting as long as in Control.

8) In both slice experiments and the model, the GABA-A synapse was blocked by bicuculline. However, in-vivo condition would involve the interaction of both GABA-A and GABA-B channels. While this may be difficult to examine in the slice work, it is feasible to examine in the computational model. It would be important to understand how the GABA-A fast acting time scale impacts the effect of slower changing GABA-B synapse, especially in network simulations.

9) Based on the fitting of dynamic-clamp data to the model, the authors' predict high density of T-channel in dendrites and A-channel density in soma. While the rest of the work, especially the computational modeling, has focused on T-channel, the changes to A-channel have not been examined. Similar to the findings in Figure 5-6 on T-channel, it would be helpful to examine the effect of A-channel (this last point is optional to addressed in a revised manuscript).

10) Given that the T-channel is involved in the generation of spindle activity during sleep, it would be helpful to include some discussion about the impact of the findings on sleep spindles.

11) While the authors have developed an interesting measure called the T-channel open channel discrepancy, I feel that the effect could be well captured using conventional phase-plane or phase-space projections. It seems that the membrane voltage, T-channel activation and inactivation variables pass through a critical region in phase-space which is required for the LTS; the failure of this dynamics leads to the lack of LTS. Possibly 3D phase-space analysis of these variables could identify the critical region?

(this point is optional to address in a revised manuscript).

12) The section "Interplay between GABAB receptors and T-type calcium channels " is extremely long and difficult to follow. I am sure some readers will follow all the details there but for many others the message may be lost. I suggest some rewriting to include summaries and highlight the main results.

13) In the network simulations (Figure 8), for TC-homogeneous case, oscillations seem to be more synchronized in control vs GAT1 or GAT3 block conditions. Why is this so?

14) The way the manuscript is written, it is somewhat difficult to learn what is actually a precise mechanism behind the differences between effects of individual GAT1/3 blocking vs dual blocking. E.g., the Summary in discussion says: "…respectively, can be recapitulated by varying the GABAB receptor activation waveform." but it does not provide a needed summary of the mechanisms discovered in the study.

15) Consequences of the GABA-B conductance properties are reflected in the behavior of individual cells, and of oscillating networks of TCR and nRT (nucleus reticularis) neurons. Details of the GABA-B conductances to be used were determined in prior publications. Experimental and numerical experiments were performed in thalamocortical slices, in dynamic clamp applied to single TCR neurons, and in simulations of 3-compartment TCR neurons alone, paired with an nRT cell, and in a TCR/nRT network. The main finding is that excessively large and prolonged GABA-B conductances interfere with LTS generation; detailed analysis shows that the specific kinetic properties of T-channels are responsible, in that LTS requires an epoch where inward current is regenerative (increasing m) at the same time as inactivation is not too developed, and the latter in turn depends on instantaneous inactivation (h) being larger than would occur at equilibrium. Such large slow conductances develop when both transport blockers are applied together, but not separately. The authors provide some discussion of possible relations to spike-wave seizures, and the effects and contrary side-effects of certain anticonvulsants.

I suspect that the mathematics involved is more general than discussed here for T-type calcium channels, and depends at heart (no pun intended) on the slowness of inactivation vs activation. An example might be the inability of a squid axon to fire at low frequencies, even in the presence of a low calcium buffer. The authors may, at their discretion, want to consider this larger issue (optional to address in a revised manuscript).

https://doi.org/10.7554/eLife.59548.sa1

Author response

Essential revisions:

While all of the reviewers thought that this was very good work, they all struggled with presentation clarity and a clear summary message about mechanisms involved. The authors are asked to distill and improve the presentation of their work which is very heavy in places.

It was thought that if authors have a clear understanding of principles and mechanisms, they should be able to express it in a clear way in the Abstract/Introduction/Discussion.

For some of the points raised by the reviewers it was decided that it would be at the authors' discretion whether to include them as part of the present revisions (referred to as optional below). The authors could consider them in subsequent work.

1) Experiments: There seem to be washout effect differences with the 2 different blockers. This would seem to potentially play a role on the synergistic abolition? Did you examine oscillations for longer with GAT1 or GAT3 blockers separately as done for the combination? Also, based on the raw recordings shown in Figure 1(iv), the “elimination” could potentially be lower amplitude oscillations? I wondered whether spectral/frequency analyses would be helpful and/or show something additional to PSTHs.

The effects of washing out the drugs were not that different for NO-711 vs SNAP-5114 (during the washout of either drug, the duration of evoked oscillations decreased), although the washout was generally faster for NO-711. To prevent confusion, we have selected a more representative example for SNAP-5114. In Beenhakker and Huguenard, 2010, GABAB IPSCs reached a stable plateau level for the duration of drug application within 10 minutes of applying either GAT1, GAT3, or dual GAT1+GAT3 blockers. Similarly, we observed in this study that oscillation durations stabilized within 15 minutes of applying either GAT1 or GAT3 blockers. We thus did not examine oscillations for longer than 40 minutes of either NO-711 or SNAP-5114 application because oscillations had stabilized and the prolongation in duration was robust (Figure 1D). Therefore, if wash-in periods were extended, we would not anticipate an eventual abolishment of oscillations following either single GAT1 or GAT3 blockade, unlike what we observed during dual GAT1+GAT3 blockade (Figure 1F).

We interpret “remnant low amplitude oscillations” as meaning the possibility of oscillatory spiking activity that is buried in the noise. The reviewer’s comment is well taken and we have now performed the requested spectral/frequency analysis for the recordings (Author response image 1). As we expected, the analysis failed to resolve any spiking activity in the recording of the dual block condition, while activity in the other conditions is readily apparent. Therefore, it seems quite unlikely that activity during the dual blockade condition is robust, but unresolvable with our spike detection algorithm. The original raw recordings in Figure 1B(iv) were selected for a slice where complete oscillation abolishment took longer than 40 minutes. To avoid further confusion, we have now selected a more representative example of a raw recording during dual blockade where complete oscillation abolishment occurred.

Author response image 1
Individual GAT1 or GAT3 blockade strengthens thalamic oscillations, but dual GAT1+GAT3 blockade abolishes oscillations.

Spectrograms for example evoked epileptiform oscillations at baseline and 40 minutes after perfusing with (i) control (no drug added), (ii) 4 µM NO-711 (GAT1 blocker), (iii) 100 µM SNAP-5114 (GAT3 blocker) or (iv) dual 4 µM NO-711+100 µM SNAP-5114. These are the same example oscillations as in Figure 1B. A bin width of 100 ms was used. Importantly, these analyses do not rely on any spike detection algorithms.

2) Single cell experiments: The dynamic clamp GABAB conductance “approach” would circumvent the experimental variability referred to, it does introduce a different issue in considering the biology – i.e., spatial aspects since one is only looking somatically. So, the comparison can be more “robust” but then the single neuron effect should perhaps be tempered? The authors should comment on how they expect that this would affect their results and interpretation. That is, the GAT3 blockage has an extended temporal effect which is further extended with the combination. I think this should be unpacked a bit more for interpretation purposes.

Since the GABAB receptor-mediated IPSCs in Beenhakker and Huguenard, 2010, were recorded from the soma, it is natural to introduce the corresponding conductance waveforms in the soma through dynamic clamp. Although this approach may lead to a delayed voltage response in the dendrites, the fact that voltage attenuation from soma to dendrite is small should alleviate this concern. In fact, our analysis of thalamocortical neuron current pulse responses yielded an average electrotonic length of 0.71 ± 0.08 space constants (range: 0.1-1.5 space constants). This means that voltage changes at the distal dendrite would be on average at least 50% of that in the soma at any given point in time. Furthermore, we saw the same rebound bursting differences across dynIPSCs as the amplitude was scaled between a range of 50-800% in our dynamic clamp experiments (Figure 2E). These amplitude scales should address most voltage attenuation concerns. Nonetheless, we have included a cautionary statement in the text.

3) Once one sees Figure 2B, the rest of the aspects of the paper are pretty straightforward and then is mainly about heterogeneous responses moving forward. I am not sure if it is fair to so clearly link the thalamic oscillations directly to the rebound burst aspect (single cell) for these “bidirectional differences”. It then started to get a bit confusing about interpretation of it all (see above point). That is, I think that the authors are assuming an underlying mechanism about TC cells rebound and thalamic oscillations – presumably it is based on previous network modeling studies that they build on later – but the heterogeneous responses and 1 and 2 points above give me a bit of pause. At this stage, they have shown a correlation of network oscillations and single cell responses with GAT blockers.

We agree that we may have overstated the link between T channel gating and oscillations. The phrase "in agreement" has been changed to “correlated”. We have also cited studies that linked thalamocortical neuron rebound bursting with thalamic oscillations in the beginning of the "Single neuron recordings" section. Throughout the text, we have now clarified that the link between T channel gating and oscillations is at this point only correlational. We have also added a figure panel that shows a high correlation between our measure that predicts the production of a rebound low threshold spike in thalamocortical neurons (“open probability discrepancy”, Figure 5C) with the probability for oscillations to appear in a network (see Figure 8E).

4) Model: For the examples shown in Figure 3D, could the authors say which cases they are in 3E – presumably well-fit ones? I think that it would be helpful to show the non well-fit cases to see the difference (supplementary figures as in eLife structure could help with this presentation).

The ranks are now stated in the Figure 3 legend and labelled in Figure 3—figure supplement 1. We also provided fits for 4 additional well-fitted cases and the 2 excluded cases in the supplement.

5) While the cellular heterogeneity modeling and importance is interesting, and the authors' use of 3 compartment models seems reasonable, the general statements about commonalities of T channel/A channel densities need to be tempered and/or explained a bit more, especially given the limited nature of the models. Figure 3F are where the various parameters are shown but hard to appreciate/understand them in light of the statements. A statement about negative correlation is mentioned as an example. The authors go on to fully explore T-type calcium channels, and bring about the temporal aspects which make sense given kinetic differences noted from the experimental work.

However, I did wonder what one would get if a thorough exploration was done with h-channels to understand the differential effects? Or was this clear already? (this latter point is optional to address in the revised manuscript).

We included general statements about commonalities of T channel/A channel densities because these findings converged on previous results in the literature. The relevant literature has now been cited where appropriate. Nonetheless, in order to focus the manuscript on the interplay between inhibition and T channels, we have now moved descriptions of other channels to Figure 5—figure supplement 1. We agree that a thorough exploration of h channels would be interesting, but it is probably beyond the scope of this work.

6) The network oscillation modeling and heterogeneous aspects claimed seemed to leave out a lot without mention (see points 1 and 2 above for example). I was also not sure why and whether it was “fair” to add and consider the trial-to-trial variability via the leak reversal potential (given the different TC models)? For example, did it matter which TC model one chose (Figure 3E) to use for the homogeneous networks to compare with heterogeneous? Was some aspect from the 2-cell network used as guidance in the choice for the larger networks? (I may have missed that?)

We apologize for the confusion. Leak reversal potential was not allowed to vary during the fitting process. All TC model neurons had a leak reversal potential of -70 mV. Therefore, we feel that it is fair to compare simulations derived from all 2-cell networks subjected to the same range of leak reversal potentials. We ran simulations using 14 different leak reversal potentials, and ran 5 different trials per leak reversal potential. The 5 trials differed by introducing a different leak conductance, and this leak conductance was randomized about a mean leak conductance that was based on the TC model for each network. It is customary to perform such randomization across trials (Sohal and Huguenard, 2003).

Each 2-cell network indeed behaved differently. To account for TC neuron heterogeneity, we examined all possible 2-cell networks using each of the well-fitted model TC neurons (Figure 8D). Similarly, we examined all possible TC-homogeneous networks using each of the well-fitted model TC neurons (Figure 9D). Due to the response heterogeneity, we also examined TC-heterogeneous networks that included all of the well-fitted model TC neurons, randomly ordered (Figure 9E). We now include a description of these analyses:

“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.”

7) Figure 1C: Why do oscillations seem much shorter "before" drug application in NO and SNAP conditions compare to Control? In fact, in NO condition oscillations seem to be lasting as long as in Control.

Oscillation durations are highly variable in acute slice experiments as the extent of the intact network in each slice is difficult to control. Examples shown in Figure 1C have been reselected to show similar baseline durations for Control, NO-711 and SNAP-5114. Furthermore, the distribution of baseline durations over all slices in Figure 1D were similar across all four pharmacological conditions (Control: 6.9 ± 2.4 seconds, NO-711: 6.6 ± 2.5 seconds, SNAP-5114: 4.9 ± 1.8 seconds, NO-711+ SNAP-5114: 7.7 ± 1.7 seconds).

8) In both slice experiments and the model, the GABA-A synapse was blocked by bicuculline. However, in-vivo condition would involve the interaction of both GABA-A and GABA-B channels. While this may be difficult to examine in the slice work, it is feasible to examine in the computational model. It would be important to understand how the GABA-A fast acting time scale impacts the effect of slower changing GABA-B synapse, especially in network simulations.

We have now performed both 2-cell network simulations (Author response image 2) and 200-cell TC-heterogeneous network simulations (Author response image 3) with GABAA receptors added to model TC neurons.

Observations on oscillations generated by networks (2 and 200) that incorporate GABAA receptors include: (1) a decreased oscillation period across most conditions, as one would expect from the briefer GABAA kinetics, and (2) a decreased oscillation probability when simControl, simGAT1-Block or simGAT3-Block waveform was used. Unlike networks containing only GABAB receptors (and unlike our experimental observations), oscillation probability increased when the simDual-Block waveform was used. We attribute this effect to the briefer decay kinetics provided by GABAA receptor mediated inhibition, a condition that promotes rebound bursting by increasing T channel open probability discrepancy (as demonstrated in Figure 6D).

Nevertheless, the oscillation differences across GAT blockade conditions derived from network simulations including GABAA receptors were the same as those derived from network simulations without GABAA receptors (cf. Figure 9E and Author response image 3C). That is, relative to using the simControl waveform, we observed an increase in oscillation probability using the simGAT1-Block or simGAT3-Block waveform and a decrease in oscillation probability using the simDual-Block waveform. Relative to using the simControl waveform, we also observed an increase in oscillation period using the simGAT3-Block or simDual-Block waveform.

However, we have not included these results in the manuscript as we do not have either voltage clamp recordings to physiologically constrain GABAA receptor parameters, or oscillation recordings to support the observations. In the simulations above, we simply applied GABAA receptor conductance kinetics from Huntsman et al., 1999, and set an amplitude ratio of 2 (relative to the GABAB conductance) estimated from Huguenard and Prince, 1994. However, these may not reflect physiological changes in GABAA receptor mediated-IPSCs under various forms of GAT blockade. Furthermore, in order to support the observations above, one would need to record oscillations in which GABAA receptors are present in TC neurons but not in RT neurons, which is not easy to do pharmacologically. One possibility is to either record oscillations in β3 subunit knockout mice (Huntsman et al., 1999) or a comparable genetically-modified rat, but we feel that is beyond the scope of this work. If the reviewers feel strongly about including these simulation findings, then we can add them as a supplement to Figure 9.

Author response image 2
With GABAA receptors present, GABAB-receptor mediated conductance waveforms still modulate oscillations produced by 2-cell model thalamic networks.

(A) Schematic of a 2-cell model network with GABAA receptors present (cf. Figure 8A). (B) Example 2-cell network responses under different GABAB receptor conditions. The same TC model neuron and simulation conditions as that in Figure 8C were applied, except for the addition of a GABAA conductance with maximal amplitude 2 times that of the maximal GABAB conductance. The green dotted line denotes the observed LTS threshold from Figure 5C. (C) Distributions of oscillation measures over the same 24, 2-cell networks as in Figure 8D. Compared to Figure 8D, 2-cell network simulations with GABAA receptors present had a lower oscillation probability and a smaller oscillation period (* p < 0.05, *** p < 0.001, Friedman’s test). The differences across simIPSCs were generally the same as in Figure 8D. (D) Figure 8D reproduced here for ease of comparison.

Author response image 3
GABAB-receptor mediated conductance waveforms modulate oscillations produced by 200-cell model thalamic networks with GABAA receptors present.

(A) Schematic of a 200-cell model network with GABAA receptors present (cf. Figure 9A). (B) Sample spike raster plots of the same TC-heterogeneous network and the same simulation conditions as that in Figure 9C, except for the addition of a GABAA conductance with maximal amplitude 2 times that of the maximal GABAB conductance. (C) Distributions of oscillation measures over the same 24 randomly ordered, TC-heterogeneous 200-cell networks as in Figure 9E (* p < 0.05, ** p < 0.01, *** p < 0.001, Friedman’s test for oscillation period and half activation time, repeated-measures ANOVA otherwise). The differences across simIPSCs were generally the same as in Figure 9E. (C) Figure 9E reproduced here for ease of comparison.

9) Based on the fitting of dynamic-clamp data to the model, the authors' predict high density of T-channel in dendrites and A-channel density in soma. While the rest of the work, especially the computational modeling, has focused on T-channel, the changes to A-channel have not been examined. Similar to the findings in Figure 5-6 on T-channel, it would be helpful to examine the effect of A-channel (this last point is optional to addressed in a revised manuscript).

Prior literature (Pape et al., 1994; Amarillo et al., 2014) has already examined the effects of A channels on low-threshold spikes, demonstrating the importance of the A current for controlling the amplitude and width, but not the initiation, of the LTS. It would be interesting to examine the effects of A channels on network oscillations, but since the effects of T channels are sufficient to explain most of our experimental findings, we feel like the investigation of A channels is beyond the scope of this work.

10) Given that the T-channel is involved in the generation of spindle activity during sleep, it would be helpful to include some discussion about the impact of the findings on sleep spindles.

As sleep spindles and absence seizures appear to rely on similar cellular- and circuit-level mechanisms, this point is reasonable to raise. However, sleep spindles have been shown to be primarily GABAA receptor-dependent, and are likely not as dependent on GABAB receptor-mediated signaling. Since our study focuses on GABAB receptor modulation of thalamic circuits, we are not comfortable translating our findings to sleep spindles.

11) While the authors have developed an interesting measure called the T-channel open channel discrepancy, I feel that the effect could be well captured using conventional phase-plane or phase-space projections. It seems that the membrane voltage, T-channel activation and inactivation variables pass through a critical region in phase-space which is required for the LTS; the failure of this dynamics leads to the lack of LTS. Possibly 3D phase-space analysis of these variables could identify the critical region?

(this point is optional to address in a revised manuscript).

In this work, we have demonstrated that the change of either membrane voltage V, T channel activation variable mT or inactivation variable hT relative to time is critical for LTS production. Notably, these derivatives are independent of each other. We also demonstrated that either dV/dt (slope of voltage) and d[log(mT2hTmT,2hT,)]/dt (slope of open probability discrepancy) have some predictive power towards the production of an LTS. Therefore, we think that as opposed to 3D phase plane analysis with V, mT and hT, a phase plane analysis with dV/dt, d(mT)/dt and d(hT)/dt might provide an opportunity to identify a "critical region required for the LTS”. However, we also expect this critical region to be cell dependent. We acknowledge that such further analysis would be interesting, but we also feel that it would add unnecessary complexity for the general reader.

12) The section "Interplay between GABAB receptors and T-type calcium channels " is extremely long and difficult to follow. I am sure some readers will follow all the details there but for many others the message may be lost. I suggest some rewriting to include summaries and highlight the main results.

We thank the reviewer for bringing this point to our attention. We hope to make this section as clear as possible. Therefore, we have shortened this section by moving less pertinent details to the Figure 5 and Figure 7—figure supplement 1. We have also rewritten several sections to simplify the concepts we address and moved some interpretation to the Discussion. Additionally, we agree that the original presentation of Figure 6 may have been somewhat overwhelming. Therefore, we have now split Figure 6 into two new figures (Figure 6 and Figure 7).

13) In the network simulations (Figure 8), for TC-homogeneous case, oscillations seem to be more synchronized in control vs GAT1 or GAT3 block conditions. Why is this so?

We looked into this curious observation and quantified a metric we call burst latency jitter. This parameter reflects the variability in burst latencies across trials (i.e., jitter = standard deviation). Burst latency jitter was significantly higher for both simGAT1-Block and simGAT3-Block across the 31 model neurons (Author response image 4). With higher burst latency jitter, oscillations become more desynchronized.

Author response image 4
Burst latencies were more variable following simIPSCs corresponding to single GAT blockade across well-fitted model neurons.

(A) Distributions of burst latency jitter (standard deviation of burst latencies across all trials) over the 31 well-fitted model neurons across GABAB IPSC waveforms shown in Figure 2B (* p < 0.05, ** p < 0.01, *** p < 0.001, Friedman’s test). (B) Burst latency jitter averaged over all 31 model neurons, across 4 different GABAB IPSC waveforms and 3 different conductance amplitude scales. Error bars denote 95% confidence intervals.

14) The way the manuscript is written, it is somewhat difficult to learn what is actually a precise mechanism behind the differences between effects of individual GAT1/3 blocking vs dual blocking. E.g., the Summary in discussion says: "…respectively, can be recapitulated by varying the GABAB receptor activation waveform." but it does not provide a needed summary of the mechanisms discovered in the study.

We have rewritten the first paragraph of the Discussion to highlight the most important results of this work, including a summary of the discovered mechanism underlying the differences between single and dual GAT blockade. We have also clarified the section in the results wherein we highlight our mechanism.

15) Consequences of the GABA-B conductance properties are reflected in the behavior of individual cells, and of oscillating networks of TCR and nRT (nucleus reticularis) neurons. Details of the GABA-B conductances to be used were determined in prior publications. Experimental and numerical experiments were performed in thalamocortical slices, in dynamic clamp applied to single TCR neurons, and in simulations of 3-compartment TCR neurons alone, paired with an nRT cell, and in a TCR/nRT network. The main finding is that excessively large and prolonged GABA-B conductances interfere with LTS generation; detailed analysis shows that the specific kinetic properties of T-channels are responsible, in that LTS requires an epoch where inward current is regenerative (increasing m) at the same time as inactivation is not too developed, and the latter in turn depends on instantaneous inactivation (h) being larger than would occur at equilibrium. Such large slow conductances develop when both transport blockers are applied together, but not separately. The authors provide some discussion of possible relations to spike-wave seizures, and the effects and contrary side-effects of certain anticonvulsants.

I suspect that the mathematics involved is more general than discussed here for T-type calcium channels, and depends at heart (no pun intended) on the slowness of inactivation vs activation. An example might be the inability of a squid axon to fire at low frequencies, even in the presence of a low calcium buffer. The authors may, at their discretion, want to consider this larger issue (optional to address in a revised manuscript).

Comparing the dynamics of low-threshold spikes to those of action potentials would be interesting indeed. It is true that the sodium channel’s inactivation time constant is also much higher than its activation time constant, just like the T channel. However, there are additional considerations that make a direct comparison difficult, such as the fact that sodium channels are not inactivated at rest, and thus do not need to be recovered. Therefore, we expect that similar dynamics would be at play in terms of a lack of channel opening in response to slow voltage changes, but the post-inhibitory aspect would be lacking for sodium channels.

References:

Huntsman, M. M., Porcello, D. M., Homanics, G. E., DeLorey, T. M., & Huguenard, J. R. (1999). Reciprocal Inhibitory Connections and Network Synchrony in the Mammalian Thalamus. Science, 283(5401), 541–543. https://doi.org/10.1126/science.283.5401.541

https://doi.org/10.7554/eLife.59548.sa2

Article and author information

Author details

  1. Adam C Lu

    Department of Pharmacology, University of Virginia, Charlottesville, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    al4ng@virginia.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1008-1057
  2. Christine Kyuyoung Lee

    Department of Neurosurgery, Massachusetts General Hospital, Boston, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Investigation, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1422-4606
  3. Max Kleiman-Weiner

    Department of Psychology, Harvard University, Cambridge, United States
    Contribution
    Software, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6067-3659
  4. Brian Truong

    Department of Pharmacology, University of Virginia, Charlottesville, United States
    Contribution
    Software
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0179-0932
  5. Megan Wang

    Princeton Neuroscience Institute, Princeton University, Princeton, United States
    Contribution
    Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8845-4936
  6. John R Huguenard

    Department of Neurology, Stanford University, Palo Alto, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Methodology, Project administration, Writing - review and editing
    For correspondence
    John.Huguenard@stanford.edu
    Competing interests
    Senior editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6950-1191
  7. Mark P Beenhakker

    Department of Pharmacology, University of Virginia, Charlottesville, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Investigation, Methodology, Project administration, Writing - review and editing
    For correspondence
    mpb5y@virginia.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4541-0201

Funding

National Institute of Neurological Disorders and Stroke (NIH grant R01-NS099586)

  • Adam C Lu
  • Brian Truong
  • Mark P Beenhakker

National Institute of Neurological Disorders and Stroke (NIH grant R01-NS034774)

  • Christine Kyuyoung Lee
  • Max Kleiman-Weiner
  • Megan Wang
  • John R Huguenard

University of Virginia (Whitfield-Randolph Scholarship)

  • Adam C Lu

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

Acknowledgements

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.

Ethics

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.

Senior Editor

  1. Ronald L Calabrese, Emory University, United States

Reviewing Editor

  1. Frances K Skinner, Krembil Research Institute, University Health Network, Canada

Reviewers

  1. Frances K Skinner, Krembil Research Institute, University Health Network, Canada
  2. Maxim Bazhenov, University of California, San Diego, United States

Publication history

  1. Received: June 1, 2020
  2. Accepted: September 8, 2020
  3. Accepted Manuscript published: September 9, 2020 (version 1)
  4. Version of Record published: October 1, 2020 (version 2)

Copyright

© 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.

Metrics

  • 980
    Page views
  • 126
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Medicine
    Homa MohammadiPeyhani et al.
    Tools and Resources

    The discovery of a drug requires over a decade of intensive research and financial investments – and still has a high risk of failure. To reduce this burden, we developed the NICEdrug.ch resource, which incorporates 250,000 bioactive molecules, and studied their enzymatic metabolic targets, fate, and toxicity. NICEdrug.ch includes a unique fingerprint that identifies reactive similarities between drug–drug and drug–metabolite pairs. We validated the application, scope, and performance of NICEdrug.ch over similar methods in the field on golden standard datasets describing drugs and metabolites sharing reactivity, drug toxicities, and drug targets. We use NICEdrug.ch to evaluate inhibition and toxicity by the anticancer drug 5-fluorouracil, and suggest avenues to alleviate its side effects. We propose shikimate 3-phosphate for targeting liver-stage malaria with minimal impact on the human host cell. Finally, NICEdrug.ch suggests over 1300 candidate drugs and food molecules to target COVID-19 and explains their inhibitory mechanism for further experimental screening. The NICEdrug.ch database is accessible online to systematically identify the reactivity of small molecules and druggable enzymes with practical applications in lead discovery and drug repurposing.

    1. Computational and Systems Biology
    2. Microbiology and Infectious Disease
    Wellington Miranda S et al.
    Research Article Updated

    Many bacteria communicate with kin and coordinate group behaviors through a form of cell-cell signaling called acyl-homoserine lactone (AHL) quorum sensing (QS). In these systems, a signal synthase produces an AHL to which its paired receptor selectively responds. Selectivity is fundamental to cell signaling. Despite its importance, it has been challenging to determine how this selectivity is achieved and how AHL QS systems evolve and diversify. We hypothesized that we could use covariation within the protein sequences of AHL synthases and receptors to identify selectivity residues. We began by identifying about 6000 unique synthase-receptor pairs. We then used the protein sequences of these pairs to identify covariation patterns and mapped the patterns onto the LasI/R system from Pseudomonas aeruginosa PAO1. The covarying residues in both proteins cluster around the ligand-binding sites. We demonstrate that these residues are involved in system selectivity toward the cognate signal and go on to engineer the Las system to both produce and respond to an alternate AHL signal. We have thus demonstrated that covariation methods provide a powerful approach for investigating selectivity in protein-small molecule interactions and have deepened our understanding of how communication systems evolve and diversify.