## Abstract

A temporal sensory code occurs in posterior medial (POm) thalamus of the rat vibrissa system, where the latency for the spike rate to peak is observed to increase with increasing frequency of stimulation between 2 and 11 Hz. In contrast, the latency of the spike rate in the ventroposterior medial (VPm) thalamus is constant in this frequency range. We consider the hypothesis that two factors are essential for latency coding in the POm. The first is GABA_{B}-mediated feedback inhibition from the reticular thalamic (Rt) nucleus, which provides delayed and prolonged input to thalamic structures. The second is sensory input that leads to an accelerating spike rate in brain stem nuclei. Essential aspects of the experimental observations are replicated by the analytical solution of a rate-based model with a minimal architecture that includes only the POm and Rt nuclei, i.e., an increase in stimulus frequency will increase the level of inhibitory output from Rt thalamus and lead to a longer latency in the activation of POm thalamus. This architecture, however, admits period-doubling at high levels of GABA_{B}-mediated conductance. A full architecture that incorporates the VPm nucleus suppresses period-doubling. A clear match between the experimentally measured spike rates and the numerically calculated rates for the full model occurs when VPm thalamus receives stronger brain stem input and weaker GABA_{B}-mediated inhibition than POm thalamus. Our analysis leads to the prediction that the latency code will disappear if GABA_{B}-mediated transmission is blocked in POm thalamus or if the onset of sensory input is too abrupt. We suggest that GABA_{B}-mediated inhibition is a substrate of temporal coding in normal brain function.

## INTRODUCTION

Substantial anatomical and physiological evidence indicates that visual (Guillery and Sherman 2002b), auditory (Suga et al. 1997; Zhang and Suga 1997), and tactual (Castro-Alamancos 2002b; Diamond et al. 1992; Nicolelis and Chapin 1994; Richardson 1973; Temereanca and Simons 2004) thalamic nuclei do not act as a relay of sensory input but rather process sensory information (Guillery and Sherman 2002a). In the rodent vibrissa somatosensory system, input from brain stem nuclei is transformed by three interacting thalamic nuclei (Fig. 1*A*); the posterior medial (POm) thalamic nucleus, the ventroposterior medial (VPm) thalamic nucleus, and the reticular (Rt) thalamic nucleus (for a review of anatomy, see Deschenes et al. 1998; Diamond 1995; Kleinfeld et al. 1999). The POm and VPm nuclei receive afferent input primarily from trigeminal nucleus principalis (Pr5) and spinal nucleus interpolaris (Sp5I) (Veinante et al. 2000) and project reciprocally to the same area of the Rt nucleus (Crabtree and Isaac 2002; Crabtree et al. 1998). The Rt nucleus, in turn, provides feedback inhibition to the thalamocortical neurons in both POm and VPm nuclei via fast GABA_{A}- and slow GABA_{B}-mediated synaptic currents (Von Krosigk et al. 1993). The possible mixing of signals between POm and VPm nuclei by Rt neurons forces consideration of the full POm-Rt-VPm circuit as a means to delineate the thalamic dynamics.

Despite the apparent symmetry of the organization of the VPm and POm thalamic pathways (Fig. 1*A*), recent experiments demonstrated a critical difference in the steady-state dynamics of the two nuclei in response to passive rhythmic movement of the vibrissae (Ahissar et al. 2000, 2001; Sosnik et al. 2001). The time course of the motion of the vibrissae during whisking, albeit not the mechanics (Szwed et al. 2003), was approximated by the application of periodic air puffs to one or two rows of vibrissae. The latency and rate of spiking for brain stem units were observed to be independent of the stimulus frequency for frequencies between 2 and 11 Hz (Ahissar et al. 2000; Deschenes et al. 2003; Jones et al. 2004; Sosnik et al. 2001), which incorporates the range of natural exploratory whisking (Berg and Kleinfeld 2003; Welker 1964). While the response latency of units recorded in VPm thalamus was also nearly independent of stimulus frequency (Ahissar et al. 2000; Hartings and Simons 1998; Hartings et al. 2003; Sosnik et al. 2001), the response latency of units recorded from the POm thalamus increased considerably with frequency for stimuli with a 50-ms, albeit not a 20-ms rise time (2000, Ahissar et al. 2001; Sosnik et al. 2001). Thus POm but not VPm thalamus can code stimulus frequency in terms of latency. In addition, there is a clear coding of stimulus frequency in the rate of VPM neurons.

We consider the hypothesis that the essential determinants of the latency coding observed in POm thalamus depend on both features of the internal circuitry, i.e., the strong facilitation as well as delayed and prolonged response of GABA_{B}-mediated synaptic transmission (Kim et al. 1997), and features of the external stimulus, i.e., the gradually increasing spike rate of the brain stem input to thalamic nuclei. Our challenge is to address this hypothesis in terms of a neural-based model and account for the differences in the response of neurons in POm versus VPm thalamus to periodic stimulation. We ask: why does the latency in the response of POm neurons to prolonged stimuli increase substantially with frequency during steady state conditions? Why is there a distinct difference in latency coding for POm versus VPm thalamus? And why is there no latency coding for brief stimuli? A distinctive feature of our approach is the use of reduced models that are amenable to analytical treatment as a means to gain intuition into the role of synaptic versus external inputs as well to gain insight into the stability and sensitivity of solutions of the model.

## MODEL

### Architecture of the POm-Rt-VPm circuit

The architecture of the model includes the three thalamic nuclei: POm, VPm, and Rt. The Rt nucleus receives fast excitatory input from both POm and VPm nuclei. In turn, the POm and VPm nuclei receive inhibitory feedback from Rt thalamus as a fast GABA_{A}- and delayed and prolonged GABA_{B}-mediated current (Fig. 1*A*). As a means to obtain insight into the neuronal dynamics and to develop a set of robust predictions, we examine a succession of models of increasing complexity. We begin with a reduced but analytically tractable circuit comprising only the POm and the Rt nuclei (Fig. 1*B*). Then, we investigate the effects of feedforward input from the VPm to the Rt (Fig. 1*C*). Finally, we study the full thalamic architecture that further includes feedback inhibition from the Rt to the VPm nucleus (Fig. 1*D*).

We note that the difference in latency coding between the VPm and POm thalamus is also reflected by neurons in cortical layers that receive axonal projection from these nuclei, i.e., layers IV and Va, respectively (Ahissar et al. 2000, 2001; Ahrens et al. 2002). We focus on feedback dynamics at the level of the thalamus. We test if GABA_{B}-mediated currents and the brain stem stimulus shape are sufficient to explain the observed differences in latency.

### Measured spike rates

The measured peristimulus time histograms (PSTHs) in brain stem nuclei Pr5 and SpVI and in the thalamic nuclei VPm and POm, based on population-average data from all the recorded multiunits (Ahissar et al. 2000, 2001; Sosnik et al. 2001), show different steady-state responses to periodic air-puff stimuli, each 50 ms in duration, delivered at frequencies of *f*_{stim} = 2, 5, 8, and 11 Hz (Fig. 1*A*). The population-averaged temporal responses for both the Pr5 and the Sp5I nuclei of the brain stem may be described as a dual-sloped ramp in which a brief, rapidly rising phase is followed by a prolonged, slowly rising phase. The neuronal activity then gradually decays after the activity reaches its maximum value. In contrast to the prompt response in brain stem, the activity of units in POm thalamus are delayed and then rise to their maximal value. Critically, the response rises more gradually to its peak value as the stimulus frequency, *f*_{stim}, is increased (cf. 2- and 11-Hz data in Fig. 1*A*), so that the rise time is, very roughly, inversely proportional to *f*_{stim}, whereas the maximal activity decreases with *f*_{stim}. In contrast to the case for POm thalamus, the activity of units in VPm thalamus begins with shorter delay and then rises rapidly to a peak. The rise time is almost independent of *f*_{stim}, but the maximal activity decreases with *f*_{stim}.

### Rate model

Our technical approach makes use of rate equations (Hopfield 1982; Wilson and Cowan 1973), where the instantaneous spike rate of a population of neurons is expressed in terms of the presynaptic activation for each synaptic current (Ermentrout 1994; Kang et al. 2003; Rinzel and Frankel 1992; Shriki et al. 2003). The activation parameters are denoted by subscripted *u*_{i}(*t*), where the subscript *i* is used to denote the presynaptic neuronal population. For inhibitory conductances, a second subscript, A or B, denotes whether the synapse is mediated by GABA_{A} or by GABA_{B}, respectively. Thus as examples, *u*_{POm}(*t*) describes activation of the excitatory current that originated from POm neurons and *u*_{Rt,B}(*t*) describes activation of the slow inhibitory current that originates from the Rt nucleus. We further note the aggregate measure of the strength of the synaptic connection of the projection from nucleus *j* to nucleus *i* by the synaptic conductance *g*_{i;j}. For example, *g*_{POm;Rt,B} denotes the GABA_{B}-mediated input from Rt to POm thalamus. The activation parameters *u*_{i}(*t*), multiplied by the synaptic conductances *g*_{i;j}, determine the total synaptic current to the postsynaptic cell, e.g., *g*_{POm;Rt,B}*u*_{Rt,B}(*t*) for GABA_{B}-mediated input to POm neurons from Rt thalamus.

##### NONLINEAR, DELAY DIFFERENTIAL EQUATIONS THAT DEFINE THE DYNAMICS.

The observable quantities in a network are the instantaneous spike rates for each population of neurons. This rate is denoted by *M*_{i} for nucleus *i*, e.g., *M*_{POm} for the POm nucleus. Without intrinsic adaptation, the presynaptic spiking rate *M*_{i} is approximated by the instantaneous input-output relation, or *f*-*I* curve, *M*_{i}(*t*) = [*I*(*t*)]_{+}, and is taken to be the rectification (linear-threshold) function, that is (1)

The total current *I*(*t*) is the sum of the synaptic currents and the brain stem input, i.e., *I*_{POm}(*t*) and *I*_{VPm}(*t*) for the POm and VPm nuclei, respectively, relative to the threshold for spiking, e.g., θ_{POm} for neurons in POm thalamus. Thus the instantaneous spiking rates of the three thalamic nuclei are (2) (3) (4) The POm and the VPm nuclei, but not the Rt nucleus, further possess multiplicative adaptation processes (Hartings and Simons 2000) defined through the dynamics of the variables *a*_{VPm} and *a*_{POm} (see following text). All the variables and parameters, except for time, are normalized and are unit-less in our formulation (Ermentrout 1994; Kang et al. 2003; Shriki et al. 2003). *Equations 2*–*4*, with nonnegative values of the conductances, correspond to the architecture of Fig. 1*D*.

The dynamics for glutamatergic fast excitatory synapses and GABA_{A}-mediated inhibitory synapses, but not GABA_{B}-mediated inhibitory synapses, are defined by pairs of delay differential equations (5) (6) and (7) (8) where τ̃_{G}, τ_{G}, and *t*_{G} are the synaptic rise, decay, and delay times, respectively, of the excitatory synapses, τ̃_{A}, τ_{A}, and *t*_{A} are the synaptic rise, decay, and delay times, respectively, of the GABA_{A}-mediated synapses, the variable *x*_{i} is an auxiliary variable, and the subscript *i* is POm thalamus or VPm thalamus. The synaptic delays depend only on the synaptic phenotype in our architecture.

Inhibitory synapses mediated by GABA_{B} are nonlinear and facilitating. In practice, they have a delay of ∼30–40 ms and respond much more strongly to a prolonged burst of spikes than to a brief burst (Golomb et al. 1996; Kim et al. 1997). We use a version of a nonlinear model for GABA_{B}-mediated synapses (Golomb et al. 1996) (9) (10) where τ̃_{B}, τ_{B}, and *t*_{B} are the synaptic rise, decay, and delay times, respectively, of the GABA_{B}-mediated synapses. The quadratic nonlinearity in the first term in the right-hand side of *Eq. 10* is responsible for the facilitating nature of the synapse.

##### ADAPTATION.

We incorporate an idealization of cellular adaptation in the rate equations for the VPm and POm thalamic nuclei with an adaptation time scale that is much slower than the spiking time scale. A thalamic neuron fires a small number of spikes in response to a single stimulus in the relay mode (Minnery and Simons 2003; Sosnik et al. 2001). In particular, cells in POm thalamus rarely fire more than one spike in response to a single stimulus (Sosnik et al. 2001). We therefore model adaptation as a process of inactivation and slow recovery of the “neuronal pool” that is capable of spiking (Eggert and van Hemmen 2000, 2001). Such a process can be described as a multiplicative dynamical process with two time scales, one of activation as a result of neuronal activity and one of inactivation when neurons are silent. The equations for the adaptation variables *a*_{i} (*Eqs. 2* and *3*) and the auxiliary variables *b*_{i}, identified with the activation dynamics, are (11) (12)

We assume that POm thalamus has stronger adaptation than VPm thalamus; this allows us to model the fast decay of POm activity at low frequencies in comparison to the more sustained activity of VPM neurons (Fig. 1*A*). The adaptation in the VPm neurons rises faster than that in POm neurons as a result of the higher spiking rate in the former nucleus. The activity in the Rt nucleus is more prolonged than the activity in thalamic relay nuclei (Hartings and Simons 2000), and therefore we do not introduce adaptation for the Rt nucleus.

##### STIMULUS SHAPE.

The input from the brain stem nuclei is considered as an external variable that monotonically tracks the stimulus. We do not discriminate between inputs from the Pr5 and SP5I nuclei as the average responses of neurons to air-puff stimuli are very similar in the two areas (Ahissar et al. 2000). In both dorsal thalamic nuclei, the input from the brain stem begins after a short delay that represents the onset of the stimulus and the time that the stimulus drives the thalamic nuclei. We model the “double ramp” shape of the brain stem input to VPm thalamus, *I*_{VPm}, as a piecewise linear function. The brain stem input to POm thalamus is proportional to the brain stem input to the VPm, and it is delayed to allow larger latency in POm even for low stimulus frequencies (13)

The parameters of the brain stem stimuli to the two thalamic nuclei are given in appendix a.

##### PARAMETERS.

We use the following parameters throughout the analysis unless stated otherwise: thresholds: θ_{VPm} = θ_{POm} = θ_{Rt} = 0; excitatory synapses: τ̃_{G} = 1 ms, τ_{G} = 2 ms, *t*_{G} = 0 (Golomb and Amitai 1997); GABA_{A} synapses: τ̃_{A} = 1 ms, τ_{A} = 10 ms, *t*_{A} = 3 ms; GABA_{B} synapses: τ̃_{B} = 40 ms, τ_{B} = 150 ms, *t*_{B} = 35 ms (Golomb et al. 1996; Kim et al. 1997); synaptic conductances: *g*_{VPm;Rt,A} = 0.5, *g*_{VPm;Rt,B} = 2, *g*_{POm;Rt,A} = 0.16, *g*_{POm;Rt,B} = 3.5, *g*_{Rt;VPm} = 1.2, *g*_{Rt;POm} = 1; and adaptation: *k*_{a,POm} = 0.33 ms^{−1}, *k*_{a,VPm} = 0.1 ms^{−1}, τ_{b,VPm} = 10 ms; τ_{b,POm} = 30 ms, *k*_{b,POm} = *k*_{b,VPm} = 0.05 ms^{−1}, and τ_{a,VPm} = τ_{a,POm} = 100 ms.

### Derived quantities

##### SPIKE NUMBER.

The function *M*_{POm}(*t*) (*Eq. 3*) is proportional to the instantaneous spiking rate of POm neurons. A measure that is of further utility in our analysis is the total number of spikes that are fired during a time interval by these neurons. This number is proportional to *N*_{POm}^{spikes}, where (14) where *T*_{integ} is the temporal integration window and *N*_{POm}^{spikes} has units of time because *M*_{POm}(*t*) is dimensionless in our formulation (*Eq. 3*). We choose the value *T*_{integ} = 90 ms to be smaller than the period of the highest stimulus frequency considered, i.e., (11 Hz)^{−1} = 91 ms, to enable a fair comparison between the total responses to different frequencies. A similar equation is used for the spiking rate of VPm neurons (15)

##### LATENCY.

The onset latency, *t*_{0}, is computed in our analytical treatment. We further define the time between the start of the stimulus and the time that the instantaneous spiking rate, *M*_{i}(*t*), reaches half of its maximal value at the midpoint latency, as *t*_{0.5}. This second definition is consistent with that used to quantify the latency in the experimental observations (Ahissar et al. 2000, 2001; Sosnik et al. 2001).

## RESULTS

### Qualitative expectations

GABA_{B}-mediated inhibition has slow kinetics and facilitates with increased presynaptic spiking. Under steady-state conditions, the slow kinetics cause a rhythmic input to lead to a delayed inhibition of thalamic cells. For a rhythmic stimulus, the extent of the delay will depend on the rate of rise of the brain stem input and the frequency of the stimulus. If the rise is gradual, facilitation will transform an increasing rate of stimulus-induced spiking into an increased latency in the response of thalamic cells because higher spike-rates lead to stronger inhibition. The model enables us to assess the consequences and stability of this mechanism and the specific role of the three thalamic nuclei.

Our first goal is to gain qualitative insight into the role of GABA_{B}-mediated inhibition and the stimulus shape in the production of a frequency-dependent latency. To achieve this goal, we simplify the rate model to obtained a reduced, analytically solvable model. The analysis is extended to circuits with increasing level of complexity. Finally, the outcome of the reduced model is compared with simulations of the full rate model to demonstrate that the approximations used for the reduction do not disrupt the basic dynamical mechanisms. The simulations further allow us to make a detailed comparison between experimental and model results as well as to make predictions for the outcome of future experimental investigations.

### Reduced model of the thalamic circuit

A simplified model is obtained through the inclusion of four assumptions. First, the GABA_{B}-mediated decay-time τ_{B} is larger than all the other synaptic time constants in the system. Therefore with the exception of the GABA_{B}-mediated delay (which plays a special role in the dynamics) and decay, all of the synaptic processes are considered to be instantaneous. This allows us to focus on the special role of GABA_{B}-mediated currents in shaping the thalamic dynamics. In particular, we set *t*_{POm,delay} = 0 in *Eq. 13* because this time constant is much smaller than the GABA_{B} decay time so that *I*_{POm}(*t*) = α*I*_{VPm}(*t*). Further, since the time constants of the fast excitatory synapses (τ̃_{G}, τ_{G}, and *t*_{G}) are set to zero (*Eqs. 5* and *6*), we replace *u*_{VPm} by *M*_{VPm} and replace *u*_{POm} by *M*_{POm} in the expressions for the firing rate (*Eqs. 2*–*4*).

The second assumption concerns adaptation. Because adaptation does not affect the thalamic response at the onset of activity, we simply ignore it and take *a*_{POm} = *a*_{VPm} = 0.

The third set of assumptions concern the fast conductances. GABA_{A}-mediated inhibition is ignored because this inhibition has fast decay rate, which is now taken to be instantaneous. Therefore inhibition generated in response to a given stimulus decays before the next stimulus arrives and does not participate in generating the latency. The final assumption is that, for simplicity, the thresholds θ_{i} are set to 0 for all nuclei *i*.

Based on the first assumption, the GABA_{B}-mediated conductance in the Rt nucleus is simplified by taking the onset to be instantaneous, i.e., τ̃_{B} = 0 (*Eq. 9*). This allows us to substitute *x*_{Rt}(*t*) by *M*_{Rt,B}(*t* − *t*_{B}) (*Eqs. 8* and *9*). This yields (16) *Eq. 16* defines the minimal model, from which we can deduce the spiking rates of the POm, VPm, and Rt.

To realize analytical treatment, we treat only cases where *T*_{stim} = 1/*f*_{stim} >2*t*_{B}. Furthermore, because the rise time of the brain stem stimulus-evoked response is ∼50 ms (Fig. 1*A*), which is only slightly greater than the delay of GABA_{B}-mediated synapses, *t*_{B}, we take *t*_{B} to equal the stimulus duration. Further, to demonstrate the importance of this stimulus shape for the latency in a manner that avoids the algebraic complexity associated with a “double-ramp” treatment, we consider first a triangular stimulus to mimic the slowly rising phase and next a square stimulus to mimics the fast rising phase. The shapes of stimuli are, for 0 ≤ *t* < *T*_{stim} (17) (18) For triangular stimulus (*Eq. 17*), the latency to half-maximum *t*_{0.5} is given by *t*_{0.5} = (*t*_{0} + *t*_{B})/2.

The value *N*_{POm}^{spikes}, proportional to the number of spikes fired by POm neurons, is (*Eq. 14*) (19) Our normalization is such that the value of *N*_{POm}^{spikes} varies between 0, when POm is silent, and *t*_{B}, with no inhibition to the POm nucleus (*Eqs. 17* and *18*).

### Analytical results for reduced POm-Rt circuits

##### POm-Rt CIRCUIT FOR TRIANGULAR STIMULUS: ESSENTIAL FEATURES FOR CODING STIMULUS FREQUENCY BY RESPONSE LATENCY.

We consider first a minimal circuit with feedback between the POm and Rt nuclei, in addition to input from the brain stem to POm thalamus (Fig. 1*B*). The dynamics of the GABA_{B} activation variable, *u*_{Rt,B}(*t*), are determined by a single delay-differential equation with five parameters (*Eq. 16*) and the expression for the triangular stimulus (*Eq. 17*). Three of these parameters are time constants, i.e., the time period of the stimulus *T*_{stim}, the stimulus duration *t*_{B}, and the GABA_{B} decay rate τ_{B}. The other two parameters are the POm-to-Rt excitatory conductance *g*_{Rt;POm} and the Rt-to-POm GABA_{B}-mediated inhibitory conductance *g*_{POm;Rt,B}. The dynamics of activity in POm thalamus are determined by the difference between the excitatory stimulus *I*_{POm}(*t*) and the GABA_{B}-mediated inhibition *g*_{POm;Rt,B} *u*_{Rt,B}(*t*) (Fig. 2*A*). Thus the latency for activation of POm neurons corresponds to the time when the difference between the stimulus and the inhibition becomes positive, namely (appendix b, *Eq.* B*2*) (20) The activation variable *u*_{Rt,B}(*t*) decays exponentially across the entire time period except for the interval between *t*_{0} + *t*_{B} and 2*t*_{B}, when it grows in delayed response to the POm, and therefore Rt, activity. The value of *N*_{POm}^{spikes}, proportional to the number of spikes fired by POm neurons per stimulus cycle, is the time integral of the stimulus minus the inhibition (shaded area in Fig. 2*A*). This is given by the integral of the activity in POm thalamus (appendix b, *Eq.* B*10*), that is (21) and is valid for all of the architectures we consider for a triangular stimulus. The number of spikes per cycle decreases as a function of increasing latency (Fig. 2*B*). For values of *t*_{0} that are not close to those of *t*_{B}, this decrease is almost linear.

The first step of the analysis of the system dynamics is to derive a map between *u*_{Rt,B}(0), the activation at time *t* = 0 and *u*_{Rt,B}(*T*_{stim}), the value of the GABA_{B} synaptic activation after one period. Steady-state network activity that has the same periodicity as the stimulus corresponds to a fixed point (steady-state solution) of that map. To compute this activity, we set *u*_{Rt,B}(*T*_{stim}) = *u*_{Rt,B}(0). The map in the vicinity of the fixed points has the form (appendix b, *Eq.* B*6*) (22) Because the map (*Eq. 22*) is a nonlinear function of *u*_{Rt,B}(0), the fixed point with a period of a single stimulus cycle can become unstable via a period doubling (PD) bifurcation for (23) (Strogatz 1994). We examine the fixed point and its stability graphically (Fig. 2*C*). For small values of *u*_{Rt,B}(0), the corresponding value of *u*_{Rt,B}(*T*_{stim}) decreases with increasing values of *u*_{Rt,B}(0) as a result of heightened GABA_{B}-mediated inhibition and leads to a concomitant decrease in the response of POm neurons to the stimulus. For the case of large values of *u*_{Rt,B}(0), the inhibitory feedback dominates and activity of POm thalamus is suppressed. The fixed point occurs when activity in POm thalamus is not fully suppressed. This point is stable below a critical value of GABA_{B}-mediated conductance, e.g., *g*_{POm;Rt,B} ≤ 3.6 for the parameters of Fig. 2*C*. Last, although the reduced model has multiple parameters, the value and stability of the fixed point depends only on the products *g*_{Rt,POm}^{2}*g*_{Pom,Rt,B}, *g*_{POm;Rt,B}*u*_{Rt,B}(0), and *g*_{Rt,POm}^{2}/*u*_{Rt,B}(0). Thus an increase in *g*_{Rt;POm} can be offset by a decrease in *g*_{POm;Rt,B} and a rescaling of the value of the activity at the fixed point.

The latency *t*_{0} is a monotonic function of *g*_{POm;Rt,B} up to the point of period doubling (Fig. 3*A*). The number of spikes per cycle *N*_{POm}^{spikes} thus decreases with increasing stimulus frequency (*Eq. 21*; Fig. 2*B*). For values of *g*_{POm;Rt,B} that permit period doubling, i.e., 3.6 < *g*_{POm;Rt,B} < 7.1 for the parameters of Figs. 2*C* and 3*A*, the latency alternates with a period of two stimulus cycles. More complex behavior emerges for still larger values of *g*_{POm;Rt,B}. In general, we can evaluate the stability of the solution with a single stimulus cycle for arbitrary values of stimulus frequency, *t*_{B}, and *g*_{POm;Rt,B} (Fig. 3*B*). This solution is stable for all frequencies for small *g*_{POm;Rt,B} and is unstable for stimulus frequencies above, roughly, *f*_{stim} = 2 Hz for large values of *g*_{POm;Rt,B}. The dependence of the latency *t*_{0} on the stimulus frequency *f*_{stim} and *g*_{POm;Rt,B} is shown in Fig. 3*C*. For values of *g*_{POm;Rt,B} that do not lead to period doubling, the latency increases with *f*_{stim} and can be even somewhat larger than *t*_{B}/2 but not close to *t*_{B}. Thus the minimal model captures the essential scaling of increased latency with increased stimulation frequency, albeit for a constrained set of parameter values.

##### POm-Rt CIRCUIT FOR RECTANGULAR STIMULUS: LACK OF LATENCY CODING WITH AN ABRUPT STIMULUS ONSET.

To demonstrate the effect of the stimulus shape on the dependence of latency on stimulation frequency, we study the POm-Rt circuit with a brain stem input with a rectangular shape (*Eq. 18*). The POm activity with a period of a single stimulus cycle is stable if *g*_{POm;Rt,B} is below a critical value, corresponding to a period-doubling bifurcation. Below this critical *g*_{POm;Rt,B} value, and even somewhat above it, the latency *t*_{0} is zero; an example with *f*_{stim} = 8 Hz is shown in Fig. 4*A*. At values of *g*_{POm;Rt,B} larger than the critical value, the latency *t*_{0} alternates between 0 and a positive value. The regime of stability of the fixed point with a period of a single stimulus cycle is qualitatively similar to the corresponding regime of stability for triangular stimulus (cf. Figs. 3*B* and 4*B*). In contrast to the case of the triangular stimulus, the latency *t*_{0} for the rectangular stimuli is zero in the parameter regime in which the state with a period of a single stimulus cycle is stable (Fig. 4*C*). Thus as a principle, a gradually increasing brain stem input is essential to code stimulus frequency by spike latency via GABA_{B}-mediated synaptic feedback from Rt to POm thalamus.

##### POM-RT-VPM CIRCUIT: VPM INPUT TO RT LEADS TO MORE ROBUST DYNAMICS.

The reduced architecture of a POm-Rt circuit (Fig. 1*B*) with only GABA_{B}-mediated inhibition and a ramp-like input accounts for both the frequency-dependent latency and spike rate that is observed in POm thalamus (Fig. 1*A*). However, too small a value of *g*_{POm;Rt,B} leads to only latencies that are significantly smaller than the stimulus duration *t*_{B}, while too large a value leads to period doubling (Fig. 3), a phenomenon not apparent in the data (Ahissar 1998) (appendix c). This behavior lies in the coupled activation of the Rt and POm nuclei. We consider how excitation of Rt by VPm thalamus, rather than solely through POm thalamus, can decouple the activity of the Rt and POm nuclei, lead to more robust dynamics, and enable large latencies up to *t*_{B}. We begin with an analytically tractable model that incorporates two additional features over the POm-Rt circuit: solely feedforward excitatory connections from VPm to RT thalamus, with conductance *g*_{Rt;VPm} (Fig. 1*C*), and brain stem input the shape of which is identical for the two nuclei but the amplitude of which for POm thalamus is taken to be weaker than that for VPm thalamus, that is (24) A full analysis of this model is given in appendix d.

To gain insight into the role of VPm input to Rt thalamus, and to examine latency as a function of the coupling between these nuclei *g*_{Rt;VPm}, we first consider the circuit with POm-to-Rt feedback excitation turned off, i.e., *g*_{Rt;POm} = 0. For this case, the steady state with the period of the stimulus is always stable, but the output of POm thalamus is silent above the value *g*_{Rt;VPm}*. For both triangular (*Eq. 17*) and rectangular stimuli (*Eq. 18*), the spiking activity decays gradually from a maximal value to zero; this is shown for *f*_{stim} = 8 Hz in the examples of Fig. 5, *A* and *B*. For the triangular stimulus, the latency increases gradually from 0 to *t*_{B} (Fig. 5*A*). In contrast, for the rectangular stimulus, the latency remains zero for an extended interval of values of *g*_{Rt;VPm} that satisfy *g*_{POm;Rt,B} *u*(0) ≤1; note that *u*(0) increases gradually with *g*_{Rt;VPm}. The latency then jumps steeply and reaches *t*_{0} = 1 at *g*_{Rt;VPm}*, for which *g*_{POm;Rt,B} *u*(0) = exp(*t*_{B}/τ_{B}) (Fig. 5*B*). Therefore for a rectangular stimulus, the latency increases with frequency only for a restricted range of values of *g*_{Rt;VPm} for which POm thalamic activity is very small. This implies that a gradual increase of the stimulus activity with time is necessary to code stimulus frequency by latency for this purely feedforward architecture.

We now investigate the circuit that includes the POm-to-Rt feedback excitation. As we have now shown that a gradually increasing brain stem stimulus is necessary for coding frequency by latency, we present the results only for the triangular stimulus. The value and stability of the fixed point were studied as a function of *g*_{Rt;VPm} and the GABA_{B} conductance *g*_{POm;Rt,B} (Fig. 6*A*, *1* and *2*). We found that period doubling is avoided for a sufficiently large value of *g*_{Rt;VPm}, i.e., *g*_{Rt;VPm} > 0.34 for the parameters of Fig. 6*A*, which corresponds to the case where Rt thalamus is driven primarily by VPm as opposed to POm thalamus. The values of *g*_{Rt;VPm} that are needed to generate a substantial frequency-dependent latency effect in POm thalamus are about an order of magnitude smaller than the concomitant values of *g*_{Rt;POm} (Fig. 6*B*, *1* and *2*; *f*_{stim} = 8 Hz). For example, the ratio *t*_{0}/*t*_{B} = 0.75 is achieved for *g*_{Rt;VPm} = 0.6 when *g*_{Rt;POm} = 0, as opposed to *g*_{Rt;POm} = 5 when *g*_{Rt;VPm} = 0. The cost of this robust solution is that the neurons in POm thalamus can be completely silent for *g*_{Rt;VPm} values above *g*_{Rt;VPm}* (Fig. 6, *A1* and *B1;* appendix d, *Eq.* D*6*). Thus feedforward connections from VPm to Rt thalamus obviate period doubling and lead to a robust frequency-dependent latency, regardless of whether feedback connections from POm to Rt thalamus are present.

##### POM-RT-VPM CIRCUIT: RT INHIBITORY FEEDBACK TO POM SHOULD BE STRONGER THAN THAT TO VPM.

The incorporation of GABA_{B}-mediated inhibition from Rt to VPm thalamus completes the full pattern of known connectivity (Figs 1, *A* and *D*). In principle, this can lead to a frequency-dependent latency in the activity of VPm as well as POm neurons. The relative values of the latency are determined by the ratio between the average level of GABA_{B}-mediated inhibition and the average strength of the brain stem input (*Eq. 24*), i.e., *g*_{VPm;Rt,B} versus α^{-1}*g*_{POm;Rt,B} for VPm versus POm thalamus, respectively. The former value should be the smaller of the two to ensure that the latency in VPm thalamus is essentially constant.

### Numerical results for the full POm-Rt-VPm circuit

The analytical theory suggests an explanation for the latency coding. Still, several questions still remain open. First, are the analytical results valid despite the approximations that were made? Second, what is the effect of the dual-sloped rising phase of the brain stem stimulus, as observed experimentally (Fig. 1*A*), as opposed to just a triangular shape? Third, can we explain the detailed shape of the PSTHs recorded in the POm and the VPm nuclei? Specifically, why does the POm activity decay rapidly with time after the initial onset at low stimulus frequency, whereas such decay is seen neither in the POm in response to high-frequency stimulation nor in the VPm?

To answer the preceding questions and to draw a close comparison between responses calculated for the model and those observed in experiments, we have carried out numerical simulations of the full rate model (appendix a) with the architecture of Fig. 1*D*. We incorporate inhibitory feedback that satisfies *g*_{POm;Rt,B} > *g*_{VPm;Rt,B}; adaptation in the activity of POm thalamus, and weaker adaptation in VPm thalamus, to account for the reduction in the population of neurons that are available to spike after each stimulus (Eggert and van Hemmen 2001 2000); and a dual-sloped rising phase for the stimulus (Fig. 7, *A* and *B*).

##### COMPARISON OF THEORY WITH EXPERIMENT: WAVEFORMS.

There is little activation of GABA_{B}-mediated inhibition at low frequencies, i.e., *f*_{stim} = 2 Hz (Fig. 7*C*). For such low frequencies, the activity in POm thalamus consists of an initial peak that rapidly attenuates (Fig. 7*E*) primarily as a consequence of intrinsic adaptation, although GABA_{A}-mediated inhibition (Fig. 7*D*) also contributes. The activity in VPm thalamus also shows an initial peak that is followed by an attenuated response (Fig. 7*F*). However, there is now a second peak that results from the waning attenuation and the continued rise of the stimulus. The second peak occurs in VPm thalamus, but not POm thalamus, as a consequence of the relatively weaker adaptation.

GABA_{B}-mediated inhibition by Rt thalamus becomes strongly activated with increased values of the stimulation frequency above *f*_{s}_{tim} = 2 Hz (Fig. 7*C*). The excitatory stimulus minus the GABA_{B}-mediated inhibition is positive only near the peak of the stimulus, where the stimulus rises gradually, as in the case of a triangular stimulus. As a result, the activity in POm thalamus is strongly suppressed just after the onset of the stimulus (Fig. 7*E*). With increasing time, the brain stem input overcomes the GABA_{B}-mediated inhibition and the activity in POm rises until the stimulus ends. For times ∼50–90 ms after the stimulus onset, the GABA_{B} synaptic activation is balanced by the decreasing adaptation and the activity in POm thalamus is essentially independent of *f*_{stim} (Fig. 7*E*). These dynamics result in a calculated instantaneous spike rate whose latency to peak and overall amplitude compares favorably with those in the experimental data (cf. Figs. 7*E* with 1*A*).

With regard to VPm thalamus, the slow GABA_{B}-mediated inhibition reduces the overall activity with increasing stimulus frequency (Fig. 7*F*). Yet, as a result of a larger input from the brain stem (α = 0.6 in *Eq. 13*) and the weaker GABA_{B}-mediated inhibition compared with the case for POm thalamus, the excitatory stimulus minus the GABA_{B}-mediated inhibition is positive everywhere except just after the onset of the stimulus. The reduction of activity by GABA_{B}-mediated inhibition decreases the firing activity just following the stimulus onset but is insufficiently strong to affect the latency of the response. Thus there is negligible change in the latency to peak, with increasing stimulus frequency. As in the case of POm thalamus, the time dependence of the instantaneous spike rates in our simulations compares favorably with those in the data (cf. Figs. 7*F* with 1*A*).

##### COMPARISON OF THEORY WITH EXPERIMENT: LATENCY AND SPIKE COUNT.

The calculated and observed latencies were expressed in terms of *t*_{1/2}, which is the time for the instantaneous rates of spiking, *M*_{POm}(*t*) and *M*_{VPm}(*t*), to reach one-half of their maximum value. The values of *t*_{1/2} computed from the model well approximate the observations for both the POm and VPm nuclei (Fig. 7, *G* and *H*). Further, the calculated increase in the value of *t*_{1/2} with increasing stimulus frequency for POm thalamus is a robust, monotonic function of *g*_{POm;Rt,B} up to the frequency where POm neurons are silenced by inhibition (Fig. 7*I*). In contrast to the close match of theory and observation for the latency effect, the calculated values of the number of spikes per cycle in POm thalamus decreases more sharply with frequency than in the observations (Figs. 7*G*). This deviation results from adaptation at low frequencies of stimulation that is relatively weak in the model. Nonetheless, the overall match between calculated and observed values (Fig. 7, *E*–*H*) lends support to the idea that the heightened frequency-dependent latency and reduced spike-count in POm versus VPm thalamus results from a stronger GABA_{B}-mediated inhibitory feedback, a weaker brain stem input for neurons in POm thalamus (*Eq. 13*), and the double-ramp, monotonically increasing brain stem stimulus.

##### EFFECT OF STIMULUS DURATION AND SHAPE.

We seek to account for the effect of stimulus width on the latency effect in POm thalamus and further predict the nature of the spiking response in the POm and VPm nuclei for arbitrary stimuli. Intuitively, the strongly facilitating nature of the GABA_{B}-mediated currents should make the rate of spiking sensitive to the duration of the stimulus.

We first consider variations in the total duration of the stimulus (Fig. 8*A*). The period of the initial rise was kept fixed but the duration of the second rising phase and the falling phase were co-varied. The duration was parameterized by a multiplicative factor, denoted β, where β = 1 is the value for a stimulus with a total rising phase of 50 ms (Fig. 8*A1*). The change in time course of *M*_{POm}(*t*) and *M*_{VPm}(*t*) as a function of frequency is relatively weak for the case of a brief stimulus, i.e., β = 0.4 for a rise-time of 23 ms (Fig. 8*A*, *2* and *3*), compared with a long stimulus (Fig. 7, *E* and *F*). The calculated time course compares favorably with those recorded with a 20 ms stimulus (Fig. 7 in Sosnik et al. 2001). In general, the range of the frequency-dependent latency for spiking in POm thalamus is predicted to increase as the duration of the stimulus increases toward its full width for β = 1 (Fig. 8*A4*).

We next consider changes to the shape of the dual-sloped rising ramp of the stimulus (Fig. 8*B*). The slope of the first ramp and the total duration of the stimulus remained fixed, while the slope of the second ramp was systematically varied by changing the duration of the first ramp (Fig. 8*B1*). For the particular case of a flat-top stimulus, i.e., a duration of the first ramp of 9.4 ms, the instantaneous spike rate for both POm and VPm nuclei exhibit latencies that do not appreciably increase with increasing frequency (Fig. 8*B*, *2* and *3*). In general, the latency for spiking in POm thalamus increases with frequency only if the slope of the second ramp is positive (Fig. 8*B4*) and if the GABA_{B}-mediated inhibition is sufficiently strong, so that the difference between the stimulus input and the inhibitory currents increase gradually over time. Our simulation results, which are consistent with our analytical results, provide an explanation for the absence of a latency effect for mechanical ramp-and-hold versus air-puff stimulation of the vibrissae (Diamond et al. 1992; Sosnik et al. 2001) as the former method typically involves a much more rapid initial movement.

## DISCUSSION

We have shown that GABA_{B}-mediated inhibition can play a critical computational role in neural coding of stimulus frequency by modulating the latency of the spike response in POm thalamus. Our analysis suggests that feedback among only two nuclei, POm and Rt thalamus, is sufficient to generate this effect provided that the stimulus rises gradually with time (Figs. 2–4). A full architecture that further incorporates reciprocal connections with VPm thalamus leads to particularly robust dynamics (Figs. 5 and 6) and offers contrast between the spiking properties in POm versus VPm thalamus (Fig. 7). We calculated the output from the full model under the assumption that Rt thalamus inhibits POm thalamus more strongly than it inhibits VP_{m} thalamus and that the brain stem input to POm thalamus is weaker than that to VPm thalamus. The corresponding values for the latency of spiking and the number of spikes per cycle are in good agreement with the observed values for both thalamic areas over a large range of stimulus frequencies (Fig. 7, *G* and *H*).

### Comparison with experimental results

Analysis of our model reveals that the latency increases with frequency in POm thalamus but not in the VPm thalamus if the brain stem input to the VPm is larger than the brain stem input to the POm and the GABA_{B}–mediated inhibition is stronger in the POm than in the VPm. The first condition is clearly seen in in vivo experiments because the firing rate of VPm neurons in low stimulus frequencies is considerably higher in the VPm than in the POm (Diamond 1995; Diamond et al. 1992; Sosnik et al. 2001). Furthermore, in recent in vitro experiments, C. E. Landisman and B. W. Connors (unpublished data) found that VPm neurons exhibit higher maximal firing rates from POm neurons when stimulated with injected current. These researchers also found that VPm cells were 2.5 times less likely than POm cells to have a GABA_{B} receptor-like component as part of their inhibitory responses. Hence these recent experimental results support our view of the role of GABA_{B}–mediated inhibition in coding frequency by latency in the POm but not in the VPm.

We find that strong feedback connections between POm and Rt thalamus destabilize the steady-state activity with the period of the stimulus via a period doubling bifurcation. Analysis of the experimental data (Sosnik et al. 2001) (appendix c) did not reveal any indication for activity with double period. Recent observations (Webber and Stanley 2004) have revealed period doubling in cortical activity in response to stimuli with duty cycles of 0.5 and frequencies between 4 and 8 Hz. This behavior may reflect similar behavior in the thalamic input to cortex. The stimuli used by those authors are more prolonged than the stimuli analyzed here and may generate a strong GABA_{B}-mediated inhibitory response.

The GABA_{B} kinetics used in the model (*Eqs. 9* and *10*) are based on the results from dual recordings (Kim et al. 1997). These observations show that a single presynaptic neuron should fire a high-frequency burst of action potentials to evoke a substantial GABA_{B}-mediated response from the postsynaptic neuron. Indeed, Rt cells fire at a high rate in comparison with thalamic relay nuclei (Hartings et al. 2000). Furthermore because GABA_{B}–mediated synapses are metabotropic and involve the cooperative activation of G proteins (Destexhe and Sejnowski 1995; Golomb et al. 1996), near-synchronous spiking of several presynaptic cells may act cooperatively to facilitate the stronger GABA_{B}-mediated response.

### Consequences of approximations in our model

Single-cell properties may well affect the latency and adaptation effects in the POm nucleus. Although there is little data on the cellular properties of neurons in POm thalamus, the lateral posterior (LP) nucleus is an analogous nucleus in the visual stream (Peters 1985) and the nonlemniscal nuclei are analogous nuclei in the auditory stream (Yu et al. 2004). In vitro intracellular recordings from neurons in the LP nucleus revealed significant potassium A-type and calcium T-type currents (Li et al. 2003). The A current can further delay the occurrence of the first spike in a stimulus cycle, whereas the combination of a T current and inhibitory input may contribute to the production of a burst of spikes followed by an afterhyperpolarization that terminates neuronal activity (Sherman 2001; Steriade et al. 1993). These cellular properties are likely to affect the details of our analysis of the full model (Fig. 7), but not our conclusions regarding the latency effect (Figs. 2–6).

The difference in latency coding between the VPm and POm thalamus is also reflected by neurons in cortical layers that receive axonal projection from these nuclei, i.e., layers IV and Va, respectively (Ahissar et al. 2000, 2001; Ahrens et al. 2002). Cortical feedback will act to supplement the brain stem input to thalamic relay cells (Deschenes et al. 1998; Diamond et al. 1992), and especially to excite Rt thalamus (Golshani et al. 2001). Recalling that thalamocortical (Gil and Amitai 1996) and even corticothalamic (Swadlow 1994) synaptic connections are much faster than the decay time of GABA_{B}-mediated synapses, we anticipate that the introduction of cortical effects into the present thalamic model can, to a first approximation, be accomplished by a variation of the strength of the synaptic connections within the present architecture (Fig. 1*D*). As a test of this conjecture, we examined a circuit that included cortical feedback such that the Rt thalamus received all its input from excitatory cortical cells (appendix e). We found that the dynamics of this circuit was similar to that of the solely thalamic circuit (Fig. 7) in support of our reductionism. Other factors, such as GABA_{B}-mediated synapses within the neocortex (Garabedian et al. 2003), may influence the latency effect in POm thalamus at high frequencies of stimulation.

The POm nucleus is innervated by GABAergic neurons from zona incerta (ZI) (Bartho et al. 2002; Lavallee et al. 2005; Power et al. 1999; Trageser and Keller 2004) and from the anterior pretectal (APT) nucleus (Bokor et al. 2005). The ZI and the APT nuclei are innervated by trigeminal inputs, and the ZI is also innervated by layer 5 cortical projections. It is not clear whether ZI or APT can play the role of Rt thalamus in shaping the latency since there is currently no experimental evidence for a GABA_{B}-mediated connection from these nuclei to thalamic nuclei. From the perspective of temporal dynamics, the GABA_{A}-mediated inhibition from neurons in ZI and in APT nucleus (Bokor et al. 2005) plays the same role as the GABA_{A}-mediated inhibition from Rt neurons, because the brief delay or the brief advance of inhibition is negligible in comparison with the relevant time scale of order of tens of milliseconds.

The analytically solvable model was developed from the full model (*Eqs. 1*–*13*) using several approximations: the GABA_{B}-mediated decay-time, τ_{B}, is larger than all the other synaptic time constants in the system; adaptation is ignored; GABA_{A}-mediated inhibition is ignored; and the thresholds, θ_{i}, are set to 0. We numerically simulated the full model to show that the results of an analytically solvable simplified model did not depend on the preceding approximations. Specifically, the circuit produces an increase in latency with increasing frequency when GABA_{A}-mediated synapses and adaptation are included. The simulations show that GABA_{A}-mediated inhibition and adaptation bring the temporal structure of the solution of the model closer to the observed structure.

### Adaptation and depression

Adaptation is incorporated in the full model to account for the fast decay of spiking activity in POm but not VPm thalamus at low stimulus frequencies. The adaptation process in the model (*Eqs. 2*, *3*, *11*, and *12*) is multiplicative. Qualitatively, this process takes into account the fact that a cell that has just fired cannot fire anymore until the adaptation effects vanish. Biophysically, this process may be based on the activation of the T-type Ca^{2+} current that leads to a burst of one or more spikes followed by a prolonged refractory period. Interestingly, POm neurons are more prone to burst mode activity than VPm neurons (Ramcharan et al. 2004), consistent with the stronger adaptation they exhibit in vivo. Adaptation may also result from the activation of slow K^{+} currents. In previous rate models (Hansel and Sompolinsky 1998; Tabak et al. 2000), adaptation was modeled as a subtractive process. However, we could not mimic the strong adaptation at low *f*_{stim} with such a process.

In principle, the role of GABA_{B}-mediated synapses proposed here can be replaced by synaptic depression of the brain-stem-to-POm excitatory synaptic connections or by intrinsic adaptation with the proper time scale, mediated, for example, by slow Ca^{2+}-activated K^{+} currents or other long-lasting cellular adaptation processes. We are not aware of any experiment in which the dynamic properties of the brain-stem-to-POm synapses were quantified. A depression-based mechanism will likely require that the depression of the brain-stem-to-POm synapses be considerably stronger than the depression of the brain-stem-to-VPm synapses. A mechanism based on intrinsic adaptation may require a second adaptation process in addition to the one that is responsible for the gradual decline of POm activity at low frequency. Because mechanisms based on depression or adaptation do not include an interaction between the POm and VPm nuclei, it is expected that the steady state with the period of the stimulus will undergo a period doubling bifurcation if the latency at high frequencies is large enough.

### Predictions from the model

Our analysis yields several predictions for future experiments. First, the complete blockade of GABA_{B}-mediated synaptic transmission in POm thalamus will eliminate the increase of latency with increasing stimulus frequency (Fig. 7*I*). Second, barring adaptation effects, the blockade of GABA_{B}-mediated synaptic transmission in either the POm or VPm nuclei should eliminate the decrease in spike rate with increasing frequency. In partial support of this claim, Castro-Alamancos (2002a) blocked both GABA_{A}- and GABA_{B}-mediated inhibition in VPm thalamus and found that the spiking activity no longer decreases with the frequency of stimulation. Third, partial blockage of GABA_{B}-mediated inhibition will preferentially eliminate the increase in latency but not the decrease in spiking with increasing frequency. Fourth, we argued that activation of Rt thalamus via excitatory input from VPm thalamus, which acts as feedforward inhibition from the VPm to the POm nucleus, is essential for a robust latency effect. Thus blockage of activity in VPm thalamus should suppress the latency effect in POm thalamus (Fig. 6). Last, our analysis emphasized the importance of the shape of the spiking pattern in brain stem, which follows the time course of vibrissa movement (Deschenes et al. 2003; Jones et al. 2004; Lichtenstein et al. 1990), in the latency effect. We expect that the latency effect will be absent if deflection of the vibrissae occurs too quickly (Figs. 5*B* and 8*B*) or is too brief (Fig. 8*A*).

### Comparison with a phase-locked loop algorithm

The dependence of the latency in POm thalamus on the stimulus frequency has been shown to be qualitatively consistent with the behavior of a phase-locked loop (PLL) algorithm (Ahissar and Kleinfeld 2003; Ahissar et al. 2000). The implementation of the algorithm involves negative feedback to neurons in POm thalamus from neuronal oscillators in cortex or thalamus. In principle, a thalamocortical PLL circuit is similar to our reduced POm-Rt circuit (Fig. 1*B*) but with an oscillator replacing Rt. In both models, the latency of the response in POm thalamus is determined by the interaction between the inhibitory feedback signal and the brain stem input. However, the dependency of the latency on the frequency stems from different sources. In the PLL model, it arises from a phase correction of the oscillator in Rt thalamus, whereas in the present neuronal model, it arises from a gradual increase in the amplitude of the GABA_{B}-mediated inhibition. Thus for example, a rectangular brain stem stimulus is expected to yield latency coding within a PLL circuit but not with a GABA_{B}-mediated mechanism (Fig. 5*B*).

### Computational relevance of Rt thalamus

Experimental (Steriade et al. 1993; Von Krosigk et al. 1993) and theoretical (Golomb et al. 1996) studies have demonstrated the importance of the Rt nucleus in the generation of brain oscillations, such as spindles. The Rt nucleus is also a locus for epileptic seizures (McCormick and Contreras 2001). Here we suggest a computational role for the Rt nucleus in somatosensory information processing: it employs GABA_{B}-mediated inhibition to code stimulus frequency by the latency of the spike response in POm thalamus. The Rt nucleus can modulate sensory inputs using GABA_{B}-mediated synapses if the stimulus is strongly temporally-modulated in the time scale of the decay of GABA_{B}-mediated inhibition, namely the stimulus frequency should be ∼5–10 Hz. A more complex role that involves GABA_{B}-mediated inhibition, such as the detection of a change in timing or amplitude in an otherwise periodic sensory stream (Mehta and Kleinfeld 2004), remains to be determined.

## APPENDIX A: PARAMETERS OF BRAIN STEM INPUT

The brain stem activity in response to the first stimulus in a chain is given by (A1)

For *t* > *T*_{stim}, the brain stem input is *I*_{VPm}(*t*) = *I*_{VPm}(*t* − *kT*_{stim}), where *k* is the largest integer smaller than *t*/*T*_{stim}. The brain stem input to POm thalamus is proportional to the brain stem input to the VPm, and it is delayed to allow larger latency in POm even for low stimulus frequencies (A2) α = 0.6, *t*_{POm,delay} = 7 ms.

The values used for Fig. 7 are: *q*_{0} = 0, *q*_{1} = 6 ms, *I*_{0} = *I*_{1} = 0, *q*_{2} = 11 ms, *I*_{2} = 0.8, *q*_{3} = 56 ms, *I*_{3} = 1.5, *q*_{4} = 96 ms, *I*_{4} = 0. The starting times of the first ramp for the inputs to the VPm and POm nuclei, *q*_{1} and *q*_{1}+ *t*_{POm,delay}, were chosen such that the latencies *t*_{0.5} in the two nuclei, obtained from the model at low stimulus frequencies, fit the latencies observed experimentally.

The following parameters are used for Fig. 8.

*A*) The width of the second rising phase and the decaying phase together is β times the reference values of these phases, respectively. The parameters are *q*_{0} = 0, *q*_{1} = 6 ms, *q*_{2} = 11 ms, *I*_{0} = *I*_{1} = 0, *I*_{2} = 0.8, *q*_{3} = *q*_{2} + β × 45 ms, *I*_{3} = *I*_{2} + β × 0.7, *q*_{4} = *q*_{3} + β × 40 ms, *I*_{4} = 0.

*B*) The slope of the first rising phase is d*I/*d*t* = 0.16 ms^{−1}; *I*_{2} = (d*I/*d*t*) (*q*_{2} − *q*_{1}). The other parameters are: *q*_{0} = 0, *q*_{1} = 6 ms, *I*_{0} = *I*_{1} = 0, *q*_{3} = 56 ms, *I*_{3} = 1.5, *q*_{4} = 96 ms, *I*_{4} = 0. The values of *q*_{2} are 15.4 ms for Fig. 8, *B1* (solid line), *B2*, *B3*; 11 ms for *B1* (dotted line); and 19.8 ms for *B1* (dashed line).

## APPENDIX B: ANALYSIS OF THE REDUCED POM-RT MODEL

We consider the POm-Rt architecture (Fig. 1*B*), namely *g*_{Rt;VPm} = 0. *Eq. 16* becomes (B1)

Our goal is to find *u*_{Rt,B}(*T*_{stim}), the GABA_{B}-mediated synaptic activation after one period, knowing *u*_{Rt,B}(0). The variable *u*_{Rt,B} is affected by the stimulus at time *t* if *I*_{POm}(*t* − *t*_{B}) − *g*_{POm;Rt,B} *u*_{Rt,B}(*t* − *t*_{B}) > 0 (*Eq.* B*1*). We denote the onset time of activity of POm neurons by *t*_{0}. From now on, the analysis depends on the stimulus shape.

### Triangular stimulus (Eq. 17)

The POm onset time *t*_{0} is determined by the equation (Fig. 2*A*) (B2)

This equation is valid only if *t*_{0} < *t*_{B}, namely *u*_{Rt,B} (0) < u_{Rt,B;max} = . If *u*_{Rt,B}(0) > *u*_{Rt,B;max}, *u*_{Rt,B} decays exponentially during the whole period (B3)

To calculate *u*_{Rt,B}(*T*_{stim}) for smaller *u*_{Rt,B}(0), we note that *u*_{Rt,B}(*t*_{B} + *t*_{0}) = , and integrate *Eq. B1* to obtain *u*_{Rt,B}(2t_{B}) (B4) Using the equation (B5) we obtain (B6)

*Equations* B*3* and B*6* define a map from the value of the GABA_{B}-mediated conductance variable *u*_{Rt,B}(0) at the beginning of one stimulus to the value of *u*_{Rt,B}(*T*_{stim}) at the beginning of the next stimulus. The intersection between the map (*Eqs.* B*3* and B*6*) and the diagonal line *u*_{Rt,B}(*T*_{stim}) = *u*_{Rt,B}(0) is the fixed point of the map, denoted by . This fixed point corresponds to *u*_{Rt,B} being periodic with a period of a single stimulus cycle. The fixed point is stable if (Berge et al. 1984) (B7)

The fixed point undergoes a period doubling (PD) bifurcation if the derivative in *Eq.* B*7* equals −1. The value *u*_{Rt,B}(*T*_{stim}) depends on *u*_{Rt,B}(0) explicitly and also implicitly through the latency *t*_{0} (*Eq.* B*2*). Therefore (B8)

By differentiating *Eq.* B*4* with respect to *t*_{0} and using *Eqs.* B*5* and B*2,* one sees that ∂*u*_{Rt,B}/∂*t*_{0} = 0. Therefore d*u*_{Rt,B}(*T*_{stim})/d*u*_{Rt,B}(0) = ∂*u*_{Rt,B}(*T*_{stim})/∂*u*_{Rt,B}(0), and the PD bifurcation point is given by (B9) To compute the PD bifurcation line, we follow the solutions of *Eqs.* B*6*, B*2*, and B*9* using the program XPPAUT (Ermentrout 2002).

The quantity *N*_{POm}^{spikes}, proportional to the number of spikes fired by POm neurons (*Eq. 19*, Fig. 2*A*) can be calculated exactly. We find (B10)

### Rectangular stimulus (Eq. 18)

According to *Eqs.* B*1* and *18*, three behavioral regimes are determined by the initial value *g*_{POm;Rt,B} *u*_{Rt,B}(0).

First, if *g*_{POm;Rt,B} *u*_{Rt,B}(0) ≤ 1, the POm onset time is *t*_{0} = 0, and (B11)

Second, if 1 < *g*_{POm;Rt,B} *u*_{Rt,B}(0) < , the POm onset time is *t*_{0} = τ_{B} log[*g*_{POm;Rt,B} *u*_{Rt,B}(0)], and (B12)

Third, if *g*_{POm;Rt,B} *u*_{Rt,B}(0) ≥ , *u*_{Rt,B} decays exponentially during the whole period according to *Eq.* B*3*.

Numerical simulations show that the fixed point of the map *u*_{Rt,B}(*T*_{stim}) versus *u*_{Rt,B}(0) is given by the intersection of the line *u*_{Rt,B}(*T*_{stim}) = *u*_{Rt,B}(0) with *Eq.* B*11*. The PD bifurcation point is given by (B13)

## APPENDIX C: TESTING FOR A PERIOD-2 BEHAVIOR IN THE DATA

Do experimental recording from POm neurons exhibit doublet behavior? To answer this question, we re-examined 18 multiunit recordings from the POm thalamus (Sosnik et al. 2001). For comparison, we also re-examined 23 multiunit recordings from the VPm thalamus (Sosnik et al. 2001). If the multiunit recording exhibits period-2 behavior, the average absolute value of the difference in the number of spikes between the response to a first stimulus and the response to the subsequent stimulus should be large in comparison to the absolute value of the difference in the number of spikes between the response to a first stimulus and the response to the second subsequent stimulus. To quantify this behavior, we denote the number of spikes fired by the population in response to the *i*th stimulus by *r*_{i}. In practice, each train of stimuli lasts for 3 s, and we drop the response for the first second to avoid transient effects. In the remaining 2 s, there are *N*_{T} stimuli. For example, if *f*_{stim} = 8 Hz, *N*_{T} = 16. We define the value *S*_{i} as (Fig. C1*A*) (C1)

We define *Y* as (C2)

The trial-average (over 24 trials) and the SD of *Y* for a specific multiunit recording and a specific period are denoted by 〈*Y*〉 and σ_{Y}, respectively. This analysis is an alternative to spectral analysis.

The *Z* score of the distribution is defined as (C3)

The value of *Z* should be large compared with 1 if the unit exhibits doublet behavior. Figure C1 *B*–*G,*, shows the distributions *P*(*Z*), computed for the POm and VPm populations of units and for *f*_{stim} = 5, 8, and 11 Hz. The distributions *P*(*Z*) for both POm and VPm units for all those frequencies are centered around zero. All the values of *Z* were <1, except for one case of a POm unit and *f*_{stim} = 8 Hz, where *Z* = 1.2. We conclude that there is no statistically significant indication for doublet behavior in the data.

## APPENDIX D: ANALYSIS OF THE VPm-POm-Rt CIRCUIT WITH FEEDFORWARD VPm-TO-Rt CONNECTIONS

### Triangular stimulus (Eq. 17)

For the VPm-POm-Rt circuit without Rt-to-VPm feedback connections (Fig. 1*C*), *Eq. 16* for the GABA_{B}-mediated synaptic activation from the Rt becomes (D1)

For a triangular stimulus (*Eq. 17*), the POm term in *Eq.* D*1* is nonzero for *t*_{0} < *t* − *t*_{B} < *t*_{B}, where *t*_{0} is given by *Eq.* B*2*. Therefore (D2) and (D3) Using *Eqs.* D*2,* D*3,* and B*5* we obtain (D4)

The stability of the fixed point *u*_{Rt,B}^{FP} = *u*_{Rt,B}(*T*_{stim}) = *u*_{Rt,B}(0) is given by *Eqs.* B*7* and B*8*. Using *Eqs.* D*2* and D*3,* one shows that ∂*u*_{Rt,B}(*T*_{stim})/∂*t*_{0} = 0. Therefore *du*_{Rt,B}(*T*_{stim})/*du*_{Rt,B}(0) = ∂*u*_{Rt,B}(*T*_{stim})/∂*u*_{Rt,B}(0), and the PD bifurcation point is given by (D5)

When the value *g*_{Rt;VPm} reaches the value *g*_{Rt;VPm}*, *t*_{0} = *t*_{B}, the POm thalamus is silent, and its contribution to *Eq.* D*3* is zero. Therefore *g*_{Rt,VPm}* does not depend on *g*_{Rt;POm}, and is given by (Figs. 5 and 6) (D6)

### Rectangular stimulus (Eq. 18).

For simplicity, we limit ourselves to a purely feedforward VPm-POm-Rt circuit: *g*_{VPm,Rt;B} = *g*_{Rt,POm} = 0 (Fig. 1*C*), Using *Eqs. D1* and *18,* we obtain the fixed point value (D7)

For *g*_{POm,Rt,B}*u*_{Rt,B}^{FP} ≤ 1, *t*_{0} = 0 and (D8)

For 1 < *g*_{POm,Rt,B}*u*_{Rt,B}^{FP} ≤ and (D9)

## APPENDIX E: CORTICAL FEEDBACK

The thalamus and the cortex are mutually coupled. In this work, the thalamo-thalamo connections subsume part of the feedback role played by the thalmocortical connections. This simplification enables us to analyze the model and define the components of the circuit that lead to the experimentally observed relationship between stimulus frequency and latency. We now consider the additional of a simplified form of the thalamocortical circuit (Fig E1*A*). The excitatory (CxE) and inhibitory (CxI) cortical neurons are excited by the thalamic relay nuclei and are mutually coupled. The Rt is excited only by the excitatory cortical neurons. The system dynamics is described by *Eqs. 1*–*3* and *5*–*13*, while *Eq. 4* is replaced by (E1) The instantaneous spiking rates of the cortical populations are (E2) (E3) The relationships between the variables *u*_{CxE} and *u*_{CxI} and the rates *M*_{CxE} and *M*_{CxI} are given by equations similar to *Eqs. 5*–*8*.

The parameters that are different from those of Fig. 7 are: *g*_{Rt;VPm} = 0, *g*_{Rt;POm} = 0. Additional parameters for conductances, delays, and thresholds are: *g*_{Rt;CxE} = 1.5, *g*_{CxE;VPm} = 1, *g*_{CxE;POm} = 1, *g*_{CxE;CxI,A} = 0.2, *g*_{CxI;VPm} = 1, *g*_{CxI;POm} = 1, *g*_{CxI;CxE} = 0.5, *t*_{Rt;CxE,G} = 2 ms, *t*_{CxE;VPm,G} = 2 ms, *t*_{CxE;POm,G} = 2 ms, *t*_{CxE;CxI,A} = 1 ms, *t*_{CxI;VPm,G} = 2, *t*_{CxI;POm,G} = 2, *t*_{CxI;CxE,G} = 1 ms, θ_{CxE} = 0, θ_{CxI} = 0.

The dependencies of the waveforms of the GABA_{A} and GABA_{B} inhibition and the VPm and POm spiking activity on *f*_{stim} (Fig. E1, *B*–*E*) are very similar to those we computed for the thalamic-only circuit (Fig. 7, *C*–*F*). Similarly, the dependencies of the latencies *t*_{1/2} and the number of spikes per cycle in the POm and the VPM on *f*_{stim} are also very similar in the two circuits (Figs. E1, *F* and *G*, and 7, *G* and *H*). The reason behind this resemblance is that the dependence of the system dynamics on *f*_{stim} is governed by the slow Rt-to-POm and Rt-to-VPm GABA_{B} dynamics. These dynamics are dependent on the Rt activity during the cycle, and it does not matter whether the Rt is activated by the cortical feedback or directly by the thalamic relay nuclei.

## GRANTS

This work was supported by Israel Science Foundation Grants 311/04 to D. Golomb and 377/02 to E. Ahissar, United States–Israel Binational Science Foundation Grant 2003299 to E. Ahissar and D. Kleinfeld, National Science Foundation (Center Grant PHY0216576), National Institute of Mental Health Grant MH-59867 to D. Kleinfled, a grant from the Abramson Family Foundation and the Nella and Leon Benoziyo Center to E. Ahissar, and Human Frontiers Science Program Grant RGP43/2004 to E. Ahissar and D. Kleinfeld.

## NOTE IN PROOF

In recent work presented at the Society for Neuroscience 2005 Annual Meeting (Washington DC), Trageser, Keller and colleagues (Abstract 173.1) carried out experiments in which the steady-state responses of neurons in POm thalamus were investigated as a function of the frequency of vibrissa stimulation. These investigators report that the probability of spiking by POm neurons decreases with increased stimulation frequency, such that many neurons failed to respond when stimulated above 11 Hz. These preliminary findings are consistent with our prediction (Figs. 6 and 7*I*) that POm thalamus can be silent at high-stimulus frequencies.

## Acknowledgments

We thank B. Connors and C. Landisman for discussions on their unpublished data and D. Hill for comments on an early version of the manuscript.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 by the American Physiological Society