Volterra dendritic stimulus processors and biophysical spike generators with intrinsic noise sources

We consider a class of neural circuit models with internal noise sources arising in sensory systems. The basic neuron model in these circuits consists of a dendritic stimulus processor (DSP) cascaded with a biophysical spike generator (BSG). The dendritic stimulus processor is modeled as a set of nonlinear operators that are assumed to have a Volterra series representation. Biophysical point neuron models, such as the Hodgkin-Huxley neuron, are used to model the spike generator. We address the question of how intrinsic noise sources affect the precision in encoding and decoding of sensory stimuli and the functional identification of its sensory circuits. We investigate two intrinsic noise sources arising (i) in the active dendritic trees underlying the DSPs, and (ii) in the ion channels of the BSGs. Noise in dendritic stimulus processing arises from a combined effect of variability in synaptic transmission and dendritic interactions. Channel noise arises in the BSGs due to the fluctuation of the number of the active ion channels. Using a stochastic differential equations formalism we show that encoding with a neuron model consisting of a nonlinear DSP cascaded with a BSG with intrinsic noise sources can be treated as generalized sampling with noisy measurements. For single-input multi-output neural circuit models with feedforward, feedback and cross-feedback DSPs cascaded with BSGs we theoretically analyze the effect of noise sources on stimulus decoding. Building on a key duality property, the effect of noise parameters on the precision of the functional identification of the complete neural circuit with DSP/BSG neuron models is given. We demonstrate through extensive simulations the effects of noise on encoding stimuli with circuits that include neuron models that are akin to those commonly seen in sensory systems, e.g., complex cells in V1.


INTRODUCTION
Intrinsic noise sources are diverse and appear on many levels of a neuronal system ranging from electrical to chemical noise sources (Faisal et al., 2008;Destexhe and Rudolph-Lilith, 2012) and from single cells to networks of neurons. At the cellular and subcellular level, variability in biochemical reactions leads to stochastic transduction processes (Song et al., 2012), and ion channel fluctuations (Neher and Sakmann, 1976;White et al., 1998) result in variability in spike generation and propagation (Faisal and Laughlin, 2007). At the network level, probabilistic quantal release of neurotransmitters (Katz, 1962), background synaptic activity Jocobson et al., 2005) and variability in timing of spikes from presynaptic neurons (Faisal and Neishabouri, 2014) are sources of stochastic fluctuation of synaptic conductances (Destexhe et al., 2001) that are believed to have a major impact on spike time variability (Yarom and Hounsgaard, 2011).
The existence of sources of noise also leads to variability in the spike times even when neurons are subject to the same, repeated inputs (Calvin and Stevens, 1968;Berry et al., 1997;de Ruyter van Steveninck et al., 1997). Spikes are the primary form of carriers of information in the nervous system and their timing is thought to be relevant to the message neurons need to convey (Rieke et al., 1999). Therefore, the variability of spike timing may reduce or damage the information being transmitted. It is quite remarkable, however, that sensory systems manage to be very robust even if they are subject to interference due to noise. Visual and auditory systems are two examples in which the stimuli are highly time varying. These systems have been reported to convey information with high spike timing precision (Butts et al., 2007;Kayser et al., 2010).
Noise may be useful in facilitating signal detection (McDonnell and Ward, 2011). Still, interference due to noise poses an important limit on how well sensory systems can represent input stimuli. It is not clear how intrinsic noise sources affect the representation of sensory inputs based on spike times, and how they impact the functional identification of sensory neurons.
We study the representation of sensory stimuli using a novel neural circuit model, that extends previously proposed models Lazar and Slutskiy, 2014, in press) in terms of architectural complexity and the existence of intrinsic noise sources. Our base level circuit architecture consists of two interconnected neurons, each with two cascaded stages. The first stage comprises two types of dendritic stimulus processors. The first dendritic stimulus processor performs nonlinear processing of input stimuli in the feedforward path leading to the spike generator. The second dendritic stimulus processor performs nonlinear processing in the feedback loop whose inputs are spike trains generated by biophysical spike generators (BSGs). The BSGs constitute the second stage of the base level circuit.
Our nonlinear dendritic stimulus processors describe functional I/O relationships between the dendritic outputs in the first stage and inputs that are either sensory stimuli or spikes generated by BSGs. DSPs are modeled using Volterra series. Volterra series have been used for analyzing nonlinear neuronal responses in many contexts (Lu et al., 2011;Eikenberry and Marmarelis, 2012), and have been applied to the identification of single neurons in many of sensory areas (Benardete and Kaplan, 1997;Theunissen et al., 2000;Clark et al., 2011). Volterra dendritic processors can model a wide range of nonlinear effects commonly seen in sensory systems (Lazar and Slutskiy, in press).
Here, in addition, we introduce nonlinear interactions between neurons in the feedback and cross-feedback paths. This gives rise to interesting neural processing capabilities directly in the spike domain, e.g., coincidence detection (Agmon-Snir et al., 1998;Stuart and Häusser, 2001). The relationships described here by the Volterra model are functional and do not address the underlying circuit/dendritic tree level interactions. However, the latter have recently been subject to intense investigations (London and Häusser, 2005;Wohrer and Kornprobst, 2009;Werblin, 2011;Xu et al., 2012;Yonehara et al., 2013;Zhang et al., 2013). Conductance-based, biophysical spike generators are well established models that have been extensively used in studies of neuronal excitability and in large simulations of spiking neural networks (Izhikevich, 2007). Following Lazar (2010), we use formal BSG models to represent sensory stimuli under noisy conditions.
We formulate the encoding, decoding and functional identification problems under the neural encoding framework of Time Encoding Machines (TEMs). In this modeling framework the exact timing of spikes is considered to carry information about input stimuli (Lazar and Tóth, 2004). The separation into dendritic stimulus processors and spike mechanisms mentioned above allows us to study synaptic inputs and spike generation mechanisms separately, and hence independently model the intrinsic noise sources of each component. We incorporate two important noise sources into a general single-input multi-output neural circuit model. The first is a channel noise source that arises in spike generation (White et al., 2000). The second is a synaptic noise source due to a variety of fluctuating synaptic currents (Manwani and Koch, 1999).
Based on the rigorous formalism of TEMs, we show how noise arising in dendritic stimulus processors and in biophysical spike generators is related to the measurement error in generalized sampling. Dendritic stimulus processing and spike generation can then be viewed as a generalized sampling scheme that neurons utilize to represent sensory inputs . Contrary to traditional sampling where the signal amplitude is sampled at clock times, neurons asynchronously sample all stimuli.
We systematically investigate how the strength of noise sources degrades the faithfulness of stimulus representation and the quality of functional identification of our proposed class of neural circuits. Furthermore, since the representation is based on spike timing, it is natural to investigate how spike timing variability affects the precision in representing the amplitude information of sensory stimuli.
The work presented here requires a substantial amount of investment in the mathematical formalism employed throughout. There are a number of benefits in doing so, however. Formulating the problem of stimulus encoding with a neural circuit with intrinsic noise sources as one of generalized sampling, i.e., of taking noisy measurements is of interest to both experimentalists and theoreticians alike. Understanding that the problem of neural decoding and functional identification are dual to each other is key to building on either or both. Finding how many repeat experiments need to be performed for a precise quantitative identification of Volterra kernels is of great value in neurophysiology. A further qualitative insight of our work is that for neural circuits with arbitrary connectivity, feedforward kernels are typically easier to estimate than feedback kernels. Finally, our finding that some key nonlinear neural circuits are tractable for detailed noise analysis suggests a wide reaching analytical methodology.

MODELING NONLINEAR NEURAL CIRCUITS, STIMULI, AND NOISE
We present in Section 2.1 the general architecture of the neural circuits considered in this paper. In Section 2.2 we discuss the modeling of the space of stimuli. Volterra DSPs are the object of Section 2.3. Finally, in Section 2.4 we provide models of BSGs with intrinsic noise sources.

NEURAL CIRCUIT ARCHITECTURE
The general architecture of the neural circuit considered here is shown in simplified form in Figure 1. It consists of two neurons with a common time-varying input stimulus. With added notational complexity the neural circuit in Figure 1 can easily be extended in two ways. First, multiples of such circuits can encode a stimulus in parallel (see Section 2.1 in the Supplementary Material). In this case only pairs of neurons are interconnected through the feedback kernels. Second, more neurons can be considered in the neural circuit of Figure 1; all these neurons can be fully interconnected through feedback loops. Each neuron i, i = 1, 2, receives a single time-varying input stimulus u 1 (t). The modeling of the input stimulus is discussed in Section 2.2. The output of each of the biophysical spike generators (BSGs) is a spike sequence denoted by (t 1 k ) and (t 2 l ), k, l ∈ Z. The input stimulus u 1 (t) is first processed by a feedforward Dendritic Stimulus Processor (feedforward DSP) (Lazar and Slutskiy, in press). The feedforward DSP models the aggregated effect of processing in the neural circuits in the prior stages and in the dendritic tree of neuron i = 1, 2. For example, if the neurons in the model circuit are considered to be Retinal Ganglion Cells (RGCs), then the feedforward Volterra DSP models the processing that takes place in the outer-and inner-plexiform layers of the retina as well as in the dendritic trees of an RGC (Werblin, 2011;Masland, 2012). The feedforward DSPs are modeled here as second order Volterra expansion terms (Volterra, 1930). The first order terms h 11i 1 (t) in the feedforward DSPs are linear filters typically used in modeling receptive fields. The second order terms h 11i 2 (t 1 , t 2 ) model nonlinear operations on the stimulus u 1 (t).
A second group of Volterra DSPs models the cross-feedback interactions between the two neurons. Instead of time-varying stimuli, the output spikes generated by the BSGs are the inputs to these DSPs. We therefore refer to these as feedback Dendritic Stimulus Processors (feedback DSPs). The output spikes of each individual neuron i are processed by the first order term h 2ji 1 (t), i, j = 1, 2, i = j. In addition, output spikes from both neurons interact nonlinearly through the second order terms h 2ji 2 (t 1 , t 2 ), i, j = 1, 2, i = j. The summed responses from the first order feedback DSP h 2ji 1 and the second order feedback DSP h 2ji 2 are fed back to neuron i as additional dendritic currents.
The dendritic currents consisting of the output of the DSPs with added noise are subsequently encoded by biophysical spike generators. BSGs are biophysically realistic axon hillock spike generator models that are governed by a set of differential equations with multiple types of ion channels (Hodgkin and Huxley, 1952;Izhikevich, 2007). The detailed BSG models are introduced in Section 2.4. The spike times of output spikes generated by the BSGs are assumed to be observable.
We identify two intrinsic noise sources of the proposed neural circuit. First, the feedforward DSPs and the feedback DSPs are affected by additive Gaussian white noise. This noise arises from the combined effect along the path from sensory transduction to synaptic integration and includes synaptic background noise and stochasticity in the dendritic tree (Manwani and Koch, 1999;Fellous et al., 2003;Destexhe and Rudolph-Lilith, 2012). Since the outputs of the feedforward and feedback DSPs are additively combined, we consider, for simplicity, a single source of additive Gaussian white noise. Second, the ion channels of the BSGs are intrinsically stochastic and introduce noise in the spike generators (White et al., 2000;Hille, 2001).

MODELING SIGNAL SPACES
Two signal spaces will be considered here. The first, models the space of input signals to feedforward DSPs. The second models the space of input spikes to feedback DSPs. These spaces will be formally described below.

Modeling the space of input stimuli
We model the space of input stimuli as a Reproducing Kernel Hilbert Space (RKHS) (Berlinet and Thomas-Agnan, 2004). RKHSs are versatile vector spaces for modeling signals arising in computational neuroscience, signal processing and machine learning. For example, auditory signals, olfactory signals and visual signals can readily be modeled as band-limited functions of an RKHS with a sinc or Dirichlet kernel Lazar and Slutskiy, 2013). A particular choice of RKHSs in this article is the space of trigonometric polynomials. The computational advantage of working on the space of trignometric polynomials has been discussed  and is closely related to the algorithmic tractability of the Fourier series in the digital domain. If the biological signals have unknown bandwidth with a spectrum that falls off fast enough, many Sobolev spaces might be a suitable choice of RKHS (Berlinet and Thomas-Agnan, 2004;Lazar and Pnevmatikakis, 2009). In such spaces the norm may include the derivative of the signal, i.e., the rate of change of the signal that many neurons are sensitive to Kim et al. (2011).
The space of trigonometric polynomials is defined as below.
Definition 2.1. The space of trigonometric polynomials H 1 1 is a function space whose elements are functions defined on the domain where are a set of orthonormal basis functions. 1 denotes the bandwidth and L 1 is the order of the space.
H 1 1 endowed with the inner product: is a Hilbert Space. Intuitively, the basis functions e l (t), l = −L 1 , . . . , L 1 , can be interpreted as a set of discrete spectral lines uniformly spaced in the frequency domain between − 1 and 1 . For a given signal u 1 (t), the amplitude of its spectral lines is determined by the coefficients u l , l = −L 1 , . . . , L 1 .

Remark 2.2.
Functions in H 1 1 are periodic over R with period S 1 = 2πL 1 1 . Therefore, the domain D 1 covers exactly one period of the function. Note that the u l 's are closely related to the Fourier coefficients of the periodic signal u 1 (t), and can thereby be very efficiently computed via the Fast Fourier Transform.
H 1 1 is an RKHS with reproducing kernel (RK) It can be easily verified that the RK satisfies the reproducing property Definition 2.3. We shall also consider the tensor product space H 1 2 on the domain D 2 = [0, S 1 ] × [0, S 1 ], whose elements are of the form where e l 1 l 2 (t 1 , are a set of functions forming an orthonormal basis.
H 1 2 is again an RKHS with RK Note that we use the subscript to indicate the dimension of the domain of functions, i.e., the number of variables the functions in the RKHS have, and use the superscript 1 to indicate the input space.

Modeling the space of spikes
The feedback kernels of the neural circuit in Figure 1 receive as inputs spike trains generated by the BSGs. Spike trains are often modeled as sequences of Dirac delta pulses and, consequently, the outputs of linear feedback kernels are the result of superposition of their impulse responses (Keat et al., 2001;Pillow et al., 2008;Lazar et al., 2010). Dirac delta pulses have infinite bandwidth. Spikes generated by the BSGs, however, have limited effective bandwidth. Following (Lazar and Slutskiy, 2014) spikes are modeled to be the RK of an one-dimensional Hilbert space H 2 1 at spike time occurrence. Here H 2 1 is a space of trigonometric polynomials whose order L 2 , period S 2 and bandwidth 2 may differ from the input stimulus space H 1 1 , where 2 shall be larger than the bandwidth assumed for the feedback kernel, and S 2 is much larger than the support of the feedback kernel (Lazar and Slutskiy, 2014). A spike at time t i k of neuron i can then be expressed in functional form as K 2 1 (t i k ; t), where the superscript indicates that the RK belongs to the spike input space.
Due to the reproducing property, single or pairs of input spikes have the property The operator P 2 is similarly defined to P 1 above; it denotes, however, the projection onto the space of spikes. Thus, not surprisingly, incoming spikes directly readout the projection of the feedback kernels. By letting L 2 → ∞, in L 2 norm as the RK converges to the sinc function and the RKHS becomes the space of band-limited signals . A more detailed analysis is available in Lazar and Slutskiy (2014). This formalism will be employed for solving the functional identification problem formulated in Section 4.1.

VOLTERRA DENDRITIC STIMULUS PROCESSORS
As mentioned in Section 2.1, two forms of dendritic stimulus processing appear in our model.

Feedforward Volterra dendritic stimulus processors
The feedforward DSPs are modeled as up to second order terms in the Volterra series. The feedforward DSPs take continuous signals in the stimulus space as inputs, while the output can be expressed as (see also Figure 1) where h 11i 1 ∈ L 1 (D 1 ) and h 11i 2 ∈ L 1 (D 2 ) denote, respectively, the first and second order Volterra kernels, i = 1, 2. They are assumed to be real, causal and bounded-input bounded-output (BIBO)-stable. It is also assumed that both h 11i 1 and h 11i 2 have finite memory. In addition, h 11i 2 is assumed, without loss of generality, to be symmetric, i.e., h 11i 2 (t 1 , t 2 ) = h 11i 2 (t 2 , t 1 ).
Example 2.5. We present here a Volterra DSP that is akin to a model of dendritic stimulus processing of complex cells in the primary visual cortex (V1). The difference is that the complex cells operate spatio-temporally, whereas in the example given below they operate temporally. We first consider two first order kernels based on Gabor functions, g c (t) = exp − (t − 0.13) 2 2 · 0.0005 cos 2π · 10 · (t − 0.13) , g s (t) = exp − (t − 0.13) 2 2 · 0.0005 sin 2π · 10 · (t − 0.13) .
The two filters are Gaussian modulated sinusoids, that are typically used to model receptive fields of simple cells in the primary visual cortex (V1) where the variables denote space instead of time (Lee, 1996;Dayan and Abbott, 2001). In addition, the two filters are quadrature pair in phase. Both filters are illustrated in Figure 2A.
The response of applying the input stimulus u 1 on the temporal filters with impulse response g c and g s is given by D 1 g c (t − s)u 1 (s)ds and D 1 g s (t − s)u 1 (s)ds, respectively. The responses of the two linear filters of the complex cell model are squared and summed to produce the phase invariant measure v i (Carandini et al., 2005) where h 11i 2 (t 1 , t 2 ) = g c (t 1 )g c (t 2 ) + g s (t 1 )g s (t 2 ). Therefore, the operation performed by a complex cell can be modeled with a second order Volterra kernel. h 11i 2 is shown in Figure 2B. We now take a closer look at the operation of the second order kernel. The two dimensional convolution of the second order kernel with u 2 (t 1 , t 2 ) is shown in Figure 2C.
It is important to note that, since the second order kernel has finite memory, it may not have enough support to cover the entire domain D 2 for u 2 (t 1 , t 2 ). For example, as illustrated in Figure 2C, the output of the second order feedforward DSP at time t is given by the integral of the product of u 2 (t 1 , t 2 ) and a rotated h 11i 2 with the origin shifted to (t, t) [see also (10)]. Since the shift is along the diagonal, only u 2 (t 1 , t 2 ) in the domain that is contained within the black lines is multipled by nonzero values of h 11i 2 . u 2 (t 1 , t 2 ) elsewhere in the domain is always multiplied by zero in evaluating the output. Therefore, the output of the second order filter only contains information about u 2 within the domain located in between the black lines in Figure 2C. This has implications on decoding the signal (see also Remark 3.11 in Section 3.2)

Feedback Volterra dendritic stimulus processors
As already mentioned, the feedback DSPs do not operate on stimuli directly but rather on spikes generated by BSGs. We assume that h 2ji 1 ∈ L 1 (D 1 ), h 2ji 2 ∈ L 1 (D 2 ), i = j, are real, causal, BIBOstable and have finite memory. In addition, we assume that these kernels are effectively band-limited (see also Section 2.2.2). In functional form we denote a train of spikes as k K 2 1 (t i k ; t). The output of the feedback DSP i amounts to with j = i. 2 is only shown in a restricted domain and is zero elsewhere). For t = 0.3, the output of the convolution is the integral of the product of the rotated Volterra kernel and the signal underneath. Since the convolution is evaluated on the diagonal t = t 1 = t 2 , the second order kernel shifts, as t increases, along the arrow on the diagonal. See also Supplementary Figure 5E.
In particular, the inputs to the second order term of the feedback DSPs are generated by two neurons. This allows for modeling nonlinear interactions between the two neurons in the spike domain.

Overall output from DSPs
The overall inputs (without noise) to the two BSGs in Figure 1 are The system of Equations (12) above functionally describe the post-synaptic aggregate currents that are injected into the BSG i.
There are a variety of noise sources to be considered. Synaptic variability of feedforward DSPs adds noise sources to the current input to the BSGs. These include thermal noise, synaptic background noise, etc. (Jonston, 1927;Calvin and Stevens, 1968;Manwani and Koch, 1999;Fellous et al., 2003;Destexhe and Rudolph-Lilith, 2012). Feedback DSP kernels may themselves be subject to intrinsic noise sources that may lead to variability in the spike generation process. Intrinsic variability of BSG spike times can, e.g., contribute to the variability of the aggregate current driving the axon hillock in feedback loops.
Overall, the combined effect of DSP noise sources is modeled as Gaussian white noise processes that are added to the feedforward and feedback DSP outputs. The sum total of signal and noise represents the aggregate current input to the BSGs (see Figure 1). Formal DSP noise models will be incorporated directly into the BSG model presented in the next section.

BSGs and phase response curves
We consider biophysically realistic spike generators such as the Hodgkin-Huxley, Morris-Lecar, Connor-Stevens neurons (Hodgkin and Huxley, 1952;Connor and Stevens, 1971; Morris and Lecar, 1981). The class of BSGs can be expressed in vector notation as where x i are the state variables, f i are vector functions of the same dimension, and I i are the constant bias currents in the voltage equation of each BSG. Each input current v i (t) is applied to the neuron i by additive coupling to the voltage equation, typically the first of the set of ordinary differential equations, i.e., where 0 is a row vector of appropriate size. We assume that the neuron is periodically spiking when no external input is applied. This can be satisfied by a constant bias current I i additively coupled onto the voltage equation. The use of I i is necessary to formulate the encoding for the single neuron case, and this assumption will be relaxed later in this article.
A large enough bias current induces a periodic oscillation of the biophysical spike generator. Therefore, the phase response curve (PRC) is well defined for this limit cycle (Izhikevich, 2007). We denote the PRC of the limit cycle induced by the bias current

Channel noise in BSGs
As shown in Figure 1, we consider BSGs with noise sources in the ion channels. The noise arises due to thermal fluctuations (White et al., 2000;Hille, 2001) as the finite number of ion channels in the BSGs open and close stochastically. The differential equations that govern the dynamics of the BSGs in (14) are deterministic. The set of stochastic differential equations (SDEs) below represent their stochastic counterpart (Lazar, 2010): where B i is a matrix with state dependent values, , and W i p (t), p = 2, · · · , P i , are independent Brownian motion processes. Note that P i does not necessarily have to be equal to N i , the number of state variables. The first element in the stochastic differential dZ i is the aggregate dendritic input v i dt driving the voltage equation. The other entries in dZ i are noise terms that reflect the stochastic fluctuation in the ion channels / gating variables.
Randomness is often added to BSGs by setting B i = I, where I is a N i × N i identity matrix. The later setting can be viewed as adding subunit noise . Recently, it has been suggested that a different way of adding channel noise into the BSGs may result in more accurate stochastic behavior Linaro et al., 2011;Orio and Soudry, 2012). The SDEs in (15) are of general form and do not preclude them. In fact, by setting B i to be a block matrix with blocks equal to be the square root of the diffusion matrix for each ion channel, the channel SDE model Orio and Soudry, 2012) can easily be incorporated into (15).
Finally, we note that, under appropriate technical conditions the SDE formulation applies to BSGs with voltage-gated ion channels as well as other types of ion channels. The conditions require that the BSG model can be treated mathematically as a system of SDEs of the form (15) and that the latter satisfies the assumptions of Section 2.4.1.

Overall encoding of the neural circuit model
Taking into account the dendritic input from the feedforward DSPs and feedback DSPs, the encoding by the neural circuit model under the two noise sources is given by two systems of SDEs. With the Brownian motion W i 1 modeling the DSP white noise, the encoding of neuron i, i = 1, 2, can be expressed as where with v i (t) given by Equation (12).
Note that in the system of Equations (16) the two output spikes trains t i k , i = 1, 2, k ∈ Z, are the observables. Due to the intrinsic noise sources in the DSPs and in the BSGs, spike timing jitter may be observed from trial to trial by repeatedly applying the same stimulus to the neural circuit (see Section 2.3 in the Supplementary Material).

ENCODING, DECODING, AND NOISE
In Section 3.1 we present the mathematical encoding formalism underlying the neural circuit in Figure 1. We formulate stimulus decoding as a smoothing spline optimization problem and derive an algorithm that reconstructs the encoded signal in Section 3.2. Finally, we analyze the effect of noise on stimulus decoding in Section 3.3.

ENCODING
In this section, we formulate a rigorous stimulus encoding model based on the neural circuit shown in Figure 1. The input of the circuit is a signal u 1 modeling a typical sensory stimulus as described in Section 2.2.1. The neural circuit generates a multidimensional spike train that is assumed to be observable. We establish model equations by first describing the I/O relationship (i.e., the t-transform) of a single BSG. We then provide the t-transform of the entire neural circuit model that maps the input stimulus amplitude into a multidimensional spike timing sequence.

The I/O of the BSG
In the presence of a bias current I i and absence of external inputs, each BSG in Figure 1 is assumed to be periodically spiking. Provided that the inputs are small enough, and by using the PRC, the BSG dynamics of spike generation can be described in an one-dimensional phase space (Lazar, 2010).
Definition 3.1. A neuron whose spike times (t i k ), k ∈ Z, i = 1, 2, verify the system of equations where with τ i 0, I i = 0 and x i t, I i the periodic solution to (13)  As its name suggests, the PIF projects a weighted version of the input embedded in noise and the ion channel noise associated with the gating variables (B i dZ i ) onto the PRCs of the corresponding gating variables on a time interval between two consecutive spikes. Note that the integrand in (17) is identical to the RHS of (19). τ i (t, I i ) on the LHS of (19) denotes the phase deviation and is driven by the perturbation on the RHS. The LHS of (17) represents the phase deviation measurement performed by the PIF neuron. The RHS of (17) provides the value of the measurement and is equal to the difference between the inter-spike interval and the period of the limit cycle.
The BSG and the PIF neuron with random thresholds are, to the first order, I/O equivalent (Lazar, 2010). In Lazar (2010) it was also shown that a good approximation to the PIF neuron is the reduced PIF with random threshold. The functional description of the reduced PIF is obtained by setting the phase deviation in (17) to zero.

Definition 3.2. The reduced PIF neuron with random threshold is given by the equations
where (ε i k ), k ∈ Z, is a sequence of independent Gaussian random variables with zero mean and variance For reasons of notational simplicity and without loss of generality, and unless otherwise stated, we shall assume here that B = I (N i = P i ). The reduced PIF (rPIF) with random threshold can now be written as where (ε i k ), k ∈ Z, i = 1, 2, is a sequence of independent Gaussian random variables with zero mean and variance The above analysis assumes that the inputs are weak and therefore the BSGs operate on a limit cycle. Stronger signals can be taken into account by considering a manifold of PRCs associated with a wide range of limit cycles (Kim and Lazar, 2012). By estimating the limit cycle and hence its PRC using spike times, we have the following I/O relationship for each of the BSGs.
Definition 3.3. The reduced PIF neuron with conditional PRC and random threshold is given by the system of equations is an initial bias that brings the neuron close to the spiking region in the absence of input and (by abuse of notation) ε i k , k ∈ Z, i = 1, 2, is a sequence of independent Gaussian random variables with zero mean and variance and ψ i 1 (s, b i k ) is the conditional PRC (Kim and Lazar, 2012).
The conditional PRC formulation above allows us to separate BSG inputs into a constant bias current and fluctuations around it on short inter-spike time intervals. The bias current can be estimated between consecutive spikes, making the deviation from the limit cycle small in each inter-spike interval even for strong inputs. Moreover, by considering the conditional PRCs, the assumption that BSGs oscillate in the absence of input can be nearly dropped. Thus, it is not required for BSGs to always be on a limit cycle. Only when the neuron enters the limit cycle do we consider formulating the encoding using the rPIF model with conditional PRCs.
Remark 3.4. Note that by parametrizing each of the PRCs with b i k , the variance of the error in (24) depends on the estimated PRC on each inter-spike interval. In conjunction with (23), we see that the variability of spike times depends on the strength of the input to the BSGs.

The t-transform of the neural circuit
The overall encoding by the neural circuit model can be expressed as Substituting (12) into the above, we have We arrived at the following.
Lemma 3.5. The model of encoding in Figure 1 is given in operator form by and i k , k ∈ Z, are independent random variables with normal distribution N (0, 1) and j = 1, 2, j = i. Equation (26) is called the t-transform (Lazar and Tóth, 2004) of the neural circuit in Figure 1. Remark 3.6. The t-transform describes the mapping of the input stimulus u 1 into the spike timing sequence (t i k ), i = 1, 2, k ∈ Z. Thus, the t-transform shows how the amplitude information of the input signal is related to or transformed into the time information contained in the sequence of output spikes generated by the neural circuit.
We provide here further intuition behind the Equations (26). By the Riesz representation theorem (Berlinet and Thomas-Agnan, 2004), there exists functions φ i 1k ∈ H 1 1 such that , for all u 1 ∈ H 1 1 , Therefore, (26) can be rewritten in inner product form: Recall that inner products are projections that are typically interpreted as measurements. In the Equation (27) above, the signals u 1 and u 2 are projected onto the sampling functions φ i 1k and φ i 2k , respectively. We also note that traditional amplitude sampling of a bandlimited signal u 1 at times (t n ), n ∈ Z, can be expressed as where sinc (t) = sin ( 1 t) πt is the impulse response of the ideal low pass filter bandlimited to 1 or in other words, the kernel of the RKHS of finite-energy band-limited functions (Lazar and Pnevmatikakis, 2009). Thus, the neural encoding model described by the Equation (27) can be interpreted as generalized sampling with noisy measurements with sampling functions φ i 1k and φ i 2k . The formulation of the encoding model can easily be extended to the case when M neural circuits encode a stimulus in parallel. This is shown schematically in Supplementary Figure 1. A left superscript was added in the figure to each of the components to indicate the circuit number.

DECODING
In the previous section, we showed that the encoding of a signal u 1 by the neural circuit model with feedforward and feedback DSPs and BSGs can be characterized by the set of t-transform Equations (26). We noticed that the Equations (26) are nonlinear in u 1 due to the second order Volterra term. However, by reinterpreting the second order term as linear functionals T i 2k on the higher dimensional tensor space H 1 2 , (26) implies that the measurements taken by each of the neurons are the sum of linear measurements in two different vector spaces [see also Equations (27)].
In this section we investigate the decoding of signals encoded with the neural circuit in Figure 1. The purpose of decoding is to recover from the set of spike times the original signals, u 1 (t) and u 2 (t 1 , t 2 ), that respectively belong to the two different vector spaces H 1 1 and H 1 2 . We formulate the decoding problem as the joint smoothing spline problem û 1 ,û 2 = argmin where n i + 1 is the number of spikes generated by BSG i = 1, 2. Theorem 3.7. The solution to (28) is of the form T is the solution of the system of linear equations where q = q 1 1 , · · · , q 1 n 1 , q 2 1 , · · · , q 2 n 2 T , and Proof: Proof of the theorem follows the Representer Theorem (Berlinet and Thomas-Agnan, 2004) and is given in detail in Appendix.
Remark 3.8. When λ 1 = λ 2 , the solution c amounts to where I is an identity matrix of appropriate dimensions.
Remark 3.9. Although (29) solves (28), in practice a minimum number of spikes is needed to obtain a meaningful estimate of the original signal. A minimum bound for the number of measurements/spikes can be derived in the noiseless case. Clearly, the bound has to be larger than the dimension of the space. This may require the signal to be encoded by a circuit with a larger number of neurons than the two shown in Figure 1 (Lazar and Slutskiy, in press).
A number of such neural circuits in parallel can be used to encode input stimuli as shown in the Supplementary Figure 1. Theorem 3.7 can be easily extended to solving the smoothing spline problem û 1 ,û 2 = argmin where m = 1, 2, . . . , M, denotes the circuits number in Supplementary Figure 1. In addition, if the circuits consist of only first order feedforward kernels, then only u 1 (t) can be reconstructed. Similarly, if the circuits are comprised of only the second order feedforward kernels, then u 2 (t 1 , t 2 ) can be reconstructed but not u 1 (t).
Remark 3.11. The reconstruction of u 2 (t 1 , t 2 ) on D 2 strongly depends on the support (in practice the finite memory) of the kernels h 11i 2 , i = 1, 2 (see also Figure 2C). In the reconstruction example of the Supplementary Figure 5, we show thatû 2 approximates u 2 well in the restricted domain where h 11i 2 is nonzero. Outside this restricted domain, h 11i 2 vanishes and u 2 is not well recovered as suggested by the large error in the Supplementary Figure 5E.

EFFECT OF NOISE ON STIMULUS DECODING
In this section, we investigate the effect of noise sources (i) on spike timing of the reduced PIF neuron, and (ii) on the decoding of stimuli encoded with a neural circuit. We will also present the effect of an alternative noise source model on both spike timing and stimulus decoding.

Effect of noise on measurement and spike timing errors of the reduced PIF neuron
As suggested by (22), the variance of the measurement error of the reduced PIF neuron is directly related to the PRC of the associated limit cycle. We first characterize the variance of the measurement error due to each individual noise source parametrized by the bias current I i . We then evaluate the spike timing variance between the spike trains generated by the Hodgkin-Huxley neuron and the reduced PIF neuron again as a function of the bias current I i . We start with a brief description of the key elements of Hodgkin-Huxley neuron and the PIF neuron. We consider the stochastic Hodgkin-Huxley equations where f i is defined as in Section 2.2 of the Supplementary Material with additional normalization such that the unit of time is in seconds instead of milliseconds and the unit of voltage is in Volts instead of milivolts as conventionally used. Z i (t) takes the form Here W i n (t) are independent standard Brownian motion processes and σ i n , n = 1, 2, 3, 4, are associated scaling factors. The variance of the measurement error of the reduced PIF neuron due to each Brownian motion process W i n , n = 1, · · · , 4, is given by [see also Equation (22)] We show in Figure 3A the variance of the measurement error in (32) associated with each source of noise of the reduced PIF neuron for the unitary noise levels σ i n = 1, n = 1, 2, 3, 4. The variances given by (32) are plotted as a function of the bias current I i . Clearly, the noise arising in dendritic stimulus processing (W i 1 ) induces the largest error, and together with noise in the potassium channels (W i 2 ), these errors are about two magnitudes larger in variance than those induced by the noise sources in the sodium channels (W i 3 , W i 4 ). The above analysis is based on the analytical derivation of the measurement error in (32) for the rPIF neurons. The measurement error is closely related, however, to the spike timing variation of the BSGs subject to noise sources. A variance of 10 −6 in Figure 3A corresponds to a standard deviation of 1 ms in spike timing. In practice the error between the spike times of the Hodgkin-Huxley neuron and the reduced PIF neuron can be directly evaluated.
In order to do so, we randomly generated a weak bandlimited dendritic input. All evaluations were based on encoding a signal with the Hodgkin-Huxley neuron model described above with internal noise sources and bias current I i . The spike times (t i k ) of the Hodgkin-Huxley neuron were recorded. Starting from each spike time t i k , we encoded the appropriate portion of the signal by the reduced PIF neuron until a spike r t i k + 1 was generated. The difference between r t i k + 1 and t i k + 1 is the error in approximating the encoding using the reduced PIF formulation. This process was repeated for each I i . We computed the variance of the errors based on some 3000-5000 spikes generated in encoding the input.
In Figure 3B, the variance of the spike timing error r t i k + 1 − t i k + 1 for σ n = 0, n = 1, 2, 3, 4, is shown. Since the reduced PIF is an approximation (even under noiseless conditions) and, although small, the error is nonzero. From Figure 3B, the variance of the spike timing error is on the order of 10 −9 . We shall evaluate the spike timing error variance of the intrinsic noise sources in a range much larger than 10 −9 .
We also tested to what extent each individual source of noise contributes to the variance of spike timing as suggested by the theoretical analysis depicted in Figure 3A. Indeed, the error variance obtained through simulations in Figure 3C follows the basic pattern shown in Figure 3A. Figure 3C was obtained by setting one of the σ n 's to a nonzero value and the rest to 0 (the nonzero values were σ 1 = σ 2 = 0.01, σ 3 = σ 4 = 0.1). Each nonzero value was picked to be large enough so that the error variance in the absence of noise ( Figure 3B) becomes negligible, and at the same time, it was small enough such that the states of the neurons did not substantially deviate from the limit cycle. To compare the with the ones in Figure 3A we normalized the error variance obtained in simulations by σ n . Next, we tested whether the variance of spike timing due to presence of multiple noise sources is truly the summation of error variances due to individual noise sources. We simulated the Hodgkin-Huxley equations with σ 1 = σ 2 = 0.005, σ 3 = σ 4 = 0.05. The total spike timing error variance shown in Figure 3D (blue curve) is very close to the sum of error variances in Figure 3C with proper scaling (red curve in Figure 3D).
As suggested by the above analysis, the reduced PIF neuron with random thresholds largely captures the encoding of stimuli by BSGs subject to intrinsic noise sources.

Effect of noise on stimulus decoding
In order to quantitatively explore how noise impacts signal decoding, we recovered from spikes the signal encoded by the noisy neural circuit of Supplementary Figure 1. We started with the base-level noise-less case described in Section 3.2 of the Supplementary Material (M = 4) and proceeded to introduce individual noise terms with a range of scaling factors. For example, we set σ i 2 = σ i 3 = σ i 4 = 0 and varied σ i 1 . We also tested the case when 10σ i 1 = 10σ i 2 = σ i 3 = σ i 4 for the aggregated effect on stimulus recovery. We choose to use σ i 3 and σ i 4 10 times larger than σ i 1 and σ i 2 so that each noise source introduced a similar error.
In all simulations, the Euler-Maruyama scheme (Kloeden and Platen, 1992) was used for the numerical integration of the SDEs. We performed 20 encoding and decoding experiments. Each time, the input stimulus was generated by randomly picking from a Gaussian distribution the real and imaginary parts of the coefficients u l in (1). We further constrained the stimuli to be real-valued. (An example is given in Supplementary Figure 5.) For each noise level, the input signal was encoded/decoded. The mean Signal-to-Noise Ratio (SNR) across 20 experiments is reported for each noise level. The SNR for the reconstruction of u 1 was computed as SNR = 10 log 10 u 1 where u 1 is the original signal andû 1 is its reconstruction. Note that the spike time occurrences generated for the same signal are different for each noise level. Since the sampling functions are spike time dependent, the number of measurements/spikes may not be the same for each noise level. Moreover, at times, the sampling functions may not fully span the stimulus space. To reduce the uncertainty caused by the stimulus dependent sampling we averaged our SNR data over 20 different signals. Figure 4A shows the SNR of the reconstruction of signal u 1 (t) against different noise strength. Figure 4B shows the SNR of the reconstruction of signal u 2 1 (t) = u 2 (t, t). The reconstruction SNR in Figure 4A largely matches the inverse ordering of noise strength of each of the individual noise sources shown in Figure 3A. DSP noise sources degrade the reconstruction performance most strongly while noise sources associated with gating variables m and h have a much smaller effect for the same variance level. Since the variance of measurement error is the sum of error variance in each variable, the case when 10σ 1 = 10σ 2 = σ 3 = σ 4 = σ exhibits the lowest performance.

Effect of an alternative noise model on spike timing and stimulus decoding
Biologically, the effect of channel noise on the operation of the BSGs is due to the ON-OFF activity of a finite number of ion channels. The Hodgkin-Huxley equations and the noise terms used in Section 3.3.2 correctly capture the dynamics in the limit of infinitely many channels. Recent research, however, suggests that the model equations may not correctly model the ion current fluctuations for a finite number of channels . We consider here an alternative stochastic formulation of the Hodgkin-Huxley model that more precisely captures the ion channel kinetics. By using a finite number of ion channels the strength of noise amplitude becomes directly related to the actual number of ion channels. Therefore, the free variables are only the number of potassium and sodium channels that are both biologically meaningful. The successful use of an alternative noise model as described below also suggests that our analysis can be applied to a wide range of stochastic formulations of BSGs based on SDEs.
We shall construct here stochastic ion channels using conductance noise rather than subunit noise as in the previous Sections . This stochastic Hodgkin-Huxley system is simulated using a diffusion approximation following (Orio and Soudry, 2012). The system of SDEs can be expressed as where Y i has 14 state variables and the full system can be found in Section 3.3 of the Supplementary Material. Here i = 1 for simplicity.
The variance of the measurement error is now given by (20). We can decompose the variance into three terms as where ε i kV , ε i kK , ε i kNa are measurement errors associated with the noise in the DSP, in potassium channels and in sodium channels, respectively.
As ε i kV is quantitatively the same as that in Section 3.3.2, we focus here on ε i kK and ε i kNa . The variance of the errors can be respectively expressed as Note that b np , n = 1, · · · , 14, p = 2, 3, · · · , 15, are functions that dependent on either the number of potassium channels N Na or the number of sodium channels N K , and the states of the neuron.
We first evaluate (E ε i kNa 2 )(I i ) using the PRCs. The PRCs are obtained by letting N Na = N K = ∞ and thereby making the system deterministic. Since the measurement error variance for fixed I i is proportional to (N Na ) −1 , it is shown in Figure 5A as a function of the bias current I i for N Na = 1. Similarly, the variance of the measurement error E ε i kK 2 (I i ) for N K = 1 is shown in Figure 5A, and it is proportional to (N K ) −1 for a fixed I i . We notice that, when the number of channels is the same, the measurement error due to the sodium channels is of the same order of magnitude with the measurement error due to the potassium channels. However, the number of sodium channels is Frontiers in Computational Neuroscience www.frontiersin.org September 2014 | Volume 8 | Article 95 | 13

FIGURE 5 | The variance of the measurement and spike timing error associated with the sodium channels (blue) and the potassium channels (red) of the Hodgkin-Huxley equations with alternative noise sources parametrized by the bias current I. (A)
The variance of the measurement error computed from PRCs of Hodgkin-Huxley equations, with N Na = 1 and N K = 1. Actual variance with different number of ion channels is inversely proportional to N Na and N K , respectively. (B) Spike timing variance obtained in simulations by comparing the spike times generated by the Hodgkin-Huxley neuron with channel noise and the spike times generated by its reduced PIF counterpart. Blue curve is obtained by using N Na = 5 × 10 4 , N K = ∞, and normalized to 1 sodium channel. Red curve is obtained by using N K = 5 × 10 4 , N Na = ∞, and normalized to 1 potassium channel.
typically 3-4 times larger than the number of potassium channels. Therefore, in contrast to the previous case, the error induced by sodium channels shall be larger than that induced by potassium channels. We also analyzed in simulations the difference between spike times generated by the alternative stochastic formulation of the Hodgkin-Huxley equations and those generated by the corresponding reduced PIF neuron. We used in simulation N Na = 5 × 10 4 , N K = ∞, to obtain the variance E ε i kNa 2 (I i ) and scaled it by N Na to compare it with Figure 5A. Similarly, we used N K = 5 × 10 4 , N Na = ∞, to obtain the variance E ε i kK 2 (I i ).
The spike timing variances of error across different I i are shown in Figure 5B The pattern of similarity between variances in Figures 5A,B suggest that the reduced PIF with random threshold associated with this formulation of Hodgkin-Huxley equations is highly effective in capturing the encoding under internal noise sources.
We now show how ion channel noise sources affect the decoding of the input signal. We varied the number of sodium channels N Na and fixed the number of potassium channels to be N K = 0.3N Na , a ratio typically used for Hodgkin-Huxley neurons with the alternative noise source model. By decoding the input stimulus we show how increasing the number of ion channels improves the faithfulness of signal representation. The SNR of the reconstruction of u 1 (t) and u 2 1 (t) are depicted in Figure 6. SNR goes down to about 4 dB when 600 sodium channels and 180 potassium channels are used. This corresponds to a membrane area of about 10 µm 2 with a density of 60 µm 2 in sodium channels and 18 µm 2 in potassium channels . We also tested the reconstruction for the case when one type of ion channels is infinitely large, i.e., deterministic, while varying the number of ion channels of the other type. The result is also shown in Figure 6. The noise from the dendritic tree shall have similar effect on the representation since the voltage equation is the same as in Section 3.3.2.

FUNCTIONAL IDENTIFICATION AND NOISE
In Section 4.1 we show how to functionally identify the feedforward and feedback DSPs of the circuit in Figure 1 under noisy conditions. We assume here that the BSGs have already been identified using a methodology such as the one developed in Lazar and Slutskiy (2014). In Section 4.2 we discuss the effect of noise parameters on the quality of DSP identification.

FUNCTIONAL IDENTIFICATION
In the decoding problem analyzed in Section 3.2, we reconstructed unknown input stimuli by assuming that the neural circuit in Figure 1 is known and the spike trains are observable. In contrast, the objective of the functional identification problem investigated in this section is to estimate the unknown circuit parameters of the feedforward and feedback DSPs from I/O data.
The I/O data is obtained by stimulating the circuit with controllable inputs and by measuring the time occurrences of the output spikes. This basic methodology has been a standard practice in neurophysiology for inferring the function of sensory systems (Hubel and Wiesel, 1962). We assume here that either the BSGs are known in functional form or the family of PRCs associated with the BSGs have already been identified (Lazar and Slutskiy, 2014).
Although decoding and functional identification are seemingly two different problems, they are closely related. By exploiting the commutative property of linear operators, we can rearrange the diagram of the neural circuit model of Figure 1 into the form shown in Figure 7. We notice that the outputs of Figure 7 and those of Figure 1 are spike time equivalent, as long as the RKs K 2 1 and K 2 2 have large enough bandwidth. In what follows we will evaluate the four Volterra terms, i.e., the four dendritic currents feeding the BSG of Neuron 1 in Figure 7.
Formally, the first order (feedforward) Volterra term can be written as (Lazar and Slutskiy, in press Similarly, the second order (feedforward) Volterra term amounts to The above equations suggest that the projections of the feedforward kernels can be re-interpreted as inputs, whereas the signals u 1 and u 2 can be treated as feedforward kernels. In Section 2.2.2 we introduced two RKHSs, H 2 1 and H 2 2 , for modeling two different spaces of spikes. The elements of H 2 1 are functions defined over the domain [0, The period S 2 is large enough to ensure that any spike that arrives supp{h 2ji 1 } seconds prior to the arrival of t i k , or earlier, will not affect the output of the feedback kernel on the inter-spike time interval [t i k , t i k + 1 ]. Thus, such spikes will not introduce additional error terms in the t-transform evaluated on the inter-spike time Note that the domain [0, S 2 ] of the functions in H 2 1 may not be the same as the domain of the input space H 1 1 . However, such a domain can be shifted on a spike by spike basis to . This is important for mitigating the practical limitation of modeling the stimuli as periodic functions in H 1 1 . The response of the first-order feedback term to its spiking input on the inter-spike time interval [t i k , t i k + 1 ] in Figure 7 amounts to (i = j) FIGURE 7 | Diagram of the neural circuit that is spike timing equivalent with the one in Figure 1 highlighting the duality between neural decoding and functional identification. Note that the input stimuli and the DSP projections are reordered to reflect that the unknowns are the DSP projections. The input stimuli u 1 1 (t), u 1 2 (t 1 , t 2 ), and the kernel representation of spikes (see also Section 2.2.2) are intrinsic to the neural circuit. The DSP projections are interpreted as inputs.

Frontiers in Computational Neuroscience
www.frontiersin.org September 2014 | Volume 8 | Article 95 | 15 It is clear from Section 2.2.2 that if 2 is at least larger than the effective bandwidth of h 2ji 1 and L 2 → ∞.
Similarly, the response of the second-order feedback kernel to its spiking input on the inter-spike time interval [t i k , t i k + 1 ] amounts to if 2 is large enough and L 2 → ∞. Combining (34), (36), (36), and (37), for each spike interval [t i k , t i k + 1 ], the aggregated output current of the DSPs of Neuron i in Figure 7, shall converge to the DSP aggregated output current of Neuron i in Figure 1 for large enough 2 . A direct consequence of this equivalence is that, under the same additive Gaussian white noise and channel noise in the BSGs, the t-transform of the circuit in Figure 7 and in Figure 1 are identical.
Note that the outputs of the feedforward kernels are always equivalent; the equivalence of the outputs of the feedback kernels requires, however, the use of large enough bandwidth 2 . Otherwise, the equivalence in the t-transform is invalid and an additional noise term appears in the t-transform of the Neuron 1 in Figure 7.
The projections of the Volterra DSP kernels of Figure 7 are interpreted as inputs, while the input stimuli and the train of RKs at spike times replace the impulse response of the corresponding filters. Therefore, the functional identification problem has been transformed into a dual decoding problem, where the objects to decode are the set of projections of Volterra DSP kernels and the neural circuit is comprised of "stimulus DSP kernels" and "spike DSP kernels" with the same BSGs and noise sources. The only difference is that, instead of a Single-Input Multi-Output decoding problem, the identification was transformed into a Multi-Input Multi-Output decoding problem. In addition, multiple trials using different stimuli are needed; this procedure is illustrated in block diagram form in Figure 8. By stimulating the neural circuit with multiple stimuli in the functional identification setting, multiple neural circuits effectively encode the projections of the DSP kernels.
We are now in the position to derive the t-transform of Neuron 1 in Figure 7. Assuming that m = 1, · · · , M, trials are performed for identification, the t-transform (26) can be written as for i, j = 1, 2, i = j, k ∈ Z. Here m L 1i 1k : H 1 1 → R, m L 1i 2k : H 1 2 → R are bounded linear functionals associated with the feedforward DSP kernels, and m L 2i 1k : H 2 1 → R, m L 2i 2k : H 2 2 → R are bounded linear functionals associated with the feedback DSP kernels for each trial m. The above functionals are defined as and m i k , i = 1, 2, k ∈ Z, m = 1, · · · , M, are independent random variables with normal distribution N (0, 1).
The functional identification of the neural circuit in Figure 7 can then be similarly defined to the decoding problem. We formulate the identification of the noisy neural circuit again as two smoothing spline problems, one for each neuron, where m n i is the number of spikes generated by Neuron i in trial m.
The solution can be obtained in a similar way as in Theorem 3.7.
Since each of the kernel projections may be in a different RKHS, and their domain may also be different, the identification of all filters resemble that of the multi-sensory Time Encoding Machines. Recall that multi-sensory TEMs encode within the same circuit time-varying and space-time varying sensory signals while decoding remains tractable (Lazar and Slutskiy, 2013). The solution to (41) can similarly be obtained as the solution to (40) above.
Note that we are only able to identify the projection of the Volterra kernels. This is because, in practice, we can only probe the system with signals in a bandlimited space. By increasing the bandwidth of the elements of the Hilbert space, the projection of the kernels will converge to their original form (Lazar and Slutskiy, 2012).

Remark 4.2.
It is important to note that in order to have a good estimate of the kernels, stimuli must fully explore all input spaces. This can be quite easily achieved for the feedforward DSP kernels by using many (randomly generated) signals that cover the entire frequency spectrum. However, to properly identify the feedback DSP kernels, spike trains must be diverse enough to sample its different frequency components. This may not be easy to realize in practice. For first order feedback kernels, spike trains with constant spike intervals are, for example, undesirable. The Fourier transform of regular Diracdelta pulses is a train of Dirac-delta pulses in the Fourier domain. This means that only certain frequency responses of the DSP kernels are, for example the DC component, sampled. The rest of the frequency components are essentially in the null space of the sampling functions m φ ik , i = 1, 2, m = 1, · · · , M. Similar effect applies to the space of trigonometric polynomials. If the spike intervals exhibit small variations, many of the frequency components may be sampled but the energy at DC may be too dominant. In this case, noise may contaminate more severely the measurement of non-DC components and may yield unsatisfactory identification. This effect may pose even more stringent requirements on the identification of the second order feedback kernels, as it requires the interaction between two spike trains.

EFFECT OF NOISE ON IDENTIFICATION
In order to evaluate the effect of noise on the identification of the neural circuit in Figure 1 we included intrinsic noise into the example neural circuit discussed under noiseless conditions in Section 4.1 of the Supplementary Material. Randomly generated signals were used in the identification examples given here. Chosen in the same way as in the decoding example in Section 3.3.2 all these signals are used here to identify the neuron in question. Therefore, in this section, multiple signals are used in repeat experiments to identify the parameters of a neural circuit. By contrast in Section 3.3.2, multiple neurons are used to encode a single signal.
First, we evaluated the effect of noise on the quality of identification of DSP kernels of Neuron 1 in Figure 7 with a BSG modeled by the SDE (31) with 10σ i 1 = 10σ i 2 = σ i 3 = σ i 4 = σ . Figure 9 shows the SNR of the identified DSP kernels in Figure 7 across several noise levels σ . As expected, the general trend for all four kernels is decreasing SNR with increasing noise levels. We notice that the identified feedforward DSP kernels have similar shape as the original kernel, even at high noise levels. However, the feedback DSP kernels undergo a change in shape at high noise levels. We can see that the time instance of the peak amplitude in the first order feedback kernel is shifted to an earlier time instance.
Second, we investigated the identification of DSPs for a BSG noise model already described in Section 3.3.3. Figure 10 shows the SNR of the identified DSP kernels across a different number of sodium channels N Na while N K = 0.3N Na . The SNR plots suggest that the identification quality increases as more ion channels are present in the BSGs.
Additionally, as discussed in Remark 4.2, BSG noise sources may degrade severely the identification of feedback kernels when the spike trains generated in trials are not sufficient for exploring the two spike input spaces. We show an example of the later in Figure 11. The two BSGs have higher bias currents and lower input current magnitude. The later was achieved by scaling down the magnitude of the DSP kernels. The combined effect results in regular spiking intervals in both neurons. The identification result under noiseless conditions is shown in Figure 11. Note that since the t-transform of the Hodgkin-Huxley neuron is not exact, a small error is introduced even if intrinsic noise is not present. We see that the feedforward DSP kernels can be identified quite well, yielding SNRs of around 17 dB. However, the feedback DSP kernels are not well identified. In particular, the identified secondorder feedback kernel has a wide spread, similar to the high noise case in Figure 9D. This suggest that the spike pattern is not sufficiently exploring the space of feedback kernels. A large number of frequency components are only weakly sampled and they can be very easily contaminated by noise. The presence of both intrinsic noise sources can exacerbate the condition further. This effect is confirmed with a simulation of the circuit using Integrate-and-Fire (IAF) neurons. Since the t-transform for the IAF neuron is exact (Lazar and Tóth, 2004), both feedback kernels can be identified even if the generated spikes only weakly explore certain frequency components. However, by artificially adding a small measurement error to the t-transform of the circuit with IAF neurons, similar results to those in Figure 11 can be obtained (data not shown).

DISCUSSION
In this paper, we introduced a novel neural circuit architecture based on a neuron model with a biophysical mechanism of spike generation and feedforward as well as feedback dendritic stimulus processors with intrinsic noise sources. Under this architectural framework, we quantitatively studied the effect of intrinsic noise on dendritic stimulus processing and on spike generation. We investigated how intrinsic noise sources affect the stimulus representation by decoding encoded stimuli from spikes, and quantified the effect of noise on the functional identification of neural circuits. We argued that a duality between stimulus decoding and functional identification holds. Therefore, the encoding framework based on the neural circuit architecture studied here can be applied to both the reconstruction of the encoded signal and the identification of the dendritic stimulus processors. We systematically showed how the precision in decoding is affected by different levels of stochastic variability within the circuit. These results apply verbatim to the functional identification of dendritic stimulus processors via the key duality property mentioned above.
Our theoretical framework highlights two indispensable components of modeling signal representation/processing in a neural circuit-dendritic stimulus processing and spike generation. Such a divide and conquer strategy is ubiquitous in engineering circuits and leads to a separation of concerns. Recent experimental studies also showed that interesting nonlinear processing effects take place first in the dendritic trees rather than in the axon hillock (Yonehara et al., 2013).
We presented here two types of nonlinear dendritic stimulus processors. The first type are feedforward DSPs that respond to continuous input sensory stimuli. The second type are feedback DSPs that respond to spiking inputs. We quantitatively demonstrated how intrinsic noise sources would affect the identification quality of all these DSPs. The examples in Section 4.2 suggest that in identification feedback kernels are more vulnerable to internal noise sources than feedforward kernels. In particular, the overall shape of the identified feedback kernels differs significantly from that of the underlying kernels when the strength of noise sources becomes large. Meanwhile the identified feedforward kernels are qualitatively preserved, albeit not accurately.
Most of the single neuron models (LIF, QIF) in the literature have focused on the spike generation mechanism. The encoding capability of these models is typically investigated based on rate encoding (Eliasmith and Anderson, 2003;Lundstrom et al., 2008;Ostojic and Brunel, 2011). For both decoding and identification we used here the occurrence times of spikes generated by spiking neuron models. Most importantly, the BSG models discussed here were characterized by a PRC manifold (Kim and Lazar, 2012) in the presence of noise, while many simplified models (such as the LIF) can be effectively described with a single PRC. Other sensory neuron models, e.g., GLM (Pillow et al., 2011), usually rely on a rate-based output or Poisson spike generation that do not take into account key advances in dynamical systems-based spiking neuron models.
As already mentioned before, we investigated how intrinsic noise sources affect the stimulus representation by decoding encoded stimuli from spikes. We are not suggesting, however, that the decoding algorithm considered here is implemented in the brain. Rather, we argue that decoding is effective in measuring how well information is preserved in the spike domain. The decoding formalism allowed us to investigate how noise affects the fidelity of signal representation by a population of neurons by reconstructing stimuli and comparing their information content in the stimulus space.
While decoding can serve as an "oscilloscope" in understanding stimulus representation in sensory systems, functional identification serves as a guide in experiments to functionally identify sensory processing. Based on spike times, the identification algorithm presents a clear bound on the number of spikes that are necessary for perfect identification under noiseless conditions. Phrased differently, when a certain number of spikes are acquired from a neuron of interest, the identification algorithm places a constraint on the maximum DSP kernel bandwidth that can perfectly be recovered. In more practical terms, we advanced two important applications of the circuit architecture considered in this paper. The first one considers dendritic stimulus processors that process information akin to complex cells in V1. The second one adapts the widely used Hodgkin-Huxley model known in the context of neural excitability (Izhikevich, 2007) and analysis of neuronal stochastic variability to stimulus encoding in the presence of noise.
Based on the rigorous formalism of TEMs (Lazar and Tóth, 2004), we extended our previous theoretical framework  and argued that spike timing is merely a form of generalized sampling of stimuli. By studying sampling (or measurements) in the presence of intrinsic noise sources, we showed to what extent neurons can represent sensory stimuli in noisy environments as well as how much noise the identification process can tolerate while preserving an accurate understanding of circuit dynamics.
The reconstruction and identification quality are certainly not only related to the strength of noise, but also the strength of the signal. In particular, when the signal strength is small, two facts may affect the quality of reconstruction. First, neurons may not produce enough spikes that have valid t-transforms. Second, they may be contaminated by even weak noise, i.e., the signal-to-noise ratio is low. It is well known, however, that neural systems use gain control to boost the relevant signal (Shapley and Victor, 1978;Wark et al., 2007;Friederich et al., 2013). Such strategy may be useful for increasing the signal strength relatively to the strength of the noise. Gain control may also suppress large signals to fit into the range of operation of the BSGs. The gain control itself, maybe considered as a type of Volterra feedforward DSP kernel (Lazar and Slutskiy, in press) and the interaction with feedback loops driven by spikes. The lack of spikes may be compensated by adding other neurons that are sensitive to other features in the input stimuli.
A key feature in our neural circuit model is the nonlinear processing in the feedforward and feedback paths. Nonlinear interaction between feedforward DSPs and feedback DSPs have not been considered here. However, they are of interest and could be addressed in the future. Self-feedback was not included in the model for clarity, but can be considered within the framework of our approach. Self-feedback was introduced to add refractoriness to phenomenological neuron models (Keat et al., 2001;Pillow et al., 2008). Our BSG models, on the contrary, are conductance-based models that have a refractory period built in.
Throughout this paper we assumed that the BSGs themselves have been perfectly identified. The intrinsic noise in the BSGs may degrade the identification quality of conditional PRCs. This may result in a lower identification quality as shown in the examples. It is beneficial to investigate in the future a method that can identify the entire circuit at once so that intrinsic noise in the circuit only affects the identification process a single time.
The theoretical results obtained here suggest a number of experiments in the early olfactory system of fruit flies. The glomeruli of the antennal lobe can be modeled using the Volterra DSPs discussed here and the projection neurons in the antennal lobe are accessible by patch clamping (Lazar and Yeh, 2014). Functional identification of DSPs can then be carried out for studying olfactory stimulus processing in an accessible circuit with intrinsic noise sources (Masse et al., 2009).