A multi-compartment model for interneurons in the dorsal lateral geniculate nucleus

GABAergic interneurons (INs) in the dorsal lateral geniculate nucleus (dLGN) shape the information flow from retina to cortex, presumably by controlling the number of visually evoked spikes in thalamocortical (TC) cells, and refining their receptive field. The dLGN INs exhibit a rich variety of firing patterns. Depolarizing current injections to the soma may induce tonic firing, periodic bursting or an initial burst followed by tonic spiking, sometimes with prominent spike time adaptation. When released from hyperpolarization, some dLGN INs elicit rebound bursts, while others return more passively to the resting potential [1-3]. 
 
A full mechanistic understanding that explains the function of the dLGN on the basis of neuronal morphology, physiology and circuitry is currently lacking. One way to address this question is to develop and investigate mathematical models of the involved cells and their interactions. While detailed models of the well studied TC cells exist, previous models of the less studied dLGN INs are simplified and use only passive dendritic properties. Several ion channels are present in the dendrites of dLGN INs [4-6], and may be of particular importance in these neurons, as their dendrites have not only postsynaptic contacts for excitatory retinal input, but also presynaptic terminals for inhibitory output to TC-cell dendrites [7]. 
 
We here present a detailed compartmental model of a dLGN IN with a detailed morphology and active dendritic, as well as somatic, conductances. The simulation tool NEURON (http://neuron.duke.edu/) was used for the model. The conductance values were constrained by somatic voltage recordings from two dLGN INs under 8 experimental conditions. The model reproduces qualitative features in the experimentally recorded response patterns, as well as quantitative data on the firing frequency as a function of injected current. We show that relative differences in conductance values, rather than differences in ion channel composition, explain the differences in the response properties of the two neurons, and can account for the various response patterns listed above. The model lays the foundation for studying the function of the dLGN IN and its interactions with TC cells within a network context, and does also in its own right give new insights in the properties of dLGN IN activity.


Introduction
The dorsal lateral geniculate nucleus (dLGN) receives input from retinal ganglion cells and transmits processed information to visual cortex. About 75-80% of the neurons in the dLGN are thalamocortical (TC) neurons, also called relay neurons, as they relay information from the retina to the cortex. Local GABAergic interneurons (INs) constitute the remaining 20-25%, and are responsible for most of the intra-nuclear connections [1]. By providing feed-forward inhibition from retinal ganglion cells to TC neurons, INs control the number of visually evoked spikes in TC neurons, and refine the receptive fields of TC neurons (e.g. [2][3][4][5]). The INs are also important for synchronizing thalamic oscillations [6][7].
Comprehensive mathematical network models will likely be important for a comprehensive understanding of the key functional features of early sensory processing [8][9][10][11]. Due to the relatively high abundance of experimental data, the early visual system has attracted particular interest from theoretical neurobiologists. Several mechanistic network models aimed at mimicking responses of neurons in primary visual cortex, and neurons in the dLGN that provide the feed-forward input to visual cortex, have been developed [12][13][14][15][16][17][18][19][20][21]. Such network models require (i) detailed information about network connectivity and (ii) mathematical neuron models that capture the salient physiological properties of the individual neurons types.
For dLGN, a host of physiological and anatomical experiments have provided detailed information about the neuronal connectivity, as well as the morphology and physiology of TC neurons and INs (see e.g. reviews in [22][23][24]). Different electrophysiological characteristic of TC neurons have been captured in a series of modeling works [25][26][27][28][29][30][31], and the accumulated insight has been incorporated in a high-resolution model which comprises a detailed description of the cell morphology, a set of different ion channels and their distribution over the somatodendritic membrane [32]. For the INs, the situation is more problematic: models are few and less detailed [33][34][35][36]. Until now, network models for the early visual system have either omitted INs entirely, or represented them in a very simplified manner (but see [37]). A satisfactory theoretical understanding of the computational properties of the dLGN circuit, and thus also the input to visual cortex, will likely require network models incorporating more detailed IN models. The development of such models is the topic of the present paper.
INs may exhibit a rich variety of firing patterns, including (i) initial ''sags'' by hyperpolarizing current injections, (ii) rebound bursts when released from hyperpolarization, (iii) tonic firing of action potentials (APs) by depolarizing current injections, (iv) initial bursts by depolarizing current injections, (v) spike-time adaptation during depolarizing current injections, and (vi) periodic bursting during depolarizing current injections [36,[38][39][40][41][42][43][44]. Most previous models focus on aspects of passive signal propagation in INs [33][34][35]. To our knowledge, the only currently available model that includes a variety of the active mechanisms in INs was primarily developed in order to study the mechanisms behind the rebound bursts [36], whereas other properties such as dendritic conductances and the relationship between somatic current injections and action potentials frequency (hereby referred to as the I/O curve) were not taken into account.
The morphology and distribution of dendritic ion channels are crucial for the integration of synaptic input (reviewed in [45]), and can even influence the neuron's response to somatic current injections [46][47][48][49]. Several active conductances have been identified in the dendrites of INs [38,[50][51][52][53][54]. Dendritic ion channels are likely of particular importance in INs, as their dendrites have not only postsynaptic contacts for excitatory retinal and cortical input, but also presynaptic terminals for inhibitory output to TC dendrites [2]. An understanding of how INs provides feed forward inhibition to TCs thus requires models that incorporate the electrically active processes in the dendritic tree. No previous IN model includes such properties.
We here propose a new and detailed multi-compartment model of the IN, which advances previous models in several aspects. Firstly, it includes a detailed 3D reconstruction of IN morphology, and a set of somatic and dendritic active conductances which could reproduce some of the key response patterns of INs (including (i)-(vi) listed above). In this way it lays the foundation for simulating active dendritic signaling on a fine spatiotemporal scale. Secondly, we present two different parameterizations (P1 and P2) of the model, which were constrained by current-clamp data from two example neurons (IN1 and IN2). This approach allowed us to do a comparative study, and relate differences in the parameterizations P1 and P2 to differences in the firing patterns between IN1 and IN2. The two parameterizations contain the same set of ion channels such that differences in the responses of the two neurons could be explained solely in terms of relative differences in the peak values of the conductances. Thereby, we demonstrate that relative differences in conductance values of the included ion channels may account for the substantial variations in firing patterns observed between INs.
Finally, we were able to explain the experimentally observed responses in different neurons (IN1 or IN2) under 8 different input conditions, obtaining quantitative agreement with the experimental I/O-curves over the entire input range studied. This is a significant advance compared to previous models of similar cells. Since the overall input to any given cell in a network generally varies with time over a wide dynamic range, our model will likely allow a more accurate representation of the integration and computational properties of these neurons when included in a realistic dLGN network.
Preliminary results from this model have previously been presented in abstract form [55]. The model will be publicly available on ModelDB (http://senselab.med.yale.edu/modeldb).
Whole-cell voltage-or current-clamp recordings were made from INs in dLGN. Neurons were visualized using DIC optics and infrared video microscopy. INs were identified by expression of GFP which was specifically expressed in GABAergic neurons under control of the endogenous GAD67 promoter in the GAD67-GFP knock-in mice [56] we used. Recordings were obtained with borosilicate glass electrodes (4-6 MV) filled with (mM): 115 potassium gluconate, 20 KCl, 10 HEPES, 2 MgCl 2 , 2 MgATP, 2 Na 2 ATP, 0.3 GTP (pH adjusted to 7.3 with KOH). For morphology reconstruction, biocytin (0.25%; Sigma-Aldrich, St Louis, USA) was included in the intracellular solution.
Current traces were recorded and filtered at 3 kHz with a HEKA EPC 9 amplifier (HEKA Electronik, Lambrecht, Germany), while voltage traces were recorded and filtered at 10 kHz with an Axoclamp 2A amplifier (Molecular Devices, Palo Alto, CA, USA).
After recordings, slices were fixed by 0.1 M phosphate buffer containing 4% paraformaldehyde and kept there for at least 24 h.

Author Summary
The dorsal lateral geniculate nucleus (dLGN) is a part of the visual thalamus. This region contains two types of neurons: thalamocortical neurons and local interneurons. Thalamocortical neurons receive information from the retina and transmit information to visual cortex. The interneurons regulate the activity of thalamocortical neurons through inhibitory connections. This regulation is not properly understood, but it is believed to promote contrast enhancement and other vital visual functions. A powerful tool for development of a mechanistic understanding of dLGN functions is computer models that include the involved neurons, their interconnections and their interactions. Quite sophisticated models are available for thalamocortical neurons, but previous interneuron models are too simple for adequate mechanistic understanding of the functional properties of interneurons. We here present a detailed compartmental interneuron-model based on experimental data. The typical response patterns vary between different interneurons, but also within a given neuron, depending on the stimulus it receives. The model identifies a set of ionic mechanisms that can explain this diversity of activity patterns. In addition to being a useful building block for future network simulations of the dLGN, the model gives useful insight into the operating principles of dLGN interneurons.

Computer modeling
Morphology. In this study, we use our own 3D reconstructions of the morphology of mouse INs. The model was based on the selected morphology shown in Figure 1A, which is used in all simulations unless otherwise stated. The total surface area was 9864 mm 2 , the total length of dendrites was 5771 mm, the longest dendrite was 673 mm, and the mean somatodendritic diameter was about 0.5 mm. The model morphology contained 104 sections that were split into a total of 330 segments. Additional test simulations were run using morphologies of (a) a similar (9566 mm 2 ), (b) a smaller (7071 mm 2 ), and (c) a larger (14336 mm 2 ) total membrane area compared to the original morphology (o). Axons in LGN INs are generally very thin, and could not be identified in the morphology data. Given the small surface area of the axon, it is not expected to affect somatic input/output-data significantly. All simulations were carried out using the NEURON simulation environment [57]. The model files are available for public download from the ModelDB section of the Senselab database (http://senselab.med.yale.edu).
Passive properties. The axial (cytoplastic) resistivity (R a ) was taken from the literature, R a = 113 Vcm [35], which also was close to the value used in earlier models of (passive) dLGN dendrites (R a = 100 Vcm in [33,34,36]). The reason for using a fixed value for R a is that this parameter is not well constrained by somatic voltage recordings, as the response has a low sensitivity to this parameter [34,58]. The remaining passive properties, the membrane capacitance (C m ), the membrane resistance (R m ), and a leakage current specified by its conductance (g pas = 1/R m ) and its reversal potential (E pas ) were estimated for each neuron by a manual trial and error process, to obtain a qualitative fit of all responses to hyperpolarizing or small depolarizing current injections in the soma. The fit procedure also included the hyperpolarization-activated cation current (I h ) and a lowthreshold calcium channel (Ca T ), which are active around rest and can influence the input resistance (R IN ) and membrane time constant (t m ) of the neuron [35]. We found that the two neurons (IN1 & IN2) had different passive properties. The final set of parameters are listed in Table 1, and are within the physiologically relevant ranges suggested earlier [59]. Note that E pas is different from the resting potential (V rest ), due to active conductances that are nonzero around rest.
Ion channel kinetics. Seven different active ion channels were included in our model, the selection of which was based on previous literature and our own simulations (see Results). The ion channels include the traditional Hodgkin-Huxley sodium-and delayed-rectifier potassium channels (Na and K dr ), a hyperpolarization-activated cation channel (I h ), a low-threshold, T-type calcium channel (Ca T ), a highthreshold, L-type calcium channel (Ca L ), a medium-duration, calciumdependent afterhyperpolarization channel (I AHP ), and a long-lasting calcium-activated non-specific cation channel (I CAN ). Ion channel kinetics is summarized in Figure 1B-H. Ion channels were inserted in the somatic and dendritic membrane, and modeled at a temperature of 36uC. Standard Hodgkin-Huxley formulation was used [60], but with the Goldman-Hodgkin-Katz formulation [61] for calcium channels instead of a reversal potential-based current. The included ion channels are described below. Experiments have characterized a slow hyperpolarization-activated cation conductance (I h ) in rat INs [41]. The kinetics of I h in that study was too slow to account for observations in our current-clamp data ( Figure 2). We therefore performed voltage clamp experiments to measure I h , and used these measurements to determine the activation kinetics (see Results). We assumed a reversal potential of -44 mV, as was earlier found in rat INs [41]. APs were generated using standard Hodgkin-Huxley-type sodium (Na) and potassium (K dr ) channels. The channels we used were originally adapted for hippocampal neurons [62], but have since then been used successfully in modeling a variety of cell types [63]. In order to be consistent with the threshold for AP initiation in the data sets, the voltage dependence of activation and inactivation was shifted compared to the original values [62]. As IN1 and IN2 had different firing thresholds, the Na kinetics was shifted +10.4 mV (i.e. in the positive direction along the voltage axis) in P1, and +12.7 mV in P2, while the K dr kinetics was shifted +11.8 mV in P1 and +15 mV in P2. The kinetics for P1 is shown in Figure 1B-C. The low-threshold, T-type calcium channel (Ca T ) needs hyperpolarization to lift the inactivation. Therefore, it is most active shortly after hyperpolarization, and then generates a calcium-dependent depolarization envelope upon which a (rebound) burst of APs may ride [36,38,43,44,64]. Ca T kinetics was adapted to recent voltage clamp data for INs [44], and corrected for temperature differences (24-36uC) using Q10 values of 3.0 and 1.5 respectively, for activation and inactivation [36]. Activation kinetics had to be shifted +8 mV in order to prevent excessive Ca T active at rest, and in order to agree with the observed bursting behavior in both example neurons (IN1 and IN2). A high-threshold, L-type calcium channel (Ca L ) opens mainly during APs. Ca L regulates tonic firing, mainly by increasing intracellular calcium levels that trigger calciumdependent potassium currents such as I AHP [36], but may also be important for dendritic signal propagation, and for triggering dendritic GABA-release [53]. Ca L kinetics was adapted from [36], with activation kinetics shifted +7 mV in order to avoid substantial Ca L activation at potentials below the AP-firing threshold. A voltageinsensitive medium-fast afterhyperpolarization current (I AHP ) was used to modify spike frequencies for depolarizing current injections. Medium-fast I AHP (reviewed in [65,66]) is calcium dependent, reaches half activation in the 400-800 nM range, has a relatively fast time to peak (1-5 ms), and decays with a time-course dependent on the amount of calcium influx [65,66]. The kinetics for I AHP was modified from [67], so that I AHP reaches half activation in the appropriate range ([Ca] i,half = 435 nM) with a time constant of about 3 ms ( Figure 1G). A long-lasting calcium-activated nonspecific cation current (I CAN ) is known to provide afterdepolarization in other thalamic neurons [68,69]. I CAN is present in INs [36,52], where it may be involved in prolonged bursts [36]. The kinetics for I CAN was taken from [36].
Ion channel distribution. Quantitative data on the subcellular ion channel distribution in INs is only available for Ca T and Ca L [43,50,51]. For the remaining ion channels we assumed a uniform distribution over the dendrite (with one exception), generally with a different density from that in the soma. The ratios between dendritic-and somatic ion channel density (g dend / g soma ) were based on relevant literature, as specified below. Experiments suggest that I h channels play an important role in dendrites of INs [52]. We make the same assumption as in a previous model for TC neurons [32], that I h channels are uniformly distributed over the soma and dendrites (g dend /g soma = 1).

Recent experiments have shown that dendritic AP propagation in
INs is supported by Na and K dr channels [53,70]. However, increases in the rise time and duration of the AP waveform during propagation, indicate that dendritic densities may be lower than those in the soma [70], which is also the case in TC neurons [71]. We assumed that dendritic Naand K dr densities were 10% of those in the soma (g dend /g soma = 0.1). The higher somatic densities may also reflect contributions from the initial axonal segment which was not modeled explicitly. In line with experimental data, we assumed that the density of dendritic Ca Lchannels was 25% of that in the soma, and that the dendritic density did not decrease with distance from soma [51]. Ca Tchannels in the dLGN are preferentially expressed in the dendritic membrane [43,50]. Their distribution may be important for determining the threshold value of the lowthreshold calcium spike [48,72]. According to experiments, the density of Ca T -channels increases linearly with distance from the soma, and reaches 239% of the somatic value at 60 mm [50], which is approximately 4%/mm. We therefore assumed a distance-dependent Ca T -density, given by g dend (x) = g soma ?(1+0.04?x), where x [mm] denotes the distance from soma. Medium-fast I AHP channels are mainly influenced by calcium entering through high-voltage-activated calcium channels [66], and due to their relatively short time to peak, it has been argued that I AHP channels are located close to the point of calcium influx [65]. In Hippocampal neurons, I AHP channels are selectively coupled to Ca L channels [73], as has also been suggested for INs [36]. We obtained such a functional coupling by assuming that I AHP and Ca L both have their peak density in the somatic region, using the distribution g dend /g soma = 0.1 for I AHP . As Ca T channels are preferentially expressed in the dendrites, they will have a smaller impact on I AHP . I CAN conductances were not apparent in dissociated cells, suggesting a dendritic location for these ion channels [36]. Studies have suggested that I CAN is mainly activated by calcium entering through Ca T channels [36,43,74]. This may be explained by the two channel types having similar spatial distribution. As for Ca T , we thus assumed an I CANdistribution given by g dend (x) = g soma ?(1+0.04?x), where x [mm] denotes the distance from soma. Intracellular calcium dynamics. The intracellular calcium level (relevant for the gating of I AHP and I CAN ), was modeled as a leaky integrator [27,75,76]: The factor a converts the calcium currents I CaT and I CaL to a calcium-concentration increase in a small volume immediately inside the membrane. We specified a = 0.0155 mmol/(cm?C) in all compartments, which gave rise to calcium levels between 50 nM and 100 nM around resting potential (depending on resting potential and Ca T conductance), and up to about 2 mM during rapid spiking/bursts, which is within the concentration ranges obtained in similar treatments (e.g. [28,76,77]). The various extrusion and buffering mechanisms are described collectively by a first-order decay process. We chose a basal calcium concentration [Ca] rest = 50 nM and a decay rate t Ca = 50 ms, which give rise to similar calcium decay as in the previous IN model [36].
Synaptic input. About 1/4 of the total synaptic input to INs is GABAergic, most likely from other INs and from interneurons in the perigeniculate nucleus and the thalamic reticulate nucleus [78,79]. The GABAergic responses in INs are predominantly mediated by GABA A -receptors [40,79]. GABAegic terminals on INs form equal numbers of synapses on dendrites and dendritic appendages, and somewhat fewer in the soma region [78]. We therefore implemented a physiologically plausible hyperpolarizing input through GABA A -receptors that were randomly distributed in the dendritic tree. The synapses were modeled as a sum of exponentials with rise and decay times of 0.5 ms and 5 ms, respectively, based on values found for GABA A -receptors in other neurons [67,80,81]. For the reversal potential, we used -82 mV, which has been measured for the GABA A -response in INs [79].

Results
Empirical data for constraining the model IN1 had a resting potential of -63 mV ( Figure 2A1). The onset of a strong hyperpolarizing current injection (-150 pA) made the membrane potential drop rapidly to a hyperpolarized peak value, before it increased to a less hyperpolarized plateau value. This initial sag is a trademark of the I h current. IN2 had a resting potential of -69 mV ( Figure 2A2). The initial sag for hyperpolarizing current injections was less pronounced in IN2 than in IN1.
Depolarizing stimuli gave rise to an initial, transient response (characterized by a high AP-firing frequency), followed by regular AP-firing at lower frequency. In IN2 the initial response tended to be distinct and burst-like (see e.g. Figure 2A2, 40 pA), whereas IN1 tended to have a more gradual transition between the initial response and the regular AP-firing (see e.g. Figure 2A1, 70 pA). We did not quantify these differences, but use the term initial burst for all initial responses that can be distinguished from the later, more regular AP firing. IN1 needed stronger depolarizing input ($55 pA) to initiate AP-firing than IN2 ($40 pA). However, IN1 had the steepest I/O-curve (see below), so that both neurons had about the same firing frequencies when stimulated at 70pA.
In one experiment, small current injections were used to hold the membrane potential at a depolarized value (-57 mV in IN1 and -58 mV in IN2). In this case, strong hyperpolarizing current pulses (-150 pA) were followed by a rebound burst in IN1 ( Figure 2B1), but not in IN2 ( Figure 2B2).

Voltage dependence of I h -activation
The I h -current in rat INs has no calcium dependence [82] and has been measured and thoroughly described [41]. A comparison of the time course of the initial sag for recordings from rat INs [36,41] with ours from mouse INs, especially IN1 ( Figure 2A1), indicated that I h may differ between the two species. We therefore recorded the I h response to different hyperpolarizing voltage steps in our mouse INs ( Figure 2C), and used these data to estimate the I h kinetics. We derived I h -activation curves for mouse INs using the following procedure [41]: With I max denoting the maximum observed amplitude throughout the trial (e.g., about -180 pA in the dataset shown in Figure 2C), and I inf denoting the steady-state value of the response for a given command potential (V), the ratio I inf /I max was plotted against V for three data sets ( Figure 3A). The data points were then fitted by a Boltzmann curve ( Figure 3A; full line), given by: I inf /I max = (1+exp((V-shift)/stp)) 21 , where shift is the potential at half-inactivation, and stp determines the steepness of the activation curve. Using the inbuilt optimizer fminsearch in Matlab, we obtained the optimal parameters shift = -96 mV and stp = 10 mV. The corresponding values estimated for rat INs were -79 mV and 7.4 mV [36,41].
In order to determine the voltage dependence of the time constants for steady-state activation, the current traces (interval from 0-1000 ms after stimulus onset in Figure 2C) were fitted by exponential curves on the form I h = (1-exp(-t/t h ))?I inf . Only a single exponential term was used, as it yielded a good fit. In this way we estimated the time constant (t h ) at each command potential as plotted in Figure 3B. The data points for t h were in turn fitted with a bell shaped curve (as in [41]) of the form: t h (V) = exp[(V+a1)/a2]/ (1+exp[(V+a3)/a4]), with V measured in mV and t h in ms. Using fminsearch, we obtained the optimal fit for [a1, a2, a3, a4] = [250, 30.7, 78.8, 5.78]. Note that the maximum value of t h is about 200 ms, while the corresponding value found for rat INs was about 1000 ms [36,41]. The faster kinetics was in good agreement with the time course of the initial sag (see Figure 4).    conditions (Figure 2A-B). As an additional constraint we assumed that the two neurons contained the same types of ion channels and had the same axial resistivity (R a ). This means that we aimed to explain the differences between IN1 and IN2 in terms of differences in three passive parameters (R m , E pas , C m ), seven parameters representing the density of ion channels (g h , g Na , g Kdr , g CaT , g CaL , g AHP , g CAN ), and two parameters (Sh Na , Sh Kdr ) representing shifts in the activation/inactivation curves of Na and K dr along the voltage axis. Shifts in the kinetics of the APgenerating currents were allowed to be free parameters as we observed clear differences in spiking threshold between neurons IN1 and IN2. In order to limit the number of free parameters, the voltage and calcium dependence of the remaining ion channels were assumed to be identical in the two neurons.

Models of two INs
The experimental protocols for the two example neurons (Figure 2A-B) were replicated in the simulations shown in Figure 4, where we show that many significant features of the experimental results are reproduced by the model using the two sets of parameters (P1 and P2, Table 1). P1 had resting potential -63 mV, and P2 had resting potential -69 mV (as in the empirical data sets). Stimulus intensities between -150 pA and 20 pA gave rise to subthreshold responses in both models. The higher response amplitudes in P2 were due to the significantly higher membrane resistance found for this neuron. Initial sags in P1 and P2 for strong hyperpolarizing current injections (-150 pA) were due to I h . Simulated action potentials had a width (at half max response) of about 0.4 ms in P1 as well as in P2 ( Figure 5A), which is typical for INs [39,40], and agreed well with the AP width of IN1. The APs elicited by IN2 were somewhat broader, and had a width of about 0.6 ms. We did not change the kinetics of Na and K DR channels to account for the variability in AP shapes. The mechanisms behind other characteristics in the response patterns are discussed in the following subsections.

I/O-curves
As in the data sets (Figure 2A-B), depolarizing current injections to the model (55-70 pA to P1 and 40-70 pA to P2) gave rise to an initial burst (or a few APs with short intra-spike intervals) followed by regular activity with a lower AP firing frequency ( Figure 4A).
The slope of the I/O curve was mainly regulated by the interplay between calcium entering through Ca L channels and a single, calcium activated potassium channel (I AHP ). As the high-voltage activated Ca L channels open during APs, the intracellular calcium concentration will accumulate during high firing frequencies, so that also I AHP increases with firing frequency. In this way the Ca L /I AHPmechanism flattens the I/O curves, as has been well described in earlier modeling studies (e.g. [83]). Without this regulatory mechanism, the firing frequency (as resulting from the Na and K dr channels) was generally too high, and the I/O curves were too steep compared to the data in Figure 2A (results not shown).
With parameters as in Table 1, the I/O curves of P1 and P2 agreed well with the data sets, being steeper in P1, as the conductance values for both these channels (g CaL and g AHP ) were higher in P2 ( Figure 5B). In order to investigate the impact of I AHP , we interchanged the conductance values (g AHP ) between P1 and P2, leaving all other parameters at their original values. The high g AHP resulted in a lower spike frequency in P1 ( Figure 6B1), and in a threefold reduction in the slope of the I/O-curve (from the original ,0.8 spikes/pA to ,0.25 spikes/pA (results not shown)). Correspondingly, the low g AHP increased the spiking frequency in P2 ( Figure 6B2), and made the I/O curve three times steeper (from the original ,0.3 spikes/pA to ,0.9 spikes/pA).

Initial bursts
The initial bursts for depolarizing current injections were well reproduced by the model (Figure 4A), and depended mainly on Ca T . The I h conductance had negligible impact on these initial bursts (results not shown), but the bursts vanished when g CaT was set to zero ( Figure 5C), and became stronger when g CaT was increased ( Figure 5D). Despite P1 having the highest g CaT , the initial bursts were most pronounced in P2. This may seem counterintuitive, but we found that this is mainly due to differences in resting membrane potential. The lower resting potential in P2 means less inactivation of Ca T at rest, and compensates for the lower g CaT . Note that our Ca T kinetics ( Figure 1D) was shifted +8 mV compared to the empirical data set for rats [44]. Using the Ca T activation from that empirical study would thus give even stronger bursts (Figure 2A).

Transition between initial bursts and regular firing
Due to the high intra-burst firing frequencies, initial bursts gave rise to strong I AHP activation, resulting in a period of afterhyperpolarization which distinctly separated initial bursts from the regular AP firing. The afterhyperpolarization was most pro-nounced in P2, as P2 had the higher g AHP (Figure 4A), but became most pronounced in P1 when g AHP was interchanged between the two parameterizations ( Figure 6B).
Previous experiments have shown that some INs respond to depolarizing input by periodic bursting [42,52]. Although periodic bursting was not observed in our experiments (IN1 and IN2), we ran test simulations to see if the mechanism that explained the initial bursts and subsequent afterhyperpolarization in IN1 and IN2 also could give rise to periodic bursting. Using the parameter sets P1 and P2 as a starting point, an increase in g CaT (by a factor 2) increased the intensity of the initial burst, but also the overall AP firing frequency, suggesting a nonzero Ca T -activity throughout the stimulus period. In P2, the increased g CaT also gave rise to a periodic firing of pairs of APs (Figure 6D2), indicating a periodic interplay between Ca T -driven bursts and I AHP -driven afterhyperpolarizations. By also increasing g AHP (by a factor 2), both the afterhyperpolarization (triggering the bursts) and the bursts (triggering the afterhyperpolarization) became more intense, and periodic bursting was obtained in both P1 and P2 ( Figure 6E). The model thus predicted that the interplay between Ca T and I AHP could explain the periodic bursting observed in some INs [42,52], and that periodically bursting neurons have high Ca T and I AHP conductances relative to neurons that do not show this behavior (e.g. IN1 and IN2).

Rebound bursting during physiological conditions
Small current injections (22 pA in P1 and 12 pA in P2) were used to shift the membrane potential from rest to the holding potentials of -57 mV in P1, and -58 mV in P2. When held at these depolarized potentials, strong hyperpolarizing current injections (-150 pA) were followed by a rebound burst in P1, but not in P2 ( Figure 4B).
In order to elicit bursts, a preceding hyperpolarization of the membrane potential is often required [36,53]. However, the somatic current injections used in vitro to study this effect do not occur in vivo. For example, the -150 pA current injections used in our experiments gave rise to a membrane potential much more hyperpolarized than the typical GABAergic reversal potential. We therefore simulated a more realistic, synaptic GABAergic input, to investigate whether our model could elicit rebound bursts under more realistic conditions. We found that 50 synapses of moderate strength (maximum conductances of 1 nS), activated with 10 ms intervals over a time period of 300 ms, reproduced the essential findings from the current-clamp experiments. When held at normal resting potentials, no rebound burst was elicited in either of the parameterizations ( Figure 7A). However, when held at -57 mV, the series of synaptic activations provoked a rebound burst in P1, but not in P2 ( Figure 7B).
Rebound bursts in various neurons are normally mediated by Ca T and/or I h channels [39,64,[84][85][86]. We used the simulation setup with synaptic input to explore the contributions from these two ion channels to the rebound bursts in INs. In the original parameterizations (Table 1), P1 had higher values than P2 for both g CaT and g h , which explains why P1 and not P2 elicited rebound bursts. When g h was interchanged between P1 and P2, leaving all other parameters as in Table 1, the rebound burst in P1 was reduced to include only a single AP, while the (sub-threshold) rebound response in P2 became more pronounced ( Figure 7C). Similar results were obtained when g CaT was interchanged between the two parameterizations: The rebound response in P1 fell below the AP firing threshold, while the (sub-threshold) rebound response in P2 became stronger ( Figure 7D). Finally, if both parameters (g h and g CaT ) were interchanged between P1 and P2, also the rebound responses were entirely interchanged: P1 showed a small sub- Figure 6. Mechanisms behind initial bursts and regular spiking. All panels show responses to a 70 pA depolarizing current injections. As a reference, the responses in (A) use the original parameters for P1 and P2, and are the same as in Figureô 4. In P2, the initial burst and the regular spiking are separated by a pronounced afterhyperpolarization (A1), while in P1 the transition between the initial burst and the regular spiking is more gradual (A2). These characteristics were interchanged between P1 and P2 when g AHP was interchanged (i.e., multiplied/ divided by a factor 2.03 in P1/P2) between the two parameterizations (B). Initial bursts were eliminated when g CaT was set to zero (C), and became stronger when g CaT was increased by a factor 2 (D). Increasing g AHP and g CaT by a factor 2 gave rise to periodic bursting in both neurons (E). The scale bar applies to all panels. When conductance values were changed, the resting potential was kept at the original level by small compensatory current injections. doi:10.1371/journal.pcbi.1002160.g006 threshold rebound response, while P2 elicited a pronounced rebound burst, resembling that originally seen in P1 ( Figure 7E). This suggests that Ca T and I h are of comparable importance for rebound bursts in mouse INs.

Dependence on morphology
To investigate how sensitive our results are to the IN morphology, we ran test simulations using three additional IN morphologies with (a) a similar, (b) a smaller and (c) a larger total membrane area compared to the original morphology (o). We replaced the original morphology (o) with the new morphologies (a-c), but kept all other model parameters fixed at the values in Table 1. As the ion channel densities (i.e. conductances per mm 2 ) were kept fixed, the new morphologies (a-c) corresponded to INs with (a) a similar, (b) a higher, and (c) lower input resistance compared to the original morphology (o).
The AP shape was not strongly affected by the morphology changes, except from small variations in the afterdepolarization ( Figure 5A). The I/O curve in cases (o) and (a) nearly coincided ( Figure 5B). As expected from the differences in total input resistance, morphology (b) gave rise to a I/O curve that was steeper and shifted in the negative direction along the current axis, whereas morphology (c) gave rise to an I/O curve that was flatter and shifted in the positive direction along the current axis compared to the cases (o) and (a) ( Figure 5B). Although occurring at different current input levels, the characteristic responses of P1 and P2 did not change qualitatively with morphology. Depolarizing current injections of sufficient intensity gave rise to initial bursts followed by regular AP firing, and with the initial bursts being most pronounced for the parameterization P2. Furthermore, strong hyperpolarizing current injections (-150 pA in case of morphology (o), (a) or (b) and -200 pA in case of morphology (c)) were followed by rebound bursts when using parameter set P1 at a holding potential of -57 mV, but not when using the parameter set P2 at a holding potential of -58 mV. The essential firing patterns of P1 and P2 were thus preserved under changes of morphology, and the I/O curve was always steeper for parameterization P1 ( Figure 5C).

Discussion
We have developed a compartmental IN model which stands out from previous models in several ways: Firstly, the model was constrained by a broader range of I/O data than previous models of INs. It was able to exhibit a range of the response patterns characteristic for INs, including (i) initial sags, (ii) rebound bursts, (iii) tonic AP-firing, (iv) initial bursts/spike-time adaptation, and (v) periodic bursting. The mechanisms that we used to model these features have been described earlier in different works (e.g. [36,63,64,83,84,86,87]. By constraining the conductances of the different ion channels to I/O data from single neurons under eight different experimental conditions, including quantitative data on firing frequency vs. stimulus amplitude for depolarizing stimuli (I/ O curves), we were able to make more specific predictions on the contributions of each mechanism. Secondly, we used a much more realistic morphology compared to earlier IN models that include active conductances, opening the way to future studies investigating how different dendritic regions (that generally receive input from different sources [78]) process local I/O operations. Thirdly, dendritic signal propagation also depends on active dendritic properties. As the dendrites of INs are both pre-and postsynaptic, dendritic ion channels will shape both incoming and outgoing signals, and will have an important impact on the signaling between INs and other neurons. Unlike previous IN models [33][34][35][36], our model includes a spatial distribution of ion channels over the somatodendritic membrane that is consistent with the available empirical data.

Ion channels
The set of seven ionic conductances that we used to fit the experimental findings (Table 1) is the same as in the previous model by Zhu et al. [36], but with kinetics that were updated to account for recent findings for activation/inactivation kinetics and somatodendritic distributions, and with conductance values constrained by a broader range of I/O data. Although these seven channel types were successful in reproducing the observed spiking patterns, we cannot exclude the presence of additional mechanisms, either overlapping with an included conductance type in terms of their action on the firing properties, or with a minor or negligible effect on the somatic response observed in the current clamp experiments.
Several of the included ion channels have well documented roles in INs, including the role of Ca T in burst generation [36,39,44], the role of I h in generating initial sags [41] and the role of Ca L in increasing the intracellular calcium concentration [51,53]. I CAN is involved in making the dendrites leakier in connection with Figure 7. Mechanisms behind rebound bursting. All panels show responses to strong hyperpolarizing synaptic input (50 synapses, activated with 10 ms intervals, during 300 ms). Rebound bursts were not elicited by P1 (A1) or P2 (A2), starting from their normal resting potentials of -63 mV and -69 mV, respectively. From a -57 mV holding potential, P1 (B1) elicited a rebound burst following hyperpolarization, while P2 did not (B2). Interchanging either the I h conductance (C) or the Ca T conductance (D) between P1 and P2, leaving all other parameters the same, reduced the rebound response in P1, and increased it in P2. Interchanging both the I h and Ca T conductances between P1 and P2, also interchanged their bursting abilities completely (E). The scale bar applies to all panels. When conductance values were changed, the resting potential was kept at the original level by small compensatory current injections. doi:10.1371/journal.pcbi.1002160.g007 cholinergic modulation [52]. Its influence on the somatic response pattern of INs is less clear, although a previous modeling study has suggested that I CAN may be involved in generating plateau potentials and prolong bursts [36]. Due to the high calcium sensitivity of this channel [36,88], we found that even small depolarizations of the membrane gave rise to a tonically active I CAN . In our model, the main function of I CAN was to reduce the magnitude of the depolarizing current required for the neuron to reach AP firing threshold.
In comparison to somatic voltage recordings from rat INs [39,41], our recordings from mouse INs show more pronounced initial sags for strong hyperpolarizing current injections (see Figure 2A, -150 pA stimuli). This suggests that the kinetics of I h may differ between INs in rats and mice, as is the case in CA1 pyramidal neurons [89].This we confirmed by measuring the voltage dependence of I h in three mouse INs (Figure 3). Simulations with I h kinetics based on our own measurements not only agreed better with the sag-shape in the mouse INs, but also predicted that the impact of I h on rebound burst generation was comparable to that of Ca T (Figure 6). This differs from the situation in rat INs, where rebound responses are mainly due to Ca T , with no measurable contribution from I h [39]. Conversely, in CA1 pyramidal neurons it has been shown that rebound spiking can be generated by I h alone [86]. However, a joint involvement of Ca T and I h in burst generation, were found to underlie intrinsic rhythmic bursting in subpopulations of TCs during inattentiveness in guinea pigs [64,84]. Intrinsic rhythmic bursting was not observed in our neurons.
The presence of Ca L -conductances in IN dendrites [51,53,54] makes it likely that also inhibitory, calcium-dependent mechanisms (such as I AHP ) are present. The role of I AHP in INs has not been previously documented. Our model predicted that the interplay between Ca L and I AHP conductances was sufficient for explaining the modulation of the I/O curves, although additional mechanisms could be involved. For example, a slowly activating potassium channel (K M ) with a high threshold and no inactivation gave similar results to our Ca L and I AHP mechanism, but without any calcium dependence (simulations were made using K M kinetics taken from [90], results not shown). However, no clear functional role could be assigned to K M other than that covered by I AHP . We thus did not explore this issue further, since I AHP gave a better agreement than K M with the time course of the intra-spike membrane potential and, especially, with the afterhyperpolarization following the initial bursts in P2. However, the possibility that K M and I AHP have overlapping functions in regulating the spiking frequency cannot be excluded. This could be experimentally tested by blocking the respective channels.
We presented two parameterizations of the model (P1 and P2), which reproduced the electrophysiological properties of two different INs (IN1 and IN2). Rebound bursts as those generally observed only in a small subset of INs [39] were elicited by P1 but not P2. On the other hand, P2 elicited more pronounced initial bursts than P1 when exposed to depolarizing stimuli. P2 also had a less steep I/O curve and required weaker depolarization than P1 in order to reach AP-firing threshold. Our simulations showed that differences between P1 and P2 in terms of response properties and preferred input conditions arose from relative differences in specific conductances. Our model thus supports the idea that conductances values (i.e. channel density) in different subgroups of INs may be tuned in such a way as to optimize network operation under different input conditions. Under in vivo conditions, changes in input conditions (e.g. in synaptic input and/or shifts in membrane potential) may be mediated by mGlu5-receptor activation, GABAergic input from the reticular nucleus, or cholinergic modulation [4,53,54,91,92].

Dendritic signal propagation
Data on the somatic voltage responses to somatic current injections do not uniquely determine the distribution of passive and active properties in the dendrites (see e.g. [48,93]). Assumptions on the somato-dendritic distributions of Ca T and Ca L in our model were therefore based on calcium imaging data [50,51]. This may be very useful in future studies, as the specific sub-cellular localization of different types of calcium channels may be particularly important for their specific functional role [43]. The assumed distributions of the remaining ion-channels were based on what we judged to be the most relevant literature to date. With the assumption that dendritic Na and K dr densities were 10% of those in the soma, simulations showed that APs got broader during propagation in the dendrites, whereas the amplitude did not get significantly attenuated (results not shown). This is, at least at a qualitative level, in good agreement with recent experiments using voltage-sensitive dye [70].
We cannot exclude the possibility that important aspects of dendritic signaling are not captured by our model at this stage. For instance, there is an uncertainty regarding the distribution of our dendritic Ca T channels. We based our Ca T distribution on electrophysiological data [50], which indicated that the Ca T density increases with distance to soma, while other, anatomybased studies, have indicated a uniform distribution of dendritic Ca T channels [94]. In test simulations we showed that we could obtain essentially the same results as in Figure 4 and Figure 5 also with a uniform Ca T distribution, simply by rescaling the total Ca T conductances (results not shown). However, although it may not be crucial for an IN's response to somatic current injections, the Ca T distribution will likely influence aspects of dendritic signaling, such as the probability for dendritic GABA-release.
A related source of uncertainty concerns the presence of dendritic A-type potassium channels (K A ), which, suggested by early studies, counteracted the Ca T -channels in the dendrites of INs, and suppressed bursting [38]. This mechanism has, however, not been observed consistently, and several experiments have reported bursting INs [36,39,40,43]. For thalamic reticular neurons, the ability versus inability to burst was rather explained by a varying density of Ca T -channels [95], as also fits well with our findings. K A channels are often most densely present in distal dendrites and might have an influence on backpropagating action potentials [96][97][98]. However, recent experiments on INs found that attenuation of dendritic Ca-signals was relatively small [53]. As the importance of K A is unclear, these channels were not included in our model. Test simulations, using a moderate K A density (channel kinetics from [97]) in the dendrites, did not affect the somatic response to somatic stimuli significantly (results not shown). Such channels could therefore be readily added to the model if future experiments identify a clear functional role of K A in IN dendrites.
Parts of the distal dendrites of INs form so called triadic synapses with axons from retinal ganglion cells and dendrites from TC neurons [4,99]. In these triads, the IN terminals are, at the same time, postsynaptic to retinal input, and presynaptic to TC neurons. The conditions for GABAergic release from IN dendrites are not fully known, but may depend on intracellular calcium levels which are elevated by (backpropagating) sodium spikes as well as signals evoked by local synapses [4,53,54,99]. It is known that cholinergic modulation reduces the membrane resistance of INs, through the M2-receptor mediated activation of I h , I CAN and a linear, unspecified potassium current [52]. During sleep, when the cholinergic tone tends to be low, interneurons are likely to be electronically compact. INs may then provide long range inhibition, as synaptic input at one location in the dendritic tree then result in GABA release throughout the dendritic and axonal arbors. During awake states, when the cholinergic tone is high, distal dendritic regions may become electronically isolated from each other and from the soma. During these conditions, the triads may function as independent units, being excited by the presynaptic retinal afferents and then directly inhibiting only the postsynaptic dendrites of TC neurons [33,52,100,101]. In our model, somatic action potentials successfully invaded distal dendritic region, and it is thus likely that our parameterizations (P1 and P2), correspond to conditions with a low cholinergic tone.
The present model includes dendritic sodium-, potassium-, and calcium channels, and at least some of the mechanisms that are affected by cholinergic modulators. It computes the time course of the intracellular calcium levels in each compartment, and contains essential mechanisms for addressing signaling in IN dendrites on a fine spatiotemporal scale. In a network model including the dLGN circuitry, this will be of paramount importance for simulating the interactions between TCs and INs, as well as inputs to these cells from retina, cortex, thalamic reticular nucleus and possibly modulatory input from the brain stem.