Representation of Dynamical Stimuli in Populations of Threshold Neurons

Many sensory or cognitive events are associated with dynamic current modulations in cortical neurons. This raises an urgent demand for tractable model approaches addressing the merits and limits of potential encoding strategies. Yet, current theoretical approaches addressing the response to mean- and variance-encoded stimuli rarely provide complete response functions for both modes of encoding in the presence of correlated noise. Here, we investigate the neuronal population response to dynamical modifications of the mean or variance of the synaptic bombardment using an alternative threshold model framework. In the variance and mean channel, we provide explicit expressions for the linear and non-linear frequency response functions in the presence of correlated noise and use them to derive population rate response to step-like stimuli. For mean-encoded signals, we find that the complete response function depends only on the temporal width of the input correlation function, but not on other functional specifics. Furthermore, we show that both mean- and variance-encoded signals can relay high-frequency inputs, and in both schemes step-like changes can be detected instantaneously. Finally, we obtain the pairwise spike correlation function and the spike triggered average from the linear mean-evoked response function. These results provide a maximally tractable limiting case that complements and extends previous results obtained in the integrate and fire framework.


Introduction
Intracellular recordings of multiple neurons have shown that dynamical sensory stimuli can modulate input currents to cortical neurons. For example, visual stimulation with moving gratings can lead to oscillatory modulation of the membrane potential in visual cortical neurons [1]. While it is important to identify the role such dynamic modulations play in neural coding [2][3][4], it is also important to understand how their encoding depends on physiological parameters such as the background noise or the firing rate. The incoming external signals can be encoded in the mean or variance of the synaptic current to each neuron in a cortical network. In a single neuron the maximal firing rate limits the highest faithfully encoded frequencies. Yet the spike rate is remarkably low, often below 1Hz in cortical neurons [5,6]. Therefore, the representation of perceptually important fastvarying stimuli has to emerge at the population level. In neuronal populations the frequency response function quantifies the fidelity of signal representation [7][8][9][10]. Intuitively, the frequency response function measures how well the population firing rate is modulated by the incoming signal of a specific frequency. If the amplitude of the rate modulation is zero then that frequency cannot be encoded in the population rate.
An early study by Knight showed that a population of independent (perfect) integrate and fire neurons can faithfully encode any input frequency, but if finite memory is introduced to the single neuron dynamics, the population rate is no longer a perfect copy of the stimulus [7]. Subsequent studies established that two factors play a particularly important role for the frequency response: the noise statistics and the spike generation mechanism [9][10][11][12][13]. Brunel and colleagues have shown that in the leaky integrate and fire model high-frequency mean-modulating signals are represented faithfully in the population rate only on the background of colored noise [10]. Substituting colored for white noise background, on the other hand, leads to 1= ffiffi ffi f p decay of the frequency function for input frequencies f much larger than the stationary firing rate [9,14]. The analytical complexity of the colored noise results, however, necessitates a largely numerical treatment and generally allows for explicit expressions in specific limits only, e.g. the linear regime of weak amplitudes [10,15]. One notable exception are the recent results obtained by Ostojic, Richardson and colleagues for the non-linear response functions of the exponential integrate and fire model [16,17]. Notably for variance modulations, the leaky integrate and fire model neurons can faithfully encode any high-frequency input even for white noise [9,14]. So far, only specific limits of the frequency response function could be calculated in the integrate and fire framework for a limited set of temporal correlation functions. The linear response for mean modulating signals has been obtained analytically in the limit of white or almost white Ornstein Uhlenbeck currents by Brunel and colleagues in the leaky integrate and fire model [10]. Also Brunel and Latham obtained a linear response for mean modulating signals in the limit where the correlation time is much larger or much smaller than the membrane time constant in the quadratic integrate and fire model [18]. The high-frequency limit of the response function for mean modulations has been studied in various integrate and fire type models. Brunel and colleagues obtained in the exponential and quadratic integrate and fire model the response amplitude in the high-frequency limit and showed that it decay !1=f for frequencies f beyond the cut-off determined by the inverse spike initiation time [11,15,16,18,19]. More generally, integrate and fire type models with a variable spike onset initiation time have shown that for white as well as colored noise the representation of high-frequencies in the mean channel is successively enhanced if the spike onset time is reduced [11,20]. For the variance modulation, however, the analytical results appear much more sparse in the literature and are available so far only in the white noise limit and in the linear regime of the leaky integrate and fire model [9,21] or the perfect integrate and fire model [22]. In summary, these previous model studies have shown that the structure of the noise background can fundamentally change the response properties-particularly the difference between a perfect white noise and colored noise can be profound. The sharpness of the spike onset has little effect on the low and intermediate frequencies but strongly determines the high-frequency cut-off above which the decay of the frequency response function sets in [11,12,20]. Notably, recent experimental evidence indicates that cortical neurons can indeed encode input frequencies that are tens of times faster than the firing rate of individual neurons, in both mean-and variance-encoding schemes in the presence of in vivo-like correlated background noise [23][24][25]. This suggests that a threshold-based model that is driven by different types of colored noise can be a promising starting point to understand the fundamental determinants of the frequency response function for the physiologically important intermediate frequency range up to a few hundred hertz [23][24][25].
Here, we show that an alternative threshold-based framework can be used to obtain explicit and tractable results for the linear as well as non-linear response to mean-and variance-encoded stimuli. The explicit results derived here for the frequency response function, pairwise spike correlations and spike triggered averages constitute a maximally tractable limiting case that complements and extends the results obtained in the integrate and fire framework. Importantly, this framework does not limit the accessible current correlation functions to white noise or Ornstein Uhlenbeck process, thereby allowing us to explore a wide variety of shapes and time constants that can occur in vivo [26][27][28].
The manuscript is organized as follows: we start with the introduction of the mean and variance signaling in cortical networks in Section ''Signal representation in cortical networks''. We then introduce the population firing rate response dynamics in Section ''Key definitions of dynamical population response''. Subsequently, we define the model setting and compute basic quantities such as the firing rate of individual neurons in Section ''Dynamical response in threshold model neurons''. In Section ''Population response to mean and variance oscillations in the threshold model'' we obtain the population rate dynamics in response to oscillatory changes of mean and variance. In Section ''Response to step-like input current changes'' we address the population response to step-like input changes. In Section ''Weak pairwise spike correlations'' we quantify the spike correlations in two neurons that are subject to a common fluctuating mean signal. Finally, we focus on the statistics of spike triggering events in Section ''Spike-triggering events''. In discussion section we present a discussion of our results and their relation to previous theoretical and experimental findings. A nomenclature overview can be found in Table 1.

Signal representation in cortical networks
Neurons in the mammalian cortex form an interconnected network, where each neuron receives inputs from many thousands

Author Summary
Sensory stimuli in our environment are represented in the brain as input current changes to neurons. For example, a periodic bar pattern in the visual field leads to periodic current modulations in the visual cortex. Therefore, models describing the ability of neurons to represent incoming stimuli can offer important clues about how sensory stimuli are processed by the brain. As anyone who has used an old-fashioned radio can attest, there is not just one but multiple ways to encode a signal, e.g. the familiar AM and FM channels. But what are the potential encoding channels in the cortex? A signal could modify the neuronal input current in two distinct ways: it could act either on the mean or the variance of the current. Using a minimal model framework, which can reproduce many features of neuronal activity, we find that both encoding schemes could be equally potent in transmitting slow and fast signals. This allows us to describe how input signals of any functional form give rise to collective firing rate changes in populations of neurons.  of presynaptic neurons. The excitatory and inhibitory inputs at each neuron counteract each other [29,30]. What results is an excitation-inhibition balance which is schematically illustrated in Fig. 1 A(left). In Fig. 1 A (left) the mean excitatory (grey) and mean inhibitory (black) currents are counteracting each other and result in a zero-mean net current at the soma of a neuron. Yet, the subtraction of excitatory input by inhibition is not perfect and the remaining net current has a sizable variance. What could be the benefit of operating in this way? Theoretically, it is understood that a neuronal population in such a state could encode and relay incoming signals via two channels (1) modifications of the mean synaptic bombardment and (2) modification of the synaptic fluctuation variance. These two encoding channels are schematically illustrated in Fig. 1 (A). In case (1) the signal is added to the input of neurons, in case (2) the signal modulates the variance of the background fluctuations in neurons, similar to the amplitude modulation strategy which is widely used in radio communication.
To employ strategy (1) a sensory stimulus could alter the mean current by adding an external signal to the network generated background fluctuations. On the other hand, in a cortical population where the fluctuations in the activity of excitatory and inhibitory populations accurately track each other [30,31] the effect of excitation and inhibition can be precisely balanced. In this case, any perturbation would result only in a change of input variance to each neuron and a balance of inhibition and excitation would compensate any mean current changes. To rapidly encode modulations of the input current in the population firing rate, the strategy employed by neurons needs to be susceptible to subtle changes of either mean or variance. Now let us formalize the definitions and define the mean and variance of net input currents in a cortical network in the absence of external signals. According to the calculations outlined in Methods Section ''Current mean and variance in a cortical network'' the input current mean and variance for all neurons in a cortical network can be described by: In the mean channel, the activity of excitatory neurons needs to be increased by A while the activity of inhibitory neurons is decreased by the same amount A. Fig. 1 (right, green) illustrates that this procedure modifies the mean but leaves the variance of the net current unaffected.
In a second independent encoding scheme, a signal is encoded in the variance of the net current. As demonstrated in Fig. 1 (right, red) this could be achieved by increasing the excitatory and inhibitory firing rates simultaneously by an amount A: dI mean~0 : ð8Þ It is plausible to assume that in live cortical networks, both mean and variance signaling act simultaneously. For the physiologically plausible regime of small but simultaneous mean and variance changes, the resulting population response is simply a superposition of the mean and variance responses. For large simultaneous modulations of mean and variance, the resulting population response needs to be computed on a case-by-case basis, because it can potentially depend on the stimulus form, and the amplitude ratio of mean and variance signals.
Effective independence assumption for cortical populations. To understand how interconnected cortical populations respond to mean or variance signals, two ingredients are crucial: 1) response of independent neurons and 2) the influence of topology. Before 2) can be addressed, the response of independent cortical neurons in 1) needs to be clarified. Notably, several recent studies have shown that the net input current and spike correlations observed in vivo cortical are predominantly weak, with spike correlation coefficients v0:1 in awake animals [5], see Fig. 1.
Additionally, recent studies have proposed active decorrelating mechanisms based on a precise temporal correspondence between excitation and inhibition in each neuron in a network, which result in pairwise decorrelated spiking activity [31,32]. These exceptionally weak correlations suggest that the assumption of independent neurons can be a promising starting point for the investigation of response dynamics of cortical populations. Therefore, we focus in this study on the population firing rate dynamics of independent neurons subject to mean or variance modulating signals. The fluctuating background currents, which in a live network results from a sum of inhibitory and excitatory currents, will be synthesized independently for each neuron as a random realization of a random Gaussian process and modulated either in their mean or variance. In the following Sections we introduce the threshold-based single neuron model that we use to characterize the dynamical response in such a population of independently encoding neurons.

Key definitions of dynamical population response
Before we start with the specific representation of dynamical signals in the threshold model, let us first specify modelindependent definitions for the dynamical population response, such as linear response function or spike triggering events. A nomenclature summary can be found in Table 1.
In the previous Section we established that signals in a cortical network can be encoded in the mean or in the variance of input currents. If this signal constitutes only a small perturbation of the system, as can be expected for example from thalamo-cortical projections [33], then the population firing rate n(t) is directly proportional to the signal A : exp(ivt): The proportionality factors n m,1 (v) or n v,1 (v) are referred to as the linear response functions. Each of these functions essentially describes how each signal frequency is affecting the firing rate. Note, that Eqs. 11 and 12 are model-independent, and that linear response dynamics for weak enough signals can be derived for any nonlinear system. Also note that the linear response derives its name from the fact that the output rate n(t) is linearly related to the input in Eqs. 11 and 12. Fig. 2 (top) schematically demonstrates the encoding of periodic mean signals in the presence of background noise. The resulting spikes collectively lead to a periodic modulation of the firing rate described by Eq. 11. In addition to the oscillatory changes in mean and variance another group of changes bears particular physiological significance. These are the step-like changes in mean and variance. To compute the population rate response to a step-like signal, we first formally describe that by a Heaviside h-function: . This model independent relation describes the contribution of each frequency to the steplike signal. To compute the neuronal response to the mean-encoded weak step-like stimulus, it is sufficient to consider each frequency separately, compute the respective output and sum up all contributions. The response to a step signal of size A can then be formally written as a convolution: In the case of mean channel and weak signals, the linear response function and its inverse Fourier transformn n m, are the only functions [34] we need to predict the population response to any weak dynamical current stimulus I(t): For signals which constitute a larger dynamical perturbation, successively higher order Wiener kernels [34] may be necessary. Furthermore, the linear response function determines not only the response to dynamical signals but it also shapes the weak spike correlations in two statistically identical neurons driven by a fluctuating Gaussian current [34]: n cond (t)~Ss 1 (t)s 2 (tzt)T=n, Here, C I (v) is the Fourier transform of the current correlation. The input current correlation strength is 0vr%1. Eqs. [11][12][13][14][15][16] show that the response to variance and mean modulating signals are critical to a number of phenomena, from the processing of periodic stimuli to inter-neuron synchronization.

Dynamical response in threshold model neurons
Here, we use a previously introduced threshold-based model [35][36][37][38] that identifies the spike times t j of a neuron with positive threshold crossings of a correlated, stationary Gaussian voltage V (t) as it is illustrated in Fig. 3 (left): S : T denotes the ensemble average and y 0 the spike threshold. A nomenclature overview can be found in Table 1. The voltage correlation function C(t) is characterized by its peak value C(0)~s 2 V and the finite correlation time constant t s~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The correlation time t s , which is related to the autocorrelation time [37], describes the width of the correlation function in the vicinity of zero and is proportional to the zero crossings of a parabolic fit to the correlation function in the vicinity of zero. The second derivative of the correlation function C''(t) has a variance s 2 For simplicity we will assume a linear Resistor Capacitor (RC) membrane filter such that the Fourier transforms of the voltage and current correlation functions are related through Our results, however, are not specific to this current voltage transformation, but can accommodate other filters that give rise to smooth voltage correlation functions [37,38]. Note that this model does not have fixed-reset condition, but a silence period after each spike emerges from the regularity of the voltage trajectory [37]. The firing rate of a single threshold neuron is: where y 0 is the threshold voltage, d( : ) and h( : ) are the Dirac delta and Heaviside theta functions, respectively. In this approach the firing rate n is a particularly tractable expression which depends only on two parameters: the correlation time t s and the threshold-tovariance ratio y 2 0 =s 2 V , but not on the specific functional choice of the correlation function. Unless stated otherwise, we use the correlation function C(t)~s 2 V =cosh(t=t s ) which is compatible with the power spectra of cortical neurons [26] for simulations using digitally synthesized Gaussian processes [39] or numerical integration of Gaussian integrals. Numerical simulations were implemented in Matlab 2010a (The MathWorks, USA) and analytical calculations were partially implemented in MATHEMATICA 5.2 and 8.0 (Wolfram Research Inc, USA). Let us note that the dependences on the threshold or current variance in Eq. 18 are consistent with the predictions in the leaky integrate and fire model, as we have shown in [37]. Note, that this threshold model operates only in the fluctuation driven, low firing rate nv1=(2pt s ) regime that is particularly important for visual cortical neurons [5,6,29,30]. The mean-driven regime escapes the validity regime of this model.
Population response to mean and variance oscillations in the threshold model. We compute the population rate dynamics in response to changes of the mean and variance for the full range of input frequencies. First, we start with the mean modulations. As shown schematically in Fig. 3 (B) this paradigm subjects each neuron to a current which consists of a periodic signal and a fluctuating background background noise which is unique to each neuron. To calculate the response evoked by a sinusoidal modulation A : exp(ivt) of the mean in the threshold model framework, the signal is low-pass filtered in the cell membrane: Following the steps outlined in the Methods Section ''Mean modulation'', Eq. 17 can be modified to accommodate the oscillating voltage-to-threshold distance. The complete non-linear population firing rate response reads: Note that the firing rate modulation in response to periodic mean current is independent of a particular functional form of the current correlation function, because C(t) does not enter Eq. 20. We also find that Eq. 20 can be generalized to describe the population response to any signal k(t). In this case the membrane filtered signal f (t) in Eq. 20 only needs to be replaced with the corresponding solution of the inhomogeneous differential equation in Eq. 19. From Eq. 20 we identify the linear mean evoked response: Here, we find that the linear mean response is determined by only two time constants that shape the interplay between a low-and a high-pass. The amplitude of the linear response function for mean modulating signals in Eq. 21 is finite for any input frequency regardless of the stationary firing rate, see Fig. 4 (A). In the cross-over regime where the high-pass filter of the fixed threshold (nominator in Eq. 21) starts to counteract the low-pass filter of the membrane (denominator in Eq. 21) the response function transitions to a new firing rate dependent constant level. It is conceivable that the low-pass filtering by the membrane RC-circuit carries over to the firing rate dynamics, however, in this model even the highest frequencies can be relayed almost unattenuated. When is the linear approximation valid? As demonstrated in Fig. 4 and Fig. 5 linear response can be a good approximation for very small signal-to-threshold ratios of a few percent. In Fig. 6C we observe that for a signal-to-threshold ratio of A=y 0~0 :3 the linear approximation in Eq. 21 starts to be less accurate and non-linear effects kick. For A=y 0 §0:5 we observe that the linear approximation in Eq. 21 substantially overestimates the population response to hyperpolarizing (negative) transients and underestimates the response to depolarizing (positive) transients. Now, we address the population rate response to periodic modulation of the current variance. The membrane filtered voltage signal f (t) is then given by: Using s V (t) 2 , s _ V V (t) 2 and SV (t) _ V V (t)T as given in the Methods Section in Eq. 48,49 and 51 we obtain the population firing rate: The linear response function is: Here, the t s~ffi ffiffiffiffiffiffiffiffiffi t M t I p corresponds to the width of the voltage correlation function, t I is the current correlation time, and t M is ð23Þ the membrane time constant. Note, that Eq. 24 has been derived for an Ornstein Uhlenbeck current correlation of the form exp({abs(t)=t I ). We find that Eq. 24 is characterized by an interplay of two high and low-pass functions each with an individual time constant. In the mean modulating channel, on the other hand, we identified only one high-and one low-pass, and only two effective time constants. Fig. 4 (B) shows that the response for variance modulating signals is finite for any frequency, regardless of firing rate. A comparison between Fig. 4 (A) and Fig. 4 (B) reveals that mean and variance modulation can relay both slowly and fast varying signals.
Response to step-like input current changes. Abrupt changes in input current statistics can convey the onset of a sensory stimulus. Therefore, the time it takes for a population of neurons to alter its firing rate can impose limits of the detection and operation speed of cortical rate encoding [11]. It is conceivable that the membrane low-pass filtering could carry over to the firing rate dynamics and lead to a detection time scale in the order of the membrane time constant. Here we show, however, that the threshold model defies this intuition and can signal instantaneously the onset of mean-or variance-encoded step-like signals, even if they are subthreshold. Using the linear response function for mean modulations n 1 (v)~n m,1 (v) in Eq. 21 we obtain the following population rate change in response to step-like signals of amplitude A: Here, n(t) describes the population rate transient in response to a step-like increase of mean current (zA) at t~0. Eq. 25 and Fig. 5  (A,B) demonstrate that the response dynamics consists only of two components: an instantaneous component and an exponential governed by a single time scale t M . We recognize that a similar instantaneous component has been reported for Ornstein Uhlenbeck drive in the leaky integrate and fire model [10], but not for white noise drive [10,14,15]. For the white noise drive in the leaky as well as exponential integrate and fire model the response time scale is generally slower than the instantaneous component reported here. But their response time scale can, depending on input variance and firing rate, also be fast and substantially below the membrane time constant (see Fig. 1 F in [15]). As the stimulus amplitude increases, the linear approximation in Eq. 25 breaks down and needs to be replaced by the corresponding full response dynamics in Eq. 20. Fig. 6 schematically illustrates the range of amplitude strengths for which the non-linear effects kick in. Fig. 6B shows that already at an amplitude-to-threshold ratio of 0:3 sizable deviations of *0:5Hz~0:1n from the linear approximation are to be expected. Now, we come to the population rate dynamics evoked by variance changes. Using the linear response function for variance modulations n v,1 (v) (Eq. 24) we obtain the population firing rate transient in response to a step-like increase ds of the standard deviation: where t s~ffi ffiffiffiffiffiffiffiffiffi t I t M p . Fig. 5 shows the population rate transients in response to step-like changes in mean and variance predicted by Eqs. 25 and 26 alongside simulated results. Similarly to mean current steps, we find that for any n or t s the variance evoked response dynamics consists of an instantaneous component and exponential transient. The magnitude of the instantaneous component is increasing with firing rate but is largely unaffected by the correlation time t s . Thus, the response to mean and variance steps reported here is in line with the rapid response time observed in vitro cortical neurons [14,25] and also with the instantaneous response observed in the leaky integrate and fire model driven by colored noise [9,10,21].
Weak pairwise spike correlations. The common fluctuating current component shared by two neurons in a network can be viewed as a superposition of different frequencies. As such, it can be analyzed using the same tools as oscillatory or step-like signals studied in the previous chapters. Before we start, let us briefly review why we include pairwise spike correlations among the most crucial phenomena shaped by response functions. In cortical ensembles, pairwise spike correlations are known to play an important role in influencing population encoding strategies [40] and even predicting multi-neuronal firing patterns at larger distances [41]. Here, we focus on linear, weak pairwise spike correlations of two neurons firing at the same rate n. For simplicity, we assume the same statistical structure for the common fluctuating component n c (t) and the individual noise components n i (t) which make up the total input currents I 1 (t) and I 2 (t) in the two neurons. To account for the weak input correlations [5,31,32] we assume a weak correlation strength r, r%1, between the two input currents: The current correlation function C I (t) and the corresponding voltage correlation function C(t) are the same for both noise components. We characterize the spike correlations via the conditional firing rate n cond (t): n cond (t)~Ss 1 (t)s 2 (tzt)T=n: Because the contribution of n C (t) to the overall spiking activity is small and we can express it as a convolution of the linear impulse response function n m,1 (t) as in Eq. 15 [34]: The key ingredient in the calculation above is the property of the Fourier transforms of temporal derivatives, which for any function C(t) are F (C n (v))~(iv) n C(v). The peak value of this pairwise spike correlation is given by We find that peak spike correlations increase with firing rate (Fig. 7 A) and decrease with increasing time constant t s . Furthermore, the firing rate dependence of spike correlations captured by Eq. 33 corresponds to the firing rate dependent increase reported in the leaky integrate and fire model driven by white and correlated noise [42,43]. Furthermore, we can relate the correlated activity in Eq. 32 to the response dynamics evoked by oscillatory and step-like stimuli that we presented in previous sections. As we can see in Eq. 31 the firing rate dependent high-pass filter in Eq. 21 contributed the two terms proportional to C(t) and C''(t) to the pairwise correlation function in Eq. 32. Therefore, we can conclude that the firing rate dependence of linear mean evoked response function n 1,m (v) is directly related to the firing rate dependent shape of the pairwise spike correlations, in particular their firing rate dependent peak height and width.
Spike-triggering events. The dynamical response explored in the previous chapters results from the collective spike decisions of many neurons. As such it is intimately linked to the spike times and the spike triggering events on the level of single neurons. Here we explore the link between the dynamical population response and the voltage events which lead to the spike decision in individual neurons. Let us first formally define the spike triggered average voltage STA(t), a time lag t before a spike: Considering the spike train s(t)~Ð ?
{?n n m,1 (t{s)I(s)ds and the Fourier transform of the voltage correlation function C(v) we obtain: ð26Þ The key ingredient in the calculation above is the property of the Fourier transforms of temporal derivatives, which for any function C(t) are F (C n (v))~(iv) n C(v). An alternative derivation of this result via the Gaussian integrals is also shown in the Methods Section ''Spike triggered average voltage''. Fig. 8A demonstrates STA(t) as a function of firing rate n. As expected, we find in Fig. 8 that the spike triggered average voltage is mostly increasing towards the spiking threshold as the time to the spike is reduced. We also recognize that the spike triggered average exhibits a firing rate dependent transient hyperpolarization which is more pronounced for higher firing rates. We also note that the increasingly pronounced hyperpolarizing transient emerges even in the absence of any oscillatory component in the voltage correlation function. Notably, the rate dependent hyperpolarizing transients prior to spikes have been previously observed in cortical neurons [44]. The more pronounced hyperpolarizing transient has been interpreted as the possibility that properly timed hyperpolarizing events may increase the firing probability, potentially through a reduction of sodium channel inactivation or spike frequency adaptation [44]. We find here, that the transient hyperpolarization is mediated by the second term in Eq. 37 that is proportional to C'(t). It originates from the high-pass contribution of the response function n m,1 (v Taking this into account, we find for large t that STA(t)!exp({abs(t)=t M ) which is consistent with the integrate and fire result.
The spike triggered voltage covariance (STC(t)) is an additional popular and easily accessible measure for the characterization of a neuronal model or live neurons in vitro [44,46,47]. It is related to the second order Wiener (e.g. see p.330 in [48]) and could therefore serve well in future model-model or neuron-model comparisons. Using calculations similar to those of Burak and coworkers [36] and standard Gaussian integrals we obtain STC in the threshold model: If t 1~t2 we find the spike triggered variance STV (t)~s 2

Discussion
Here, we examined the relation between the frequency response functions in the mean and variance channels and the physiological parameters of single neurons, such as the functional form of the input current correlation, firing rate and membrane time constant. The threshold-based single neuron model we considered belongs to the class of spiking models that initiate a spike instantaneously after the threshold voltage is crossed. Instead of using the Fokker Planck framework we obtained the frequency response function via a direct integration of Gaussian probability densities that were modified to accommodate the mean or variance signals. This allowed us to systematically quantify the complete as well as linear frequency response in both channels, along with a number of related quantities such as population response to step signals, pairwise spike correlation function and spike triggered average. We confirmed all analytical results in numerical simulations of corresponding spiking neurons.

Frequency response functions in the mean and variance channels
We derived the complete as well as linear frequency response function for mean modulating signals (Eq. 20,21). As can be expected from models with an instantaneous spike generation, the frequency response functions had a finite high-frequency limit. In the linear regime the mean evoked frequency response could be reduced to an interplay between a low-pass and a high-pass filter which were governed by only two independent time constants: the membrane time constant t M and an effective threshold-dependent time constant. Furthermore, both linear and non-linear response functions did not depend on the specific functional form of the voltage correlation function C(t), only its temporal width t s was important for the dynamical response. Notably, the population response to a mean current step could be described by an explicit tractable expression consisting of two components: the instantaneous component governing the immediate response and an exponential transient described by the membrane time constant (Eq. 25). We also found that the linear step-response function can be a good approximation for the population rate response if the signal-to-threshold ratio is below 0:1.
For variance modulating signals we were also able to provide the complete as well as linear frequency response function (Eq. 23,24). We observed that the frequency response function remained finite in the high-frequency limit. In the linear regime, the variance evoked response could be described by an interplay between two high-and low pass functions. As in the mean channel, we found that population response exhibited two components: the instantaneous component that occurs immediately after the step onset and a combination of exponential functions that were governed by three time constants (Eq. 26). For the pairwise spike correlations that were obtained using the mean evoked linear response, we observed that the spike correlation peak increased and the temporal width decreased with firing rate (Fig. 7). The spike triggered average voltage and current could be described by an explicit expression as a sum involving the correlation function and its second derivative. Here, we found that with increasing firing rate the spike triggered average exhibited an increasingly pronounced undershoot shortly before a spike (Fig. 8 A).

Model limits
The model framework that we used here to describe the spiking activity of each neuron in a cortical population is based on three major assumptions. The first is the confinement to the fluctuationdriven regime, second the Gaussian voltage statistics [49] and third the assumption that spikes are instantaneously generated upon the crossing of a fixed threshold voltage. The confinement of this model to the fluctuation-driven regime and Gaussian voltage statistics is in line with the recent experimental evidence that the fluctuation rather than mean depolarization driven regime is the primary operation scheme in cortical neurons. The first line of support is the remarkable cortical balance of excitation and inhibition such that neuronal firing is driven by fluctuations that transiently escape this cancellation [29,30]. The second line of evidence are the exceptionally low firing rates v1Hz [5,6]. Mean and fluctuation driven regimes can differ significantly in their spike train regularity, yet they seem to exhibit very similar spike cross correlation [42] and dynamical response properties [23]. Therefore, numerous features of dynamical population response could already be understood by studying only the fluctuation-driven regime. Even though skewed input current distributions that violate the Gaussian assumption can be of interest in specific cases, experimental evidence suggests that the Gaussian distribution can be a good match for in vivo background fluctuations in the prevalent cortical cell type of pyramidal neurons, e.g see Box 1 (top) in [26]. The third assumption, that spikes are instantaneously generated upon the crossing of a fixed threshold voltage, is motivated by the observation that time scales of spike onset can be very short, e.g. the onset rise time *1{2ms and the spike slope factor that is proportional to the radius of the curvature of the I-V curve at its minimum is *1mV [50,51]. As predicted by the previous models, the instantaneous threshold-based spike generation assumed in our model leads to a non-decaying frequency limit in the response functions [11,19,20]. These studies indicated that the frequency response function for low and intermediate frequency range are largely unaffected by the spike onset time and their properties can be could be explored in a model with an instantaneous spike generation mechanism. Alternative neuron models with more involved spike generation mechanisms such as quadratic or exponential integrate and fire indicate the possibility that the frequency response functions in a pair of neurons can depend on the model specifics [15,52]. Yet, realizing how remarkably accurate many cortical spike synchrony features and response dynamics can be modeled by a clearly barebonethreshold model, we are convinced that this model will find its place alongside the classical integrate and fire models and offer a valuable maximally tractable limiting case for future studies.

Comparison with previous theoretical studies
The linear and non-linear response functions are basic tools for the description of any physical system. In theoretical neuroscience several earlier studies have quantified the response functions in various model types, particularly in the integrate and fire models. To our knowledge, our study is the first to provide explicit expressions for the complete linear and non-linear response functions for both mean and variance channels in the presence of correlated noise, and also to derived from them tractable expressions for the step-evoked population response function, the full pairwise spike correlation function and the spike triggered averages. In the mean channel we showed for the first time that the frequency response function and all derived quantities can be independent of the functional form of the input correlations. The only important parameter determining in our formalism the frequency response function, step response dynamics or correlation gain was the width of the correlation function, but not its functional specifics.
Brunel and colleagues showed analytically in the integrate and fire model driven by the correlated Ornstein Uhlenbeck current that the linear frequency response function has a finite highfrequency limit [10]. Furthermore, the authors obtained numerically the functional form of the linear response function and the corresponding step-evoked transient and observed that two factors increase the high-frequency response and the instantaneous component of the step response: the firing rate and the temporal width of the correlation function. In our threshold model, we observe the same functional dependence of the linear response function on firing rate and temporal width. Here, however, we could show that this behavior is not unique to the Ornstein Uhlenbeck current. In fact, these effects can be observed for a wide range of colored input because the functional form of input correlation does not influence the frequency response function. This is important news, because it indicates that results obtained for the Ornstein Uhlenbeck current in the integrate and fire model [10,15,21,52] could also be valid for more general current correlation functions. Let us also note, that the instantaneous response component in the leaky integrate and fire model vanishes for the physiologically remote choice of white noise drive (see Fig. 3 (top left) in [10]). In this case, the response dynamics is a superposition of different time constants [10,15]. Therefore, the presence or absence of the instantaneous response component is model specific and some models such as the exponential integrate and fire lack an instantaneous component, but their response time is finite and can be faster than the membrane time constant (Fig. 1 F in [15]).
In the variance channel, Lindner and colleagues showed that the leaky integrate and fire model supports a faithful transmission of high-frequency inputs for white noise drive, see Eqs. 4-6 [9]. Also a number of follow-up studies confirmed this in the leaky, quadratic and exponential integrate and fire models [9,12,14,19,21]. However, these results focused on white noise background and provided the linear response function for colored noise only in specific limits, such as the infinitesimal or infinite input frequency limit [19] and did not address the full linear or non-linear response functions. Therefore, Eqs. 23 and 24 are, to our knowledge, the first to describe the complete as well as linear frequency response function in the presence of correlated noise. Also for the step response dynamics in the variance channel, only the slope of the initial response phase could be obtained so far, e.g. see Eq. 28 in [19]. Here, however, we provide the complete step-evoked response function Eq. 26 that can be generalized to the nonlinear response regime using Eq. 24.

Relation to experimental data
Starting in 1970s a number of studies have addressed the frequency response functions in cortical neurons. The majority of these studies [23,53,54] were conducted in the mean channel under a variety of different noise conditions. Bair and colleagues have shown that neurons in the medial-temporal (MT) area of the monkey cortex reliably transmitted input frequencies in the range 30{100Hz [55]. Also, three subsequent experimental studies in more controlled in vitro conditions have demonstrated that the reliably encoded frequency range can be tens of times larger than the firing rate of individual neurons. For a firing rate of *5Hz this can mean a reliably encoded frequency rage of up to 200{300Hz [23][24][25]. This is consistent with broad reliably encoded frequency range of the threshold model presented here, as well as the integrate and fire model for colored noise. Notably, the class of models with threshold-based instantaneous spike initiation are so far the only model types that enable reliable transmission of inputs that are much higher than the firing rate of individual neurons [11,19,20]. A key, prediction of our model framework is that response to high frequencies in the mean channel can be enhanced by a higher correlation time of the noise background. Notably, this correlation time dependence of the response function has been shown in vitro [25]. Here, we conclude that essential properties of the mean evoked frequency response function can be understood within the results derived from the threshold model. In fact, this model goes a step further and predicts that these properties observed in vitro for the Ornstein Uhlenbeck drive could be generalized to other type of input correlations.
In the time domain, the threshold model predicts an instantaneous response to step-like stimuli. Yet surprisingly, early experimental evidence presented by Silberberg and colleagues did not support the presence of an instantaneous response component [14]. However, this study did not quantify the response time scale and its conclusions might be biased by the use of a 4AP induced noise of unspecified correlation time or an almost white noise background, which in theoretical studies has been shown to lack an instantaneous component. The lack of rapid population response to step input, on the one side, and remarkably broad range of encoded frequencies on the other side presented an apparent contradiction. A recent study, however, has confirmed the broad reliably encoded frequency range in vitro and in vivo and has shown that in the presence of correlated noise cortical neurons in vitro can detect the step onset within milliseconds after its onset for a variety of conditions [25]. Also, a recent in vivo study by London and colleagues reported remarkably fast detection of subtle 25pA mean current pulses into a single neuron in the local cortical circuit in vivo [56]. This rapid response onset and the broad reliably encoded frequency range in vitro and in vivo that depends on the correlation time of the background can be understood using the theoretical results obtained here in a threshold model framework. Let us also note that the pairwise spike synchrony properties derived from the mean response function in the threshold framework, such as the firing rate dependence of correlation gain, are consistent with the experimental observations in vivo and in vitro [5,38,42]. For variance-encoded signals, the only experimentally study by Boucsein and colleagues reported a remarkably broad range of encoded frequencies in vitro [24]. The ability to rapidly detect the step-like change of input current variance in a population of neurons has, however, been shown only for large changes of variance (for doubling [14] or tripling [25] of the input current variance). In summary, while the existence of a fast response component in response to step-like variance changes could be confirmed in vitro it may require a larger amplitude than predicted by this threshold model. This could be due to additional cellular threshold adaptation mechanisms or voltage dependent conductance changes-effects that are not included in this threshold model framework.

Outlook and novel model predictions
Here, we have shown that an alternative threshold model framework can capture many important features of the frequency response functions in vivo and in vitro. Importantly, this framework also offers novel predictions that can be tested experimentally. One important new prediction is the independence of the mean evoked response function on the functional form of the input correlations. In Eq. 21 the mean evoked response function depends only on the width t s of the input correlation in the vicinity of zero, but not on any other functional specifics. For live neurons, this prediction means that any change in ion channel composition that affects the periphery rather than the central peak of the correlation function will not affect the frequency response function. This prediction can be directly tested in vitro with two current correlation functions which share the same t s but have otherwise very different functional forms. On the theoretical side, we provided a number of tractable expressions for the complete spike triggered average and complete linear pairwise spike correlation function, and mean as well as variance evoked response functions. This analytical tractability framework will make it a useful addition to the theoretical toolkit of integrate and fire models and facilitate future model-model as well as model-experiment comparisons.
Taking s 2 V (t) and s 2 _ V V (t) allows us to calculate the firing rate C~s Solving this integral using standard Gaussian integrals and covariances obtained above we obtain Eq. 23.

Quantification of oscillatory firing rate modulations
We use the vector strength r [8] to obtain the linear response function numerically: where t i are the spike times, N their number, T~2p=v is the period length and a is the relative amplitude of the firing modulation evoked by the signal A exp(ivt). Taking a~Ajn 1 (v)j=n we can directly identify the vector strength as r~An 1 (v)=(2n). Here n 1 (v) represents either n m,1 (v) or n v,1 (v). Fig. 9 illustrates the relation between the linear response function n 1 (v) and the vector strength r. Because the vector strength r is constructed directly from the spike times, it omits a sinusoid fit of the peristimulus time histogram, which can potentially be biased by the fitting algorithm or the bin size chosen.

Spike triggered average voltage
Here, we obtain the spike triggered average voltage STA(t) a time lag t before a spike in the threshold model via a direct calculation of Gaussian integrals: 1 n where the correlation matrix is Using standard Gaussian integrals we obtain the result in Eq. 37, which reads: