Performance Limitations of Relay Neurons

Relay cells are prevalent throughout sensory systems and receive two types of inputs: driving and modulating. The driving input contains receptive field properties that must be transmitted while the modulating input alters the specifics of transmission. For example, the visual thalamus contains relay neurons that receive driving inputs from the retina that encode a visual image, and modulating inputs from reticular activating system and layer 6 of visual cortex that control what aspects of the image will be relayed back to visual cortex for perception. What gets relayed depends on several factors such as attentional demands and a subject's goals. In this paper, we analyze a biophysical based model of a relay cell and use systems theoretic tools to construct analytic bounds on how well the cell transmits a driving input as a function of the neuron's electrophysiological properties, the modulating input, and the driving signal parameters. We assume that the modulating input belongs to a class of sinusoidal signals and that the driving input is an irregular train of pulses with inter-pulse intervals obeying an exponential distribution. Our analysis applies to any order model as long as the neuron does not spike without a driving input pulse and exhibits a refractory period. Our bounds on relay reliability contain performance obtained through simulation of a second and third order model, and suggest, for instance, that if the frequency of the modulating input increases or the DC offset decreases, then relay increases. Our analysis also shows, for the first time, how the biophysical properties of the neuron (e.g. ion channel dynamics) define the oscillatory patterns needed in the modulating input for appropriately timed relay of sensory information. In our discussion, we describe how our bounds predict experimentally observed neural activity in the basal ganglia in (i) health, (ii) in Parkinson's disease (PD), and (iii) in PD during therapeutic deep brain stimulation. Our bounds also predict different rhythms that emerge in the lateral geniculate nucleus in the thalamus during different attentional states.


Introduction
Relay neurons are found in various brain nuclei including the thalamus [1][2][3]. Experiments have suggested that the inputs to a thalamic relay neuron can be divided into two categories: driving and modulating. The driving input typically contains sensory information (e.g visual, motor) and the modulating input controls relay of this sensory information back to cortex [4]. The driving input is made up of a few synapses on the proximal dendrites whereas the modulating input comprises all other synapses [5,6] (see Figure 1 A).
For example, the lateral geniculate nucleus (LGN) in the thalamus receives the driving input from the retina and projects to the primary visual cortex. The modulating input comprises descending inputs from layer 6 of the visual cortex and ascending inputs from the brain stem. The function of the LGN is to selectively relay sensory information from the retina subject to attentional needs [4,7]. It has been observed that during different attentional needs (which translate into different relay demands), local field potentials (LFPs) in the LGN have a concentration of power in different frequency bands (D?c) [8,9]. LFPs may be reflected in the modulating input because they are believed to arise from ensemble synaptic activity [10]. This would then suggest that one mechanism that controls relay in the LGN cell is the frequency of the modulating input.
Similarly, relay neurons in the motor thalamus receive driving inputs from sensorimotor cortex, and modulating inputs from the basal ganglia (BG), specifically the Globus Pallidus internal segment (GPi) [4,11]. The driving input contains information about the actual movement via proprioception, and the modulating input facilitates/impedes relay of this information to motor cortex [12][13][14][15][16]. It has been observed that neural activity in the GPi changes its oscillatory patterns when a subject's cognitive state moves from being idle to planning a movement [17]. In particular, GPi activity has prominent beta band oscillations when the subject is idle, which then get suppressed when the subject plans to move. This suppression coincides with an emergence of gamma band oscillations. This would suggest, again, that one mechanism that controls relay in the motor thalamic cell is the frequency of the modulating input.
In this study, we set out to quantify when and how these thalamic cells relay driving inputs. Previous attempts to study relay neurons are made in [15,16,[18][19][20]. Specifically, in [18,19] invitro experiments are used to understand how background synaptic input modulates relay reliability of a thalamic neuron. These studies suggest that the neuron's reliability of relaying an incoming spike is governed by the background synaptic input (the modulating input) combined with intrinsic properties of the neuron. In particular [19], showed that if the variance of the background synaptic input increases, the transmission reliability goes down, and [18] showed that the feedback inhibition from the nucleus reticularis modulates the excitability of the thalamic cell membrane and hence gates transmission of spikes from the retina.
An attempt to analytically study relay neurons is made in [15], where in they studied the effects of BG inhibition on the thalamic relay reliability. They used a 2 nd order non-bursting model and phase-plane analysis to study relay neuron properties. However, they only considered a constant and a low frequency periodic modulating input. Additionally, only one deterministic periodic waveform was considered for the driving input. A follow up study with a similar objective is presented in [20], wherein the authors analyzed a relay neuron driven only by a driving input (no modulating input). Using Markov models, they studied how different distributions of driving pulse arrival times affect relay reliability. However, they did not present an explicit expression for the dependence of reliability upon input distributions and relay neuron properties.
The work presented here is different from the above computational studies in that we include classes of modulating and driving inputs in our analysis, and we employ systems theoretic tools to obtain explicit analytical bounds on reliability as a function of the neuron's electrophysiological properties (i.e., model parameters), the modulating input signal, and the driving signal parameters. Our analysis is applicable to any n th order model as long as the neuron does not spike without a pulse in the driving input and exhibits a refractory period. Consequently, our analysis is relevant for relay cells whose electrophysiological dynamics, including bursting, may be governed by several different ion channels and is more rigorous than previous works. Our lower and upper bounds contained reliability computed through simulation of both a second-and third-order model, and suggest, for example, that if the frequency of the modulating input increases and/or its DC offset decreases, then relay reliability increases.
The methods used here are generally applicable to understanding cell behavior under various conditions. In the discussion section, we show how our analysis shed new insights into motor signal processing in health and in Parkinson's disease with and without therapeutic deep brain stimulation. We also discuss how our bounds predict neural activity generated in the LGN during visual tasks with different attentional needs as well as during sleep. In particular, we show how our bounds predict the following observations in the LGN: (i) prominent c and b rhythms (25{100 Hz) in the LGN LFPs during high attentional tasks [9]; (ii) phase locking between a rhythm (12{15 Hz) in LFPs and spiking activity in the LGN in awake behaving cats [21]; (iii) h rhythms (3{5 Hz) in drowsy cats; and, (iv) even slower 0:2{0:5 Hz rhythms in sleeping cats [8].

Materials and Methods
In this section, first we describe a biophysical model of a relay neuron, and then use systems theoretical tools to compute bounds on relay reliability.

A Relay Neuron Model
A relay neuron receives two kinds of inputs: a driving input, r(t) and a modulating input u(t), and generates one output, V (t), as

Author Summary
In cellular biology, it is important to characterize the electrophysiological dynamics of a cell as a function of the cell type and its inputs. Typically, these dynamics are modeled as a set of parametric nonlinear ordinary differential equations which are not always easy to analyze. Previous studies performed phase-plane analysis and/or simulations to understand how constant inputs impact a cell's output for a given cell type. In this paper, we use systems theoretic tools to compute analytic bounds of how well a single neuron's output relays a driving input signal as a function of the neuron type, modulating input signal, and driving signal parameters. The methods used here are generally applicable to understanding cell behavior under various conditions and enables rigorous analysis of electrophysiological changes that occur in health and in disease. shown in Figure 1 B. The function of this type of neuron is to generate an output that relays the driving input at appropriate times. The modulating input does as its name implies i.e. it modulates the neuron's ability to relay the driving input [4]. This relay neuron model structure has been widely used to model thalamic relay neurons [15,16,[22][23][24][25][26].
We would like to understand exactly how the modulating input affects relay reliability of the neuron. To do so, we use a biophysical-based model to describe the electro-physiological dynamics of the relay neuron. We first begin with a second order model to highlight structure in the model dynamics, and then we present an n th order generalization. Recall that the output of the cell, dh dt~a In (1), C m ,I ion ,I ext ,V syn [ R are the membrane capacitance, ionic current, external current and synaptic reversal potential, respectively. I ion is composed of currents I T , which is a low threshold calcium ion current, and I L which is the neuron's membrane leakage current. I ext is a constant external current, and h(t) [ ½0,1, is an internal state of the system representing the probability that a calcium channel inactivation gate is open at a time t. a,g T ,g L [ R z are temperature correction factor, maximum calcium current and leakage current conductance, respectively. The details of h ? (V ), t h (V ) and m ? (V ) and numerical values used in our simulations are given in Tables 1 and 2. This is a simplified model of a thalamic neuron that is driven only by calcium ion and leak currents. We begin with this model because it is simple and still contains low threshold calcium currents which are shown to govern input selectivity of relay neurons, in a computational study [23]. This model has also been used to model neurons in the inferior olive for the purpose of studying subthreshold oscillations [27].
State space representation and n th order generalization. By defining a state vector x~½(V {V syn ),h T , an equivalent state space representation to (1) can be written as: Note that f : R 2 ?R 2 is a non linear, continuous and differentiable vector-valued function of x. In general, a state space representation takes the form _ x x~f(x,r,x), however, there is more structure in (2). From (2), one can see that f(x) is only a function of the system's internal states. The modulating input, u(t), multiplies the first component of the state x, while the driving input, r(t), is an exogenous input to the system. The 2nd order model (2) can be generalized to an n th order model to include more ion channels as well as more complicated spiking dynamics such as bursting. The n th order model is as follows:   Each g i is the conductance of the i th ion channel. V i [ R is the reversal potential of i th ion. k i [ N are such that k iz1 {k i are the number of gates in the i th ion channel and k 0~0 . Each a i is a temperature correction factor. h i ? () and t i () are functions similar to h ? () and t h ().
Inputs and outputs. For our relay reliability analysis, we assume that the two inputs belong to the following classes of signals: N Driving Input(r): This input represents the spiking activity from other neurons (e.g. cortical neurons), which the neuron must relay. Synapses of the driving input occur on proximal dendrites and are excitatory in nature. The driving input synapses are fewer in number than modulating input synapses. However, the magnitude of post synaptic potential of each driving synapse is larger compared to a modulating input synapse [4,6]. Therefore, we assume driving input belongs to the following class of functions: Here, I 0 ,t i ,t [ R z and i,n [ N. D(t) is a Dirac delta function [28]. The t i ' s are generated randomly such that t iz1~ti zT 0 zt', where T 0 [ R is a constant that represents the refractory period of driving input, and t' [ R z is exponentially distributed with probability density function: where l [ R z . The average inter-pulse interval is T~E(t iz1 {t i )~T 0 z1=l. Note that t i s are characterized completely by T and T 0 . A sample driving input is shown in Figure S1 A (supplementary material).
N Modulating Input(u): This input modulates the dynamics of the neuron and governs relay performance. Synapses of the modulating input are generally inhibitory and occur on distal dendrites. The magnitude of post synaptic potential of each synapse is smaller as compared to a driving synapse [4,6]. Therefore, this input is represented in the biophysical model (1) as a synaptic input and belongs to the following class of sinusoidal functions: Here c 1 ,c 2 ,v, and t [ R z and c 1 §c 2 . Since u(t) represents a conductance, we impose the constraint c 1 §c 2 to ensure that u(t) [ R z . Also, c 2 is appropriately small so that the modulating input does not make the relay neuron spike without a driving input pulse. This property of the modulating input will be useful when we linearize (1) for the analysis.
N We model the modulating input in a deterministic manner as it represents the ensemble sum of inhibitory post synaptic potentials (IPSPs). These IPSPs are generally small because inhibitory synapses activate the T-type Ca 2z channels allowing an influx of Ca 2z thereby reducing the magnitude of IPSPs at the soma. In relay cells, T-type Ca 2z channels have a higher density on distal dendrites [29], and this reduces the magnitude of the IPSPs even further. An ensemble effect [30] of these small IPSPs give rise to a deterministic u(t). Note that excitatory postsynaptic potentials of driving input will not get attenuated by the T-type Ca 2z channels as these channels get activated only when the cell is hyperpolarized.
N We choose the class of sinusoidal signals to shed insights into the mechanisms of oscillatory behavior or rhythms of LFPs which are often analyzed in experiments [8,9,21]. Note that LFPs arise from ensemble synaptic activity and hence may represent the modulating input. [10]. A sample modulating input is shown in Figure S1 B (supplementary material).
N Output: The output of the relay neuron is its membrane voltage V (t)~x 1 (t)zV syn .
Properties of f(x). The function f(x) is assumed to have the following 3 properties but is otherwise general: 1. Stable neuron: Consider the following undriven system: This system is the same as (4) where u(t)~c 1 and r(t)~0. Although, this system is nonlinear, we can study it via linearization about trajectories and/or an equilibrium point.
In general, a non-linear system may have multiple equilibria with different stability properties. But for our purposes, we choose f(x) such that (9) has only one globally stable equilibrium point, x x, for all pragmatic c 1 . Such a neuron is called a stable neuron [27]. This condition ensures that the neuron does not have any limit cycle, therefore, the neuron does not spike without a pulse in r(t).
This further implies that if a small periodic modulating input is applied to a stable neuron (4), u(t)~c 1 zc 2 sin(vt),r(t)~0, then after a sufficient amount of time the system's state vector will lie within a small neighbourhood of the equilibrium point. However, the state vector never reaches x x due to the time varying modulating input. The trajectory of the state in this neighbourhood can be solved using linearization methods and is periodic as we will show later. We define this periodic trajectory as the steady state orbit of a stable neuron, x o (t). See Figure 2 A.
Next, we define X ? as the collection of all points in the steady state orbit. If the initial state of the system x(t~0) 6 [ X ? then X ? is not achievable in finite time. Therefore, we relax our definition to the collection of all points inside a tube of E thickness around the steady state orbit, and define this tube as the set X r , i.e.
An illustration of equilibrium point x x, steady state orbit x o (t) and the orbit tube, X r is shown in Figure 2 A.
2. Threshold behaviour: To define threshold behaviour of a neuron, we first define a ''successful response''. A successful response at time t is a change in V such that V (t)w{50 mV & V (t{t)ƒ{50 mV Vt [ (0,L) ms. Note that both a single spike or a burst of spikes, with intra burst interval less than L ms, are counted as a single successful response under this definition. We use this definition so that we can extend our analysis to bursty neurons characterized by higher order models. Now, we state the following Lemma which defines the critical hypersurface.
Lemma 1: Given an n th order system (9), there exists a critical hypersurface of the system, x 1~S (x 2 , Á Á Á ,x n ) ƒ{50{V syn , such that x 1 (t)w{50{V syn if and only if Note that these parameters are defined by the undriven system (9). (C) Illustrates the critical hypersurface x 1~S (x 2 , Á Á Á ,x n ), a successful response trajectory, an unsuccessful response trajectory, and the refractory zone, X R for the undriven system (9). The time it takes for the solution to leave X R after generating a successful response is called the refractory period, T R . Note that refractory zone depends on I 0 and therefore T R also depends on x 1 (t 0 )wS(x 2 (t 0 ), Á Á Á ,x n (t 0 )) for some t 0 ƒt: That is, the neuron only generates a successful response if the voltage crosses the critical hypersurface (see Figure 2 C and 3).
We leave a formal proof to the reader. Essentially, by definition of S, one can show that the solution to (9) always moves away from S, unless it is on S (see Figure 3). This means that at least one of the eigenvalues of DF(x) Dx D ½S(x2,ÁÁÁ,xn),x2,ÁÁÁ,xn T has a positive real part.
This threshold property is also used in other studies [31].
which is a point on the critical hypersurface. Note that x th1 zV syn , is the traditional threshold voltage V th that people refer to for neurons [31][32][33]. In [31] it has been shown that spike threshold is influenced by ion channel activation/inactivation and synaptic conductance. In our case, the threshold x 1~S (x 2 , Á Á Á ,x n ) shows the same behavior as it is a function of the availability of activation/inactivation gates. The effect of time varying synaptic conductance is not captured by the hypersurface S. However, we used linearization methods from systems theory in section ''Response in x th Neighbourhood Under u(t)~c 1 zc 2 sin(vt)'' to include this effect. This yields a time varying threshold. Although we never explicitly deal with time varying threshold, it is implicit in our analysis. Finally, we define the threshold current, I th , such that x xz½I th ,0 T~x th . Note, by definition both I th and V th have the same units and hence can be added.
Illustrations of a successful response, unsuccessful response, the critical hypersurface S, I th , V th , x th are shown in Figure  3. Refractory period: Most neurons may generate a successful response when they are depolarized. However, they are unable to generate a successful response immediately after generating one. The duration for which they cannot generate a second successful response is called a refractory period [34]. This is because when a neuron returns back to its equilibrium point after generating a successful response, it becomes hyperpolarized, requiring extra depolarization to generate a new successful response. Additionally, due to inactivation of sodium and calcium ion gates, extra depolarization is required for the state to cross S and hence generate a successful response. This extra depolarization results in an unsuccessful response soon after a successful response.
We define the refractory zone, X R 5X c r as the region such that if x(t i ) [ X R , the neuron of type (4) (with u(t)~c 1 ) cannot generate a successful response on the arrival of a pulse in r(t) with height I 0 at time t i . Note that X c r is the complement of X r .The time spent in this zone after a successful response is the refractory period, T R . Note that, T R is not an absolute refractory period as a stronger depolarization event may result in a successful response even if x(t) [ X R .
In Figure 2 C, we illustrate X R for a second order system with f(x) given by (3) and c 1~0 :01. For this system, T R decreases with I 0 , as shown in Figure 2 E. Note that X R and X r are disjoint sets by definition.

Relay Reliability
Before we define relay reliability, we first define a relayed pulse. A relayed pulse is a successful response, V (t), that occurs within W ms after a pulse in the driving input, r(t). See Figure S2 (Supplementary Material). Let, then the empirical reliability is defined as: This definition of reliability is similar to the one defined in [15] and is not meaningful if V (t) spikes without a pulse in r(t). But since our neuron is a stable neuron, this will never happen. In the limit that we observe the neuron for an infinite amount of time, the empirical reliability converges to Let us define events SR ¼ D Successful Response to a driving pulse ð14aÞ I 0 . Additionally, note that the region shaded in the darker grey is also in the refractory zone, because if x(t i ) is in this region then AtvL ms such that We then see that Here we have used the total probability law and the definition of conditional probability [35] to go from (13) to (15). Because we cannot generate a spike in the refractory zone, X R , we get that For most neurons, the dynamics of the first component of the state, x 1 , are faster than the other states in the region (X r |X R ) c , see Figure 2 C. Therefore, when x(t) 6 [ X r |X R , it returns to X r only if it is close to X r , otherwise it returns to X R . The return process to X R is much faster as compared to the return process to X r , due to slower dynamics arising near X r . Therefore, when x(t) 6 [ X r |X R , it spends most of its time close to X r , and hence we assume that the Pr( , this assumption does not affect our results much. We will convince the reader that these assumptions are mild in the results section. Essentially, we will show that our reliability expressions under these assumptions match well to numerically computed curves for different relay neurons. Finally, since X r and X R are disjoint sets, we get: Although not explicitly in (17), relay reliability is a function of the driving input parameters, I 0 and T, the modulating input parameters, c 1 ,c 2 and v and the neuron's dynamics (i.e. model parameters) denoted by H. In the next sections, we compute closed-form approximations of lower and upper bounds of reliability as a function of I 0 ,T,c 2 ,c 1 ,v and H, by computing P response and bounds on P pulse .

Calculation of P response
To compute P response , we first find a solution for the orbit tube X r and then find a solution for the response to a driving pulse given the state starts in X r . This solution shows us when the neuron generates a successful response. We later use this information to compute P response .
The orbit tube: Response to u(t)~c 1 zc 2 sin(vt) and r(t)~0. Here, we examine the state vector response to a periodic modulating input when no driving input is applied.
The solution to (4) in the orbit tube is given by its steady state solution with r(t)~0. This steady state solution can be approximated using linearization (4) and linear time invariant (LTI) systems theory. Specifically, we linearize (4) about the nominal solution x x(t)~ x x given the nominal input u 0 (t)~c 1 . Now, if the input is perturbed such that u(t)~u 0 (t)zDu(t) and the initial condition is perturbed such that x(0)~ x xzDx(0), the state trajectory will also be perturbed to . When we substitute these values and perform a first order Taylor series expansion of (4) about the nominal solution and nominal input, we get: which can be equivalently written as: where The solution to (19) with Substituting the laplace transform of Du from (20) and defining we get: From (23), one can compute the steady state solution of (19) by taking inverse Laplace transform of (23) and taking the limit t? ?. This gives: . . .
Performance Limitations of Relay Neurons Here, %(a) denotes the angle of complex number a. Note that du(t)~c 2 sin(vt) for u(t) [ U. (19) approximates (4) in steady state when c 2 is small, which is always the case by definition of U. Also note that we will get the same steady state response even if dx o (0)=0. Using (24), we can write the steady state solution of (4) as: Now we can find the orbit tube using its definition. Figure 2 A, 4, plots the steady state orbit for a second order stable neuron with f(x) given by (3).
Response to r(t) pulses and u(t)~c 1 zc 2 in the orbit tube. We now examine the neuron's response to a driving input pulse when the solution is in the orbit tube. It is straightforward to see how a r(t) pulse affects the solution trajectory. Suppose that the state vector is at the x o (t i ) and at some time t i , when the driving signal generates a pulse, i.e., r(t i )~I 0 D(0). Then, the state vector ''jumps'' out of the orbit tube, to the point, x(t i )~x o (t i )z½I 0 ,0 T (see Figure 4). This is shown by direct integration of (4), on the time interval lim Dt?0 ½t i {Dt,t i . Now, three cases arise: 1. If I 0 wwI th , then x 1 (t i )wwS(x 2 (t i ), Á Á Á ,x n (t i )) and therefore the neuron always generates a successful response. 2. If I 0 vvI th , then x 1 (t i )vvS(x 2 (t i ), Á Á Á ,x n (t i )) and therefore the neuron never generates a successful response. 3. If I 0^Ith , then x 1 (t i )^S(x 2 (t i ), Á Á Á ,x n (t i )) or equivalently x(t i ) lies in the neighbourhood of x th . This case is biologically interesting as only for this case does the modulating input control relay reliability of the neuron. To determine whether the neuron generates a successful response or not in this case, we need to know the behaviour of the system in the neighbourhood of x th .
Response in x th neighbourhood under u(t)~c 1 z c 2 sin(vt). To approximate the response of the system in the neighbourhood of x th , we first linearize (4) about the nominal solution x S (t) (here S stands for the critical curve S) given the initial condition x S (0)~x th and nominal input u 0 (t)~c 1 . Now, if the nominal input is perturbed such that u(t)~u 0 (t)zdu(t) and the initial condition is perturbed such that x(0)~x th zdx(0), the state trajectory will also be perturbed to x(t)~x S (t)zdx(t). Note that in our case du(t)~c 2 sin(vtzvt i ) and dx(0)~dx o (t i )z ½I 0 {I th ,0 T . When we substitute these values and perform a first order Taylor series expansion of (4) about the nominal solution and nominal input, we get: In the neighbourhood of x th this system can be further approximated to: which can be equivalently written as: x and x th . When an r(t) pulse arrives, the solution jumps from x o (t i ) to x(t i )~x o (t i )z½I 0 ,0 T . Now, whether the neuron generates a successful response or not is governed by the local dynamics. Therefore, we linearize (4) about x th to analyze the behaviour of dx(t) for twt i . If a successful response is generated, At 0 wt i such that dx 1 (t 0 )w0 else if an unsuccessful response is generated At' 0 wt i such that dx 1 (t' 0 )v0.
We later show that this linear approximation does not significantly impact our expression for relay reliability of the neuron, as the numerically computed reliability fits the analytically derived curve well.
The solution to (28) is: Using the eigenvalue decomposition [37] of M, such that M~VDV {1 , V~½v 1 , Á Á Á ,v n , with each v i a right eigenvector, V {1~½ u 1 , Á Á Á ,u n T with each u i a left eigenvector and D is a diagonal matrix with eigenvalues l i at the diagonal arranged in descending order without loss of generality, we get that Note that for most stable neurons of interest, all the eigenvalues of matrix M are real. Therefore, we assume real eigenvalues for an easier read (a more messy expression can also be derived for complex eigenvalues). Recall, by the properties of f(x), that the trajectories of (4) divert away from x th , therefore l 1 must be positive. Now, if the neuron does not generate a successful response, dx 1 (t) will eventually become negative. On the other hand, if it generates a successful response, then dx 1 (t) will become positive after a sufficient amount of time (see Figure 4). The direction in which dx 1 (t) eventually moves is decided by the sign of the first component of the coefficient of e l1t . Therefore, the neuron generates a successful response if and only if Note that we substituted dx(0)~dx o (t i )z½I 0 {I th ,0 T and du(t)~c 2 sin(vtzvt i ) in (31) and integrated it to get (32). Now, we substitute dx o (t i ) from (24) into (32) and get: This equation can be written as where From (34), we see that the neuron generates a successful response if and only if Finally, we can use (36) to calculate P response , which is the fraction of the time in the orbit tube that the neuron spent in the interval in (36). This is the length of the interval divided by 2p. Therefore, Calculation of Bounds on P pulse In this section, we compute P pulse in (17) to ultimately obtain an expression for R. Since a driving pulse that arrives at time t i can only result in either a successful response or an unsuccessful response, we can equivalently write the definition of P pulse as: Pr Pr Here, we have used the law of total probability and the definition of conditional probability [35] to arrive at (38c). We know that after a successful response at t i{1 , the system state Similarly, if T us R denotes time spent in refractory zone after unsuccessful response, then we get: Now by combining (13), (38c), (39) and (40) we get: Since T us R has a complicated dependence on the input and model parameters, it is difficult to calculate Pr(t i {t i{1 §T us R ). However, it is certain that T us R ƒT R . This implies that Pr(t i {t i{1 §T R )ƒPr(t i {t i{1 §T us R ), by properties of cumulative distributive functions [35]. Therefore, we get the following bounds: Putting (41) and (42) together, we get: and P pulse ƒR : Now, we calculate Pr(t i {t i{1 §T R ). Recall that the inter pulse intervals of r(t), t i {t i{1~t 'zT 0 , here t' is generated from an exponential distribution and T 0 is the refractory period. Therefore: It can be easily shown that: T is the average inter pulse interval, E(t i {t i{1 ). Finally, by combining (43c) and (44) we get: Calculation of Bounds on R Now we compute bounds on relay reliability i.e R l ƒRƒR u . Recall that: Similarly, we can write lower bound on reliability as: Combining (47) and (48) we get: From (49) and (44), one can see that if TwwT R ,a?1, which makes R l ?R u ?P response . This result is intuitive because if pulses in r(t) occur at a slow rate, then the solution of (4) has enough time to return to the orbit tube after each pulse. Therefore, Another interesting case emerges if 0vT{T 0 vvT R {T 0 . In this case a?0 and R l ?0,R u ? P response 1zP response . This case has two interesting extremes: 1. T 0 ?0, making TvvT R , 2. T 0 ?T R , and both T{T 0 and T R {T 0 approach 0. In case 1, an average a T R =T? ? number of pulses occur in the T R time interval after a successful response. All of these pulses generate unsuccessful responses because the system state is inside X R during this interval. Therefore, for each successful response, we get ? unsuccessful responses making R?0~R l . However, in the second case, exactly one pulse occurs during the T R period after a successful response.
Therefore, for every successful response we get at least 1 unsuccessful response. Now, if P response~1 , we get exactly one unsuccessful response for each successful response making R? P response 1zP response~1 2 .

Results
In this section we verify our reliability bounds by simulating a second and third order model for a thalamic relay neuron.

Dependence On Model Parameters
The dependence of reliability on the cell's input parameters is explicit in our bounds. However, dependence of reliability on the model parameters is captured implicitly by the gain G, I th and T R . The refractory period, T R , is well studied in literature and depends on inactivation gate time constants [38]. Therefore, in this section we discuss how the gain G(v) and I th depends on the properties of a relay neuron membrane dynamics.
In Figure 7 A, we plot I th vs conductances g T and g L . We see that I th first decreases with increasing g L and then increases forming a parabola. Furthermore, with increasing g T , I th decreases. In Figure 7 B, we plot the dependence of the gain G(v) on g T and g L . G(v) is essentially a low pass filter whose amplitude decreases as frequency increases. Consequently, reliability increases with frequency (see (49)). From the Figure, we can see that the gain, G(v), in the high frequency range (v=2pw100Hz) increases with g L and decreases with g T . For lower frequencies, vv100Hz, G(v) has a complex dependence on g T & g L . This is an important result as we can increase/ decrease reliability of the relay neurons by increasing/decreasing T-type Ca 2z or leak channel conductances which can be further used to treat diseases such as Parkinson's disease (see discussion).

A 3 rd Order Model of a Thalamic Neuron
In this section, we will apply (49) to a third order model of a thalamic relay neuron. In this case, the parametrs I th ,G(v),T R in the equation are computed from the third order model.
We chose the 3rd order thalamic model used in [15,16,22], which is a simplification of model used in [39,40]. This model exhibits bursting activity in the hyperpolarized state and non bursty firing in the depolarized state. The two responses of the model for an oscillating modulating input and a Poisson driving input (inter-pulse interval is given by (7)) are shown in Figure 8 A and 9 A. The equations and parameters of the model are the same as those used in [15,22]: In the (50), I K~gL (0:75(1{h)) 4 (V {V L ) are the leak current, sodium and potassium current, respectively. I T~gT p ? 2 (V )r(V {V T ) and I ext are the low threshold potassium current and external current respectively. a 1~1 ,a 2~2 :5 are the temperature correction factors. All the parameters used are given in Table 3.
A thalamic neuron generates a single spike when depolarized in the relay mode [15,41]. However, it generates a burst of spikes when it receives a depolarizing input when it is in a hyperpolarized state [42]. We used I ext~{ 0:56, to model the hyperpolarized or bursty state. Whereas, I ext~0 models a single spike state of thalamic neuron. We can rewrite the (50) in the form of (4) by defining the state vector X~½V {V syn ,h,r T with: f (X )~: In Figure 8 A, we plot the time profile of the voltage for a bursty neuron along with a zoomed in view of the burst in Figure 8 B. Figure 8 C plots our reliability bounds (49) along with empirical reliability computed numerically through simulation of the 3rd order model. We see that our bounds predict reliability well even for a bursty neuron. Note that we consider a burst response to a pulse as a successful response.
In Figure 9 A, we plot the time profile of voltage for a non bursty neuron along with a zoomed in view of a successful spike in Figure 9 B. Figure 9 C plots our reliability bounds (49) along with empirical reliability computed numerically through simulation of the 3rd order model. Note that here T R vT 0 therefore R l ?R u . We see that our bounds predict reliability well in this case also.
In general, our analytical bounds are applicable as long as the model 1. does not generate a spike if there is no pulse in r(t), and 2. has a threshold behaviour as defined in Materials and Methods section, and 3. shows a refractory period. The second condition is true for most neurons that satisfy the first condition. Our analysis may also be extended to include neurons that spike without any driving input (see Discussion), but in this manuscript we neglect such dynamics.

Discussion
In this manuscript, we studied the reliability of a relay neuron. A relay neuron receives two inputs: a driving input, r(t), and a modulating input, u(t). The neuron generates one output, V (t), which relays r(t) conditioned on u(t). Our goal was to precisely determine how the modulating input impacts relay reliability. To calculate relay reliability, we used systems theoretic tools to derive the analytical bounds (49) on relay reliability as a function of different input and model parameters. Specifically, (49) implies that if the modulating input is of the form u(t)~c 1 zc 2 sin(vt), then increasing c 1 or c 2 decreases reliability. However, increasing v increases reliability. In addition, our reliability curve (see Figure 5) suggests that reliability first increases slowly with v and then increases rapidly and plateaus. (49) is powerful as it characterizes the multiple dependencies of reliability on u(t),r(t) and relay neuron model parameters. Furthermore, analytic bounds from (49) contain results obtained through simulation of the 2 nd and 3 rd order models of a relay neuron. Our bounds captured reliability under both the depolarized and hyperpolarized states of the 3rd order neuron and shows the generality of our analysis.

Comment on Spontaneous Firing in Relay Neurons
Our reliability bounds were calculated assuming that the relay neuron does not fire spontaneously. However, many relay neurons show spontaneous firing in the absence of any input. This spontaneous firing is usually periodic (period T n ) because it arises from the emergence of a limit cycle [43] and can be thought of as responses to background noise. Our analysis can therefore be extended to capture this by adding a periodic noise pulse train n(t) ¼ D P i I n D(t{iT n ) in the reference input r(t), therefore the new reference input becomes: Since a successful response to a pulse in n(t) is undesirable, we must modify our definition of reliability. To do this, we assume that the arrival of a pulse in n(t) cannot coincide with an arrival of a pulse in r(t) and thus successful responses to pulses in each signal are disjoint events. This leads us to define reliability as R~Pr(Successful Response to a r(t) pulse) |Pr(r(t) pulse) {Pr(Successful Response to a n(t) pulse) With this approach, our analysis can be extended to spontaneously firing neurons. We believe that the reliability will approximately be bounded as: The above expression is reduced to (49) in the case T n wwT i.e the noise period is much larger than the period of the driving input. In the case when T n vvT, the reliability becomes negative because noise pulses occur very frequently as compared to desirable driving input pulses. This generates undesirable successful responses making reliability negative. Note that (54) is only an approximate solution for the reliability of spontaneously firing relay neurons and we leave the exact solution to this problem for the future work. Plots I th as a function of g L ,g T B. G(v) (see (35) versus g L and g T . Note that G(v) depends largely upon g L , whereas its dependence upon g T is minimal. g T changes the maximum value of G(v) but does not effect it much in the high frequency range. doi:10.1371/journal.pcbi.1002626.g007  (48) and (47), respectively. The solid line is plots R emp calculated by running simulations of (4), and the error bars indicate +std. We estimated I th~8 :7126 as the minimum height of a r(t) pulse that makes the neuron generate a successful response. doi:10.1371/journal.pcbi.1002626.g008 Table 3. Parameters and functions for (50).

Function
Value

Motor Signal Processing
In the motor circuit, thalamocortical neurons receive a driving input from the motor cortex and a modulating input from the GPi segment in the basal ganglia (BG). See Figure 10 A. The function of the GPi input is hypothesized to enable/disable thalamic cells to relay cortical stimuli related to movement when movement is intended/not intended [14]. This is consistent with evidence that the BG both inhibits unwanted movements and enables intended movements in a timely manner [12,13]. This GPi modulated thalamic relay ultimately enables reliable transfer of information from higher cortical layers to lower layers which then command the musculoskeletal system to generate planned movements [44]. The thalamic relay hypothesis is supported by previous studies [4,16,22]. In [16,22], it is shown that relay reliability computed from a data-driven computational model of a thalamic neuron is low in Parkinson's disease (PD), and high in both healthy and when therapeutic DBS is applied to the BG in PD.
Previous works emphasize the inhibitory projections from GPi to motor thalamus [45][46][47][48]. They argue that when movements are intended/not intended, appropriate task-related GPi neurons decrease/increase their firing rates. This in turn disinhibits/ inhibits thalamus and consequently enables/disables thalamic relay, respectively. Our analysis as well as recent experimental observations show that the story is a bit more complicated. GPi firing rates alone may not be the mechanism for thalamic relay, rather, the dynamics of the GPi activity control thalamic relay. In particular, it appears that the oscillatory dynamics of GPi activity control relay. Our relay bounds predict that if one intends to move, then the GPi neurons that project to motor thalamus should initially generate LFP activity that has prominent low frequency oscillations which allows the subject to remain idle, and then generate activity that has prominent high frequency oscillations which allows the subject to plan an intended movement and then move.
We first discuss how our analysis concurs with observations obtained from a computational model of the motor circuit that characterizes neural activity dynamics in the BG and motor thalamus in health and in PD with and without therapeutic DBS. The computational model simulates neural activity when movements are planned and hence when motor thalamus should relay information from the cortex at all simulated times. We then discuss how our relay bounds accurately predict how GPi activity recorded from two healthy primates modulates during a structured  (48) and (47), respectively. Note that here T R vT 0 , therefore R l ?R u . The solid line plots R emp calculated by running simulations of (4), and the error bars indicate +std. We estimated I th~7 :0155 as the minimum height of a r(t) pulse that makes the neuron spike. doi:10.1371/journal.pcbi.1002626.g009 behavioral task that forces an idle phase, and a planning phase during each task trial.
Predicting data from a computational model. In PD, the GPi input to thalamus becomes pathological and prevents the thalamus from properly relaying information back to the cortex. In particular, people have observed pathological 10-30 Hz beta rhythms and synchronization emerging throughout the BG in PD [49][50][51][52]. High frequency DBS (HFDBS) modulates activity in the BG structures, including GPi, and may restore thalamic reliability leading to clinically observed reversal of symptoms in PD [51,53].
To better understand how HFDBS may restore relay reliability, we first consider a computational study [15] of basal-gangliathalamic neural signal processing. In [15], a biophysical-based model of multiple BG structures and motor thalamus is constructed and parameters are tuned to generate 3 states: healthy, PD and PD with HFDBS applied to the subthalamic nucleus (STN) in the BG. In Figure 10 C, we reproduce plots from this study that illustrate the simulated GPi modulating input to thalamus in the 3 states. We then discuss how our reliability bounds predict what is observed in these simulations.
N According to the computational model in [15], in the healthy case, relay reliability is high. When we look at the simulated Gpi activity, the amplitude of the GPi modulating input, c 2 , is small enough to generate reliable thalamic relay of cortical inputs in accordance to our bounds (49). As Figure 10 D shows, the orbit tube is small for such a c 2 , which results in less time spent in the unsuccessful response region, X us . Physiologically, a small c 2 may be due GPi neurons being uncorrelated so that when they add they do not produce large LFP amplitudes. Gpi neurons have been observed to be uncorrelated in healthy primates [54,55].
N According to the computational model in [15], in PD, reliability is low. When we look at the simulated GPi activity, the amplitude of the GPi modulating input is larger than in the healthy case, which leads to a lower relay reliability at the thalamus according to our bounds (49). As Figure 10 D shows, the orbit tube is large for a larger c 2 and results in more time spent in the unsuccessful response region, X us . In Figure S3 (Supplementary Material), we have plotted R versus c 2 and v. From the Figure S3, it is clear that when c 2 increases, reliability decreases, whereas when v increases reliability increases. Physiologically, a large c 2 may be due to GPi neurons being synchronized so that when they add their peaks sum producing large LFP amplitudes. Synchronization of neurons in the BG has been observed in PD patients [56] and parkinsonian primates [49,55].  [15] for the Healthy, PD and PD with high frequency deep brain stimulation (HFDBS) cases. As we can see in the healthy case, the amplitude of the BG output, c 2 , is smaller compared to the PD BG output, resulting in a higher relay reliability. HFDBS increases the frequency, v, of the BG output, resulting in a higher relay reliability. (D) Intuition of how reliability changes in the three cases. In PD, c 2 is larger, therefore, the diameter of the orbit tube is larger compared to the orbit tube for healthy. This results in more time spent in the unsuccessful response region X u s, which leads to poor reliability. In contrast, in PD case with HFDBS applied, v is larger and the gains g i (jv) decrease, which generates a smaller orbit tube. In this case, the state spends more time in the successful response region X s of the orbit tube, resulting in high reliability. doi:10.1371/journal.pcbi.1002626.g010 N According to the computational model in [15], HFDBS applied to the PD model restores thalamic relay. When we look at the simulated GPi activity under HFDBS applied to the STN, the frequency of the GPi modulating input increases and the amplitude, c 2 decreases and is more comparable to that in the healthy state. This combination (increased v and decreased c 2 ) restores relay reliability according to our bounds (49). As Figure 10 D shows, the orbit tube is small for a larger v because a larger v results in a smaller g i (v)s and hence a small G(v). Recall that G is proportional to the diameter of the orbit tube (see Figure 10 D).This results in less time spent in the unsuccessful response region, X us . Note that the frequency of DBS is not directly related to frequency of modulating input. One can see from Figure 10 C, that modulating input frequency is only 80 Hz while HFDBS frequency is 140 Hz.
The working mechanisms of HFDBS in PD are thought to be (i) suppression of pathological beta oscillations in the BG [49,51], and (ii) desynchronization of BG neurons [22,57]. Our analysis accurately predicts this, but also highlights new possible therapies. For example, as discussed in section Dependence On Model Parameters, the conductance of leak channels is critical for a relay neuron because it decides the size of the orbit tube for a given v. In particular, for smaller v, the gain G(v) decreases as we decrease g T , which results in increased reliability. This suggests that if we could pharmacologically decrease g T , a lower frequency (hence lower power) DBS signal may be therapeutic. A low power DBS option would save battery power as well as minimize side effects associated with high power stimulation [58][59][60][61]. There are many ways to regulate the conductance of T-type calcium channels, reviewed in [62]. These methods include (1) hormonal regulation by dopamine, serotonin, somatostatin, opioids, ANP, and ANG II. (2) Guanine nucleotides (3) Protein kinases (4) voltage. To be target specific, these methods may require injecting the chemical directly into the thalamus.
Predicting data from experiments. The computational model in [15] does not capture the subject's intent. It is assumed in [15] that the subject is moving and that the sensorimotor cortex sends a driving input to a thalamic cell accordingly. In reality, a subject's motor program is coordinated in time. When a subject is idle, then the activity of GPi neurons (modulating input) should have slower oscillatory patterns according to our analysis so that the thalamus does not relay information. Furthermore, when the subject plans to move, the task-related GPi neurons should then generate more high frequency oscillations to enable relay of this movement via the thalamus. This can be understood more clearly by looking at the Figure 10 A. When the subject is idle, u(t) should be a low frequency signal and when the subject plans to move, u(t) should change to a high frequency signal.
This has been observed recently, when we showed that taskrelated GPi neurons indeed exhibit a ''crossover effect'' during movement planning in two healthy macaque monkeys executing a directed hand movement task [63] (see Figure S4, Supplementary Material). Initially, when the monkey is idle, there are prominent 10-30 Hz beta oscillations in the neuronal spiking activity. Then, when a final cue is given to indicate what movement should be executed, gamma band oscillations  Hz) emerge in the spike trains of GPi, displaying the ''crossover'' (beta gets suppressed while gamma emerges) [63].
Futhermore, if GPi's mechanism in motor control is to modulate its oscillatory rhythms in a timely fashion as our relay analysis predicts, then the prominent beta oscillations observed in PD [49,52,55,56] may partially block this mechanism. That is, in PD it may be more difficult to suppress beta during movement planning as it is so prominent, leading to poor thalamic relay and poorly generated movements.
Finally, we highlight another recent study of ours that showed that when therapeutic DBS (w100 Hz) is applied to the STN of a parkinsonian and healthy primate, then the propensity of GPi neurons to spike in the gamma band increases [17]. This finding, along with the above observations, indicate that perhaps the mechanism of HFDBS is to re-enable the crossover effect in GPi (i.e. increase gamma oscillations to overcome the prominent beta oscillations) that controls thalamic relay and movements in PD.

Visual Signal Processing
As mentioned in the Introduction, neurons in the LGN receive driving input synapses from the retina and modulating input synapses from layer 6 of the visual cortex and the brain stem. The LGN then relays the driving input to visual cortex for perception. The LGN functions as a ''gatekeeper'' and allows only the relevant information to go through depending on attentional demands [7,64]. In the LGN, the spatial map of the visual field is conserved [64,65].
Here, we hypothesize that the LGN functions as a filter of the spatial map which shows a high relay reliability in spatial areas requiring high attention and lower reliability otherwise. Our analysis suggests then that LGN neurons relaying attended areas of the visual field receive higher frequency modulating inputs as compared to LGN neurons relaying areas which are ignored. Note that the modulating input represents the synaptic background activity, which is a major contributor to LFPs and EEG recordings [10]. Therefore, the frequency content of LFPs and EEG reflect the frequency of the modulating input.
This hypothesis is supported by [8], where it was shown that the frequency of the LFPs in LGN depends on the arousal state of the cat. Specifically, they showed a prominent a rhythm (12{15 Hz) in awake and naturally behaving cats, a h rhythm (3{5 Hz) in drowsy cats and a slow rhythm (0:2{0:5 Hz) during sleep. Additionally [21], showed that, in wakeful naturally behaving cats, the spiking activity of relay-mode (non-bursty) neurons in the LGN is correlated with the phase of the alpha rhythm of the LFPs. Specifically, some neurons spike more at the peaks of the alpha wave while other neurons spike more at the valleys of the alpha rhythm. Using (36), we may be able explain why such phase locking occurs. In words, this equation says that relay neuron reliably relays the driving input only during a fixed phase interval of modulating input, and this phase interval depends on neuron membrane properties [21].
Finally, during deep sleep slow delta rhythms are observed in the EEG which are believed to be of thalamic origin [66]. This may cause even lower reliability in LGN and filter out all the visual information, resulting in deep sleep.
On the other hand, high frequency b (16{24 Hz) & c (30{45 Hz) rhythms are observed during visual attentional tasks in the LFPs of cat LGN [9]. Our analysis shows that reliability increases with modulating input frequency, therefore we propose that the reliability during these tasks is greater than during natural wakeful behaviour for most LGN neurons. This results in larger relay of information which increases general productivity.
In addition to the observed relationship between the LGN LFP oscillations and attention, it has been observed that during sleep, LGN neurons become hyperpolarized [42,67]. In our model, this means that the DC offset of the modulating input, c 1 , is large which decreases reliability according to our analysis. The LGN neurons relay poorly and also exhibit a bursty behaviour (see Figure 6 A and 8). The lower reliability may result in less information relay from the LGN to the visual cortex, inducing sleep whereas the bursty behaviour may only be a by-product of hyperpolarization and may have nothing to do with information suppression. This agrees with [68] where in it is shown experimentally that although all bursts combine carry lesser information than all single spikes, individual burst is more informative than a single spike in the LGN output. The information carried in the bursty mode may be critical for waking up [42]. Figure S1 Driving and modulating inputs. (A) Driving Input. Each arrow denotes a delta pulse with height I 0 . The interpulse interval is t i , which is an exponential random variable. B Modulating Input. A sinusoidal wave with DC value c 1 , amplitude c 2 and frequency v.

Supporting Information
(TIFF) Figure S2 The relayed pulse and the non-relayed pulse. A pulse is called a relayed pulse if it generates a successful response in V (t) within a W ms window, otherwise it is called an ineffective pulse. (TIFF) Figure S3 Reliability vs c 2 and v. Note that increasing c 2 decreases reliability whereas increasing v increases reliability. (TIFF) Figure S4 Neurons in GPi show a crossover effect during the planning phase. In this experiment, two primates executed a directed hand-movement task. From the figure, we can see that the percentage of neurons displaying more power in gamma band compared to beta band increases just after the final command is given. This may be because the GPi output is the modulating input to the relay neurons in motor thalamus, and a increase in the frequency of the modulating input may allow a certain motor plan to be relayed back to cortex and downstream to brain stem to ultimately get executed. This figure has been taken from [63]. The thin solid and dotted lines are the 5% and 95% confidence bounds obtained by randomization of the spike trains. (TIFF)