Channel identification machines for multidimensional receptive fields

We present algorithms for identifying multidimensional receptive fields directly from spike trains produced by biophysically-grounded neuron models. We demonstrate that only the projection of a receptive field onto the input stimulus space may be perfectly identified and derive conditions under which this identification is possible. We also provide detailed examples of identification of neural circuits incorporating spatiotemporal and spectrotemporal receptive fields.

1. Introduction

Recently, a novel class of system identification algorithms, called Channel Identification Machines (CIMs), has been developed for identifying dendritic processing in spiking neural circuits (Lazar and Slutskiy, 2010, 2012). In the simplest setting, CIMs allow one to identify a communication/processing channel in [Filter]-[Asynchronous Sampler] circuits, where the effect of the channel on an input signal can be described by a linear filter and the output of the channel is mapped, or encoded, into a time sequence by a non-linear asynchronous sampler. Such [Filter]-[Asynchronous Sampler] circuits are known in the literature as Time Encoding Machines (TEMs) and their specific embodiments in neuroscience include neural circuit models in which an analog dendritic stimulus processor is followed by a point neuron model encoding the aggregate dendritic current produced by the processor. A few examples of asynchronous samplers include asynchronous sigma/delta modulators (ASDMs), conductance-based point neuron models such as Hodgkin-Huxley, Morris-Lecar, Fitzhugh-Nagumo, Wang-Buzsáki, Hindmarsh-Rose, oscillators with multiplicative coupling and simpler models such as the integrate-and-fire neuron or the threshold-and-fire neuron (Brown et al., 2004; Lazar and Slutskiy, 2012, 2014a).

The above-mentioned asynchronous samplers incorporate the temporal dynamics of spike generation at the axon hillock of biological neurons and allow one to consider, in particular for neuroscience applications, more biologically-plausible non-linear spike generation (sampling) mechanisms. This is in contrast to existing methods, such as the generalized linear model (GLM), which typically assumes a simplified description of the spike generation dynamics by using a static non-linearity. The non-linear contribution of a dynamical system such as the Hodgkin-Huxley neuron model is not static, but dynamic and stimulus-driven (Lazar and Slutskiy, 2014a). It changes from one spike to the next and thus affects the estimation procedure if not properly taken into account. Furthermore, since the non-linearity does not fully capture the complex spike generation dynamics, the filters fit to that non-linearity may not provide an adequate description of the neural circuit, and in particular may confound dendritic processing and encoding (Lazar and Slutskiy, 2014a).

More recently, the CIM methodology has been extended to neural circuits in higher brain centers that receive multidimensional spike trains as input stimuli instead of continuous signals and to circuits with extensive lateral connections and feedback (Lazar and Slutskiy, 2014a). Together with TEMs and Time Decoding Machines (TDMs) for decoding stimuli from spike trains, CIMs provide a unified framework for studying encoding, decoding and identification in neural systems.

The motivation for multidimensional CIMs (Lazar and Slutskiy, 2013) is provided by the concept of a receptive field that is well established in neuroscience. Introduced in 1906 by Sherrington (Sherrington, 1906) to describe an area of the body surface capable of eliciting a reflex in response to a stimulus, the term “receptive field” has been extended to many different sensory modalities and spans many different types of neurons. For example, in the visual system, the receptive field of a photoreceptor is a 3-dimensional cone in space comprising all possible directions in which light can hit the photoreceptor. In the auditory system, receptive fields can correspond to certain spectral regions of audio stimuli. More broadly, the receptive field is that part of the the sensory space that can evoke a neuronal response (Dayan and Abbott, 2005). Spatial and spatiotemporal receptive fields have been successfully used in vision to model retinal ganglion cells in the retina as well as neurons in the lateral geniculate nucleus and the visual cortex (see DeAngelis et al., 1995 for a review). Similarly, spectrotemporal receptive fields have been used to describe responses of auditory neurons (Aertsen and Johannesma, 1981), neurons in cochlear nuclei (Clopton and Backoff, 1991) and neurons in the auditory cortex (Kowalski et al., 1996).

Circuits that process multidimensional feedforward input are often encountered in biological systems (e.g., the retina, Pillow et al., 2008). From a modeling standpoint, one can also consider spiking neural circuits in which, in addition to feedforward input, every neuron may receive lateral inputs from neurons in the same layer. Furthermore, back-propagating action potentials (Waters et al., 2005), or feedback, may affect computations performed within the dendritic tree. Until now, CIMs for circuits with lateral connectivity and feedback have only been considered in the context of scalar or vector-valued signals of one variable, e.g., functions of time u₁(t), t ∈ ℝ.

In this paper we discuss multidimensional channel identification machines that allow one to identify signal transformations applied to multidimensional signals u_n(x₁, …, x_n), n ∈ ℕ, where x_n typically designates the time variable. A few examples of multidimensional CIMs include (i) spatial CIMs, where the input signal u₂(x, y) is a function of a two-dimensional space, describing, e.g., images; (ii) spectrotemporal CIMs, where the input signal u₂(ν, t) is a function of spectrum and time, describing, e.g., auditory signals; (iii) spatiotemporal CIMs, where the input signal u₃(x, y, t) is a function of space and time, describing, e.g., video signals. Signal transformations are performed by dendritic stimulus processors. They are modeled here as receptive fields whose output is encoded by spiking neurons with lateral connectivity (cross-feedback) and feedback.

The rigorous mathematical formalism employed here enables us to obtain two key results: (i) the projection of the multidimensional receptive field onto the input signal space can be perfectly identified, and (ii) the problem of identification of multidimensional receptive fields is dual to the problem of neural decoding. The duality result shows that the identification of a multidimensional receptive field is equivalent to the problem of encoding its projection with a neural circuit whose receptive field has an impulse response that is exactly the multidimensional input test signal. We provide an identification algorithm and give detailed examples arising in spiking models of audition and vision.

2. Methods

2.1. Modeling Multidimensional Stimuli and Their Processing

2.1.1. The space of input stimuli

A multidimensional communication/processing channel of interest is shown in Figure 1. An analog signal u_n(x₁, …, x_n) of n dimensions is passed through a channel with memory that describes a physical communication link or a signal processing stage. The output of the channel v is then mapped, or encoded, by an asynchronous sampler into the time sequence (t_k)_k∈ℤ. In the example shown in Figure 1, the asynchronous sampler is a (leaky) IAF neuron.

FIGURE 1

Figure 1. Multidimensional problem setting. A known multidimensional signal u_n(x₁, x₂, …, x_n), is first passed through a communication/processing channel. A non-linear sampler then maps the output v of the channel into an observable time sequence (t_k)_k∈ℤ.

We model input signals as elements of a Reproducing Kernel Hilbert Space (RKHS) (Berlinet and Thomas-Agnan, 2004). For practical and computational reasons we choose to work with the multidimensional space of trigonometric polynomials _n defined below. However, the results are not limited to this particular RKHS.

Definition 1. The space of trigonometric polynomials _n is a Hilbert space of complex-valued functions u_n = u_n (x₁, …, x_n), where

$u_n(x_1,\dots ,x_n)={\displaystyle \sum _{l_1=-L_1}^{L_1}\dots }{\displaystyle \sum _{l_n=-L_n}^{L_n}{u}_{l_1\dots l_n}}{e}_{l_1\dots l_n}(x_1,\dots ,x_n),$

defined over the domain ${\mathbb{𝔻}}_n$ = [0, T₁] × [0, T₂] × ··· × [0, T_n], where u_l1…ln ∈ ℂ and

$e_{l_1\dots l_n}(x_1,\dots ,x_n)=\frac{1}{\sqrt{T_1\cdots T_n}}\exp \text{}\left(\frac{j{l}_1{\Omega }_1{x}_1}{L_1}+\dots +\frac{j{l}_n{\Omega }_n{x}_n}{L_n}\right)\text{}.$

Here Ω_i is the bandwidth in dimension x_i, L_i is the order, and T_i = 2πL_i/Ω_i is the period, for all i = 1, …, n, and _n is endowed with the inner product 〈·,·〉 : _n × _n → ℂ

$\begin{array}{cc}\langle u_n,w_n\rangle ={\displaystyle {\int }_{{\mathbb{𝔻}}_n}{u}_n}(x_1,\dots ,x_n)\overline{w_n(x_1,\dots ,x_n)}d{x}_1\dots d{x}_n.& (1)\end{array}$

Note that given the inner product in (1), the set of elements e_{l1 … ln}(x₁, …, x_n) forms an orthonormal basis in _n. Moreover, _n is an RKHS with a reproducing kernel (RK) given by

$\begin{array}{l}{K}_n(y_1,\dots ,y_n;x_1,\dots ,x_n)=\\ {\displaystyle \sum _{l_1=-L_1}^{L_1}\text{}}\text{}\text{}\dots \text{}\text{}\text{}{\displaystyle \sum _{l_n=-L_n}^{L_n}{e}_{l_1\dots l_n}}(y_1,\dots ,y_n)\overline{e_{l_1\dots l_n}(x_1,\dots ,x_n)}.\end{array}$

In what follows we provide two- and three-dimensional signal models that arise in audition and vision. For convenience, we use a simpler and widely-used notation.

Example 1. We model spectrotemporal stimuli u₂(ν, t) as elements of an RKHS of trigonometric polynomials ₂ defined on ${\mathbb{𝔻}}_2$ = [0, T₁] × [0, T₂], where T₁ = 2πL₁/Ω₁, T₂ = 2πL₂/Ω₂, with (Ω₁, L₁) and (Ω₂, L₂), being the (bandwidth, order) pairs in the spectral direction ν and in time t, respectively. A signal u₂ ∈ ₂ can be written as

$u_2(\nu ,t)={\displaystyle \sum _{\left|l_1\right|\le L_1}{\displaystyle \sum _{\left|l_2\right|\le L_2}{u}_{l_1{l}_2}}}{e}_{l_1{l}_2}(\nu ,t),(\nu ,t)\in {\mathbb{𝔻}}_2,$

where the coefficients u_{l1 l2} ∈ ℂ and functions

$e_{l_1{l}_2}(\nu ,t)=\frac{1}{\sqrt{T_1{T}_2}}\exp \left(\frac{j{l}_1{\Omega }_1\nu }{L_1}+\frac{j{l}_2{\Omega }_{2}t}{L_2}\right)$

form an orthonormal basis for the (2L₁ + 1)(2L₂ + 1)-dimensional space ₂.

Example 2. We model spatiotemporal stimuli u₃(x, y, t) as elements of an RKHS of trigonometric polynomials ₃ defined on ${\mathbb{𝔻}}_3$ = [0, T₁] × [0, T₂] × [0, T₃], where T₁ = 2πL₁/Ω₁, T₂ = 2πL₂/Ω₂, T₃ = 2πL₃/Ω₃, with (Ω₁, L₁), (Ω₂, L₂), and (Ω₃, L₃), being the (bandwidth, order) pairs in spatial directions x and y and in time t, respectively. A visual stimulus u₃ ∈ ₃ can be written as

$u_3(x,y,t)={\displaystyle \sum _{\left|l_1\right|\le L_1}{\displaystyle \sum _{\left|l_2\right|\le L_2}{\displaystyle \sum _{\left|l_3\right|\le L_3}{u}_{l_1{l}_2{l}_3}}}}{e}_{l_1{l}_2{l}_3}(x,y,t),(x,y,t)\in {\mathbb{𝔻}}_3,$

where coefficients u_{l1 l2 l3} ∈ ℂ and functions

$e_{l_1{l}_2{l}_3}(x,y,t)=\frac{1}{\sqrt{T_1{T}_2{T}_3}}\exp \left(\frac{j{l}_1{\Omega }_{1}x}{L_1}+\frac{j{l}_2{\Omega }_{2}y}{L_2}+\frac{j{l}_3{\Omega }_{3}t}{L_3}\right)$

form an orthonormal basis for the (2L₁ + 1)(2L₂ + 1)(2L₃ + 1)-dimensional space ₃.

2.1.2. Modeling multidimensional processing

In the simplest setting, the communication/processing channel is described by a receptive field with a kernel h_n(x₁, …, x_n). The kernel is assumed to be causal in the time variable (if any) and BIBO-stable. We also assume that the kernel has a finite support of length S_i ≤ T_i in each direction x_i. In other words, each kernel h_n belongs to the space H_n.

Definition 2. The filter kernel space H_n is given by

$H_n=\left\{h_n\in {\mathbb{L}}^1({ℝ}^n)|\text{supp}(h_n)\subseteq {\mathbb{𝔻}}_n=[0,T_1]\times \cdots \times [0,T_n]\right\}.$

Since the length of the filter support is smaller than or equal to the period of the input signal in each dimension, we effectively require that for given S_i and fixed input signal bandwidth Ω_i, the order L_i of the space _n satisfies L_i ≥ S_i · Ω_i/(2π) for all i = 1, …, n.

Definition 3. The operator : H_n → _n given (by abuse of notation) by

is called the projection operator.

Since h_n ∈ _n, we have

2.2. Multidimensional [Filter]-[IAF] Tems and Their T-Transforms

In this section we analyze the general single-input single-output (SISO) multidimensional TEM shown in Figure 1 and describe in detail its I/O behavior. We then provide two specific examples of SISO multidimensional TEMs that are often encountered in neuroscience.

Without loss of generality, we assume that memory effects in the neural circuit arise only in the temporal dimension t of the stimulus and interactions in other dimensions are multiplicative in their nature. Consequently, the output v of the multidimensional receptive field is given by a convolution in the temporal dimension and integration in all other dimensions, i.e.,

$\begin{array}{l}v(t)={\displaystyle {\int }_{{\mathbb{𝔻}}_n}{h}_n}(x_1,\dots ,x_{n-1},s)\\ u_n(x_1,\dots ,x_{n-1},t-s)d{x}_1\dots d{x}_{n-1}ds.\end{array}$

The temporal signal v(t) represents the total dendritic current flowing into the spike initiation zone, where it is encoded into spikes by a point neuron model, such as the (leaky) IAF neuron of Figure 1. The mapping of the multidimensional stimulus u_n into a temporal sequence (t_k)_k∈ℤ is described by the set of equations

$\begin{array}{cc}{\displaystyle {\int }_{t_k}^{t_{k+1}}v}(t)\exp \left(\frac{t-t_{k+1}}{RC}\right)dt=q_k,k\in \mathbb{Z},& (3)\end{array}$

also known as the t-transform (Lazar, 2004; Lazar and Tóth, 2004), where

$\begin{array}{cc}{q}_k=C\delta +bRC\left[\exp \left(\frac{t_k-t_{k+1}}{RC}\right)-1\right].& (4)\end{array}$

Assuming the stimulus u_n(x₁, …, x_{n − 1}, t) ∈ _n and using the kernel representation, we have

where y = (y₁, …, y_n) and dy = dy₁ dy₂ … dy_n.

Let us now define the linear functional _k : _n → ℝ

By the Riesz representation theorem there is a function ϕ_k ∈ _n such that

We arrived at the following

Lemma 1. The SISO multidimensional TEM with a receptive field described by a kernel h_n = h_n(x₁, …, x_{n − 1}, t) provides irregular samples, or quantal measurements of the projection h_n of the kernel h_n onto the input stimulus space _n. In other words, the t-transform may be rewritten as an inner product

for every inter-spike interval [t_k, t_{k + 1}), k ∈ ℤ, where ϕ_k ∈ _n.

Remark 1. The result above has a simple interpretation. First, information about the receptive field is encoded in the form of quantal measurements q_k. These measurements can be readily computed from the spike times (t_k)_k∈ℤ. Second, the information about the receptive field is only partial and depends on the stimulus space _n used in identification. Specifically, q_k's are measurements not of the original kernel h_n but of its projection h_n onto the space H_n.

Example 3. A SISO Spectrotemporal TEM is shown in Figure 2. The signal u₂(ν, t), (ν, t) ∈ ${\mathbb{𝔻}}_2$ = [0, T₁] × [0, T₂], is the input to a communication/processing channel with kernel h₂(ν, t). The signal u₂(ν, t) may represent the time-varying amplitude of a sound in a frequency band centered around ν and h₂(ν, t) the spectrotemporal receptive field (STRF) (Kowalski et al., 1996). The output v of the kernel is encoded into a sequence of spike times (t_k)_k∈ℤ by the leaky integrate-and-fire neuron with a threshold δ, bias b and membrane time constant RC. A spectrotemporal TEM can be used to model the processing or transmission of, e.g., auditory stimuli characterized by a frequency spectrum varying in time (Kowalski et al., 1996). The operation of such a TEM can be fully described by the t-transform

$\begin{array}{cc}{\displaystyle {\int }_{t_k}^{t_{k+1}}\left[{\displaystyle {\int }_{{\mathbb{𝔻}}_2}\text{}}{h}_2(\nu ,s)u_2(\nu ,t-s)d\nu ds\right]}\exp \text{}\left(\frac{t-t_{k+1}}{RC}\right)dt=q_k,& (7)\end{array}$

with q_k given by (4) for all k ∈ ℤ.

FIGURE 2

Figure 2. Block diagram of a circuit with a spectrotemporal communication channel.

Assuming the spectrotemporal stimulus u₂(ν, t) ∈ ₂, Equation (7) can be written as

where _k : ₂ → ℝ is a linear functional. By the Riesz representation theorem (Berlinet and Thomas-Agnan, 2004), there exists a function ϕ_k ∈ ₂ such that

Example 4. A simple SISO Spatiotemporal TEM is shown in Figure 3. A video signal u₃(x, y, t), (x, y, t) ∈ ${\mathbb{𝔻}}_3$ = [0, T₁] × [0, T₂]| × [0, T₃], appears as an input to a communication/processing channel described by a filter with a kernel h₃(x, y, t). The output v of the kernel is encoded into a sequence of spike times (t_k)_k∈ℤ by the leaky integrate-and-fire neuron.

FIGURE 3

Figure 3. Block diagram of a circuit with a spatiotemporal communication channel.

A spatiotemporal TEM can be used to model the processing or transmission of, e.g., video stimuli characterized by a spatial component varying in time. The t-transform of such a TEM is given by

$\begin{array}{cc}\begin{array}{l}{\displaystyle {\int }_{t_k}^{t_{k+1}}\left[{\displaystyle {\int }_{{\mathbb{𝔻}}_3}{h}_3}(x,y,s)u_3(x,y,t-s)dxdyds\right]}\\ \exp \left(\frac{t-t_{k+1}}{RC}\right)dt=q_k,\end{array}& (9)\end{array}$

with q_k given by (4) for all k ∈ ℤ.

Assuming the video stimulus u₃(x, y, t) ∈ ₃, Equation (9) can be written as

where _k:₃ → ℝ is a linear functional. By the Riesz representation theorem, there is a function ϕ_k ∈ ₃ such that

Example 5. A special case of the SISO Spatiotemporal TEM is the SISO Spatial TEM, in which the communication/processing channel affects only the spatial component of the spatiotemporal input signal. In other words, the output of the receptive field is given by

$v(t)={\displaystyle {\int }_{{\mathbb{𝔻}}_2}{h}_2}(x,y)u_3(x,y,t)dxdy.$

Note that because only the spatial component of the input is processed, a simpler stimulus that does not vary in time may be presented when identifying this system. For example, such a stimulus may be a static image u₂(x, y). Then,

where _k : ₂ → ℝ is a functional. As before, by the Riesz representation theorem, there is a function ϕ_k ∈ ₂ such that

2.3. Identification Algorithms

In the previous section we demonstrated a relationship between the problem of identification of a receptive field and an irregular sampling problem. Namely, we showed that a projection h_n of the multidimensional receptive field h_n is embedded in the output spike sequence of the neuron as samples, or quantal measurements, q_k of h_n. A natural followup question to ask is how to reconstruct h_n from these measurements. We have the following result.

2.3.1. Feedforward multidimensional SISO CIM

Theorem 1. (SISO Multidimensional CIM)

Let {uⁱ_n | uⁱ_n ∈ _n}^N_{i = 1} be a collection of N linearly independent stimuli at the input to a [Filter]-[Leaky IAF] circuit with a multidimensional receptive field h_n ∈ H_n. If the number of signals N ≥ ∏^{n − 1}_{p = 1}(2L_p + 1) and the total number of spikes produced in response to all stimuli is greater than ∏ⁿ_{p = 1}(2L_p + 1) + N, then the filter projection h_n can be perfectly identified from a collection of input/output (I/O) pairs {(uⁱ_n, ${\mathbb{𝕋}}^i$)}^N_{i = 1} as

where h = Φ⁺ q. Here [h]_l = h_{l1, …, ln}, Φ = [Φ¹; Φ²; …; Φ^N] and the elements of each matrix Φⁱ are given by

$\begin{array}{cc}\begin{array}{l}{\left[{\Phi }^i\right]}_{kl}=\frac{RC{L}_n\sqrt{T_n}{u}_{-l_1,\dots ,-l_{n-1},l_n}^i}{j{l}_n{\Omega }_{n}RC+L_n}[e_{l_n}\left(t_{k+1}^i\right)\\ -e_{l_n}\left(t_k^i\right)\exp \left(\frac{t_k^i-t_{k+1}^i}{RC}\right)]\end{array}& (12)\end{array}$

with the column index l traversing all possible subscript combinations of l₁, l₂, …, l_n for all k ∈ ℤ, i = 1, 2, …, N. Furthermore, q = [q¹; q²; …; q^N], [qⁱ]_k = qⁱ_k and

$q_k^i=C\delta +bRC\left[\exp \left(\frac{t_k^i-t_{k+1}^i}{RC}\right)-1\right]$

for k ∈ ℤ, i = 1, …, N.

Proof: The representation (6) for stimuli uⁱ_n takes the form

with ϕⁱ_k ∈ _n. Since h_n ∈ _n and ϕⁱ_k ∈ _n, we have

and

$\begin{array}{l}{\varphi }_k^i(x_1,\dots ,x_{n-1},t)={\displaystyle \sum _{\left|l_1\right|\le L_1}\dots }{\displaystyle \sum _{\left|l_n\right|\le L_n}{\varphi }_{l_1\dots l_{n}k}^i}{e}_{l_1\dots l_n}\\ (x_1,\dots ,x_{n-1},t),\end{array}$

and, therefore,

$q_k^i={\displaystyle \sum _{\left|l_1\right|\le L_1}\dots }{\displaystyle \sum _{\left|l_n\right|\le L_n}{h}_{l_1\dots l_n}}\overline{{\varphi }_{l_1\dots l_{n}k}^i}.$

In matrix form we have qⁱ = Φⁱ h, with [qⁱ]_k = qⁱ_k, where the elements [Φⁱ]_kl = ϕⁱ_{l1 … lnk}, with the column index l traversing all possible subscript combinations of l₁, l₂, …, l_n and [h]_l = h_{l1, …, ln}. Repeating for all signals i = 1, …, N, we obtain q = Φh with q = [q¹; q²; …; q^N] and Φ = [Φ¹; Φ²; …; Φ^N]. This system of linear equations can be solved for h, provided that the rank r(Φ) of the matrix Φ satisfies r(Φ) = ∏ⁿ_{p = 1}(2L_p + 1). A necessary condition for the latter is that the total number of spikes generated by all N neurons is greater or equal to ∏ⁿ_{p = 1}(2L_p + 1) + N. Then h = Φ⁺ q, where Φ⁺ denotes a pseudoinverse of Φ. To find the coefficients ϕⁱ_{l1 … lnk}, we note that

Finally, we note that since the dendritic current v has a maximum bandwidth of Ω_t, we only need 2L_t + 1 measurements to specify it. Thus, in response to each stimulus uⁱ_n, the neuron can produce a maximum of only 2L_t + 1 informative measurements, or equivalently, 2L_t + 2 informative spikes on the interval [0, T_t]. It follows that if the neuron generates ν ≥ 2L_t + 2 spikes for each test signal, the minimum required number of signals is N = ∏ⁿ⁻¹_{p = 1}(2L_p + 1)(2L_t + 1)/(2L_t + 1) = ∏ⁿ⁻¹_{p = 1}(2L_p + 1). Similarly, if the neuron generates ν < 2L_t + 2 spikes for each signal, then the minimum required number of signals is $N=\lceil {\displaystyle {\sum }_{p=1}^n(2L_p+1)}/(\nu -1)\rceil$.

The block diagram of the identification procedure and algorithm are shown in Figure 4. Identification of the filter h_n has been reduced to the encoding of the projection h_n with a SIMO TEM whose receptive fields are uⁱ_n, i = 1, …, N.

FIGURE 4

Figure 4. Block diagram of the Multidimensional CIM. (A) Time encoding interpretation of the multidimensional channel identification problem. (B) Block diagram of the multidimensional channel identification machine.

2.3.2. SISO multidimensional CIM with lateral connectivity and feedback

Generalizing the ideas above, one can consider more complex spiking neural circuits, in which every neuron may receive not only feedforward inputs, but also lateral inputs from neurons in the same layer and back-propagating action potentials may contribute to computations within the dendritic tree. A two-neuron circuit incorporating these considerations is shown in Figure 5. Each neuron j processes a visual stimulus u^j₃(x, y, t) using a distinct spatiotemporal receptive field h^1j1₃(x, y, t), j = 1, 2. The processing of lateral inputs is described by the temporal receptive fields (cross-feedback filters) h²²¹ and h²¹², while various signals produced by back-propagating action potentials are modeled by the temporal receptive fields (feedback filters) h²¹¹ and h²²². The aggregate dendritic currents v¹ and v², produced by the receptive fields and affected by back propagation and cross-feedback, are encoded by IAF neurons into spike times (t¹_k)_k∈ℤ, (t²_k)_k∈ℤ.

FIGURE 5

Figure 5. Identifying spatiotemporal receptive fields in circuits with lateral connectivity and feedback.

Theorem 2. (SISO multidimensional CIM with lateral connectivity and feedback)

Let {[u^1,i_n, u^2,i_n] | u^j,i_n ∈ _n, j = 1, 2}^N_{i = 1} be a collection of N linearly independent vector stimuli at the input to two neurons with multidimensional receptive fields h^1j1_n ∈ H_n, j = 1, 2, lateral receptive fields h²¹², h²²¹ and feedback receptive fields h²¹¹ and h²²². Let (t¹_k)_k∈ℤ and (t²_k)_k∈ℤ be the sequences of spike times produced by the two neurons. If the number of test stimuli N ≥ ∏^{n − 1}_{p = 1}(2L_p + 1) + 2 and the total number of spikes produced by each neuron in response to all stimuli is greater than ∏ⁿ_{p = 1}(2L_p + 1) + 2(2L_n + 1) + N, then the filter projections h²¹¹, h²¹², h²²¹, h²²² and h^1j1_n, j = 1, 2, can be identified as (h²¹¹)(t) = ${\displaystyle {\sum }_{l=-L_n}^{L_n}}$ h²¹¹_l e_l(t), (h²¹²)(t) = ${\displaystyle {\sum }_{l=-L_n}^{L_n}}$ h²¹²_l e_l(t), (h²²¹)(t) = ${\displaystyle {\sum }_{l=-L_n}^{L_n}}$ h²²¹_l e_l(t) (h²²²)(t) = ${\displaystyle {\sum }_{l=-L_n}^{L_n}}$ h²²²_l e_l(t) and

Here, the coefficients h²¹¹_l, h²¹²_l, h²²¹_l, h²²²_l and h^1j1_l are given by h = [Φ₁; Φ₂]⁺ q with q = [q¹¹, …, q^1N, q²¹, …, q^2N]^T, [q^ji]_k = q^ji_k and h = [h¹; h²], where h^j = [h^1j1_{−Ln, …, −Ln}, …, h^1j1_{Ln, …, Ln}, h^{2[(j mod 2) + 1]j}_−L, …, h^{2[(j mod 2) + 1]j}_L, h^2jj_−L, …, h^2jj_L]^T, j = 1, 2, provided each matrix Φ_j has rank r(Φ_j) = ∏ⁿ_{p = 1}(2L_p + 1) + 2(2L_n + 1). The i^th row of Φ_j is given by [Φ¹ⁱ_j, Φ²ⁱ_j, Φ³ⁱ_j], i = 1, …, N, with

${[{\Phi }_j^{2i}]}_{kl}=\sqrt{T}{\displaystyle {\int }_{t_k^{ji}}^{t_{k+1}^{ji}}{t}_l^{[(j\text{mod}2)+1]i}}{e}_l(t)\exp \left(\frac{t_k^i-t_{k+1}^i}{RC}\right)dt$

and

${[{\Phi }_j^{3i}]}_{kl}=\sqrt{T}{\displaystyle {\int }_{t_k^{ji}}^{t_{k+1}^{ji}}{t}_l^{ji}}{e}_l(t)\exp \left(\frac{t_k^i-t_{k+1}^i}{RC}\right)dt,$

l = −L_n, …, L_n. The entries [Φ¹ⁱ_j]_kl are as given in Theorem 1.

Proof: Essentially the same as in Theorem 1, with an addition of lateral and feedback terms. Each additional temporal filter requires (2L_n + 1) additional measurements, corresponding to the number of bases in the temporal variable t.

3. Results

Figures 6–9, and corresponding figure legends demonstrate the performance of the multidimensional Channel Identification Machine.

FIGURE 6

Figure 6. Spectro-temporal example. Original and identified spectrotemporal filters are shown in the top and bottom plots, respectively. Ω₁ = 2π · 80 rad/s, L₁ = 16, Ω₂ = 2π · 120 rad/s, L₂ = 24.

In simulations pertaining to the spectrotemporal receptive field (see also Figure 6), we used the short-time Fourier transform of an arbitrarily chosen 200 ms segment of the Drosophila courtship song as a model of the STRF. The space of spectrotemporal signals ₂ had bandwidth Ω₁ = 2π · 80 rad/s and order L₁ = 16 in the spectral direction ν and bandwidth Ω₂ = 2π · 120 rad/s and order L₂ = 24 in the temporal direction t. The STRF appeared in cascade with an ideal IAF neuron (see Figure 2), whose parameters were chosen so that it generated a total of more than (2L₁ + 1)(2L₂ + 1) = 33 × 49 = 1, 617 measurements in response to all test signals. We employed a total of N = 40 spectrotemporal signals (which is larger than the (2L₁ + 1) = 33 requirement of Theorem 1) in order to identify the STRF.

In simulations involving the spatiotemporal receptive field (see also Figures 7, 8 we used a spatial Gabor function that was either rotated, dilated or translated in space as a function of time. The space of spatiotemporal signals ₃ had bandwidth Ω₁ = 2π · 12 rad/s and order L₁ = 9 in the spatial direction x, bandwidth Ω₂ = 2π · 12 rad/s and order L₂ = 9 in the spatial direction y, and bandwidth Ω₃ = 2π · 100 rad/s and order L₃ = 5 in the temporal direction t. The STRF appeared in cascade with an ideal IAF neuron (see Figure 2), whose parameters were chosen so that it generated a total of more than (2L₁ + 1)(2L₂ + 1)(2L₃ + 1) = 19 × 19 × 11 = 3, 971 measurements in response to all test signals. In order to identify the projection h₃ we employed a total of N = 400 spatiotemporal signals, a number that is larger than the (2L₁ + 1)(2L₂ + 1) = 361 requirement of Theorem 1).

FIGURE 7

Figure 7. Spatio-temporal example #1. Top row: Four frames of the original spatiotemporal kernel h₃(x, y, t). Here, h₃ is a spatial Gabor function rotating clockwise in space with time. Middle row: Four frames of the identified kernel. Ω₁ = 2π · 12 rad/s, L₁ = 9, Ω₂ = 2π · 12 rad/s, L₂ = 9, Ω₃ = 2π · 100 rad/s, L₃ = 5. Bottom row: Absolute error between four frames of the original and identified kernel.

FIGURE 8

Figure 8. Spatio-temporal example #1 in the frequency domain. Top row: Fourier amplitude spectrum of the four frames of the original spatiotemporal kernel h₃(x, y, t) in Figure 7. Note that the frequency support is roughly confined to a square [−10, 10] × [−10, 10]. Middle row: Fourier amplitude spectrum of the four frames of the identified spatiotemporal kernel in Figure 7. Nine spectral lines (L₁ = L₂ = 9) in each spatial direction cover the frequency support of the original kernel. Bottom row: Absolute error between four frames of the original and identified kernel.

In simulations involving the spatial receptive field (see also Figure 9), we used a static spatial Gabor function. The space of spatial signals ₂ had bandwidths Ω₁ = Ω₂ = 2π · 15 rad/s and L₁ = L₂ = 12 in spatial directions x and y. The STRF appeared in cascade with an ideal IAF neuron (see Figure 2 as a reference), whose parameters were chosen so that it generated a total of more than (2L₁ + 1)(2L₂ + 1) = 25 × 25 = 625 measurements in response to all test signals. In order to identify the projection h₂ we employed a total of N = 688 spatial signals (a number that is larger than the (2L₁ + 1)(2L₂ + 1) = 625 requirement of Theorem 1).

FIGURE 9

Figure 9. Spatial example #1. Ω₁ = Ω₂ = 2π · 15 rad/s, L₁ = L₂ = 12. A minimum of N = 625 images are required for identification. 1.1 × N = 688 images were used. (A–C) Left to right: original spatial kernel h₂(x, y), identified kernel and absolute error between the two. (D–F) Left to right: contour plots of the original spatial kernel h₂(x, y), identified kernel and absolute error. (G–I) Fourier amplitude spectrum of signals in (D–E).

Identification results for the circuit in Figure 5 are shown in Figure 10. The spatiotemporal receptive fields used in this simulation were non-separable. The first receptive field was modeled as a single spatial Gabor function (at time t = 0) translated in space with uniform velocity as a function of time, while the second receptive field was a spatial Gabor function uniformly dilated in space as a function of time. Three different time frames of the original and the identified receptive field of the first neuron are shown in Figures 10A,B, respectively. Similarly, three time frames of the original and identified receptive field of the second neuron are respectively plotted in Figures 10C,D. The identified lateral and feedback kernels are visualized in plots (e-h) of Figure 10.

FIGURE 10

Figure 10. Identifying spatiotemporal receptive fields in circuits with lateral connectivity and feedback. (A–D) Identifying the feedforward spatiotemporal receptive fields of Figure 5. (E–H) Identifying the lateral connectivity and feedback filters of Figure 5.

4. Discussion

4.1. Implications for Multidimensional Encoding and Decoding

The duality between multidimensional channel identification and stimulus decoding problems allowed us to derive identification algorithms for estimation of receptive fields of arbitrary dimensions and precise conditions under which the identification is possible. At this point it is important to pause and analyze the relationship between the dual problems. As it often turns out, looking at a dual problem can provide a tremendous insight about the primal problem.

Interestingly, previous results for video time encoding and decoding machines provided only the necessary condition of having enough spikes to decode the video (Lazar et al., 2010; Lazar and Pnevmatikakis, 2011). This condition naturally follows from having to invert a matrix in order to compute the basis coefficients of the video signal. Since the matrix needs to be full rank to provide a unique solution, and there are a total of (2L₁ + 1)(2L₂ + 1)(2L₃ + 1) coefficients involved, (2L₁ + 1)(2L₂ + 1)(2L₃ + 1) + N spikes are needed from a population of N neurons (the number of spikes is larger than the number of needed measurements by N since every measurement q_k is computed using two consecutive spikes t_k and t_{k + 1}.).

Note that a necessary condition only tells us that the number of spikes must have been greater than (2L₁ + 1)(2L₂ + 1)(2L₃ + 1) + N if we were able to recover the video signal. In order to guarantee that the video can be recovered we need a sufficient condition.

The sufficient condition can be derived by drawing comparisons between the decoding and identification problems. In identification, estimation of a receptive field from a single trial is usually not possible, even if the neuron produces a large number of spikes (Lazar and Slutskiy, 2014b). Intuitively, this is because the output of the receptive field is just a function of time. In essence, all dimensions of the stimulus are compressed into just one—the temporal dimension—and we need only (2L₃ + 1) measurements to specify a temporal function. As a result, only (2L₃ + 1) measurements are informative and we do not gain any new information if the neuron is oversampling the temporal signal. Thus, if the neuron is producing at least (2L₃ + 1) measurements per each test stimulus, we need N ≥ (2L₁ + 1)(2L₂ + 1) different trials to reconstruct a (2L₁ + 1)(2L₂ + 1)(2L₃ + 1)-dimensional receptive field. Similarly, to decode a (2L₁ + 1)(2L₂ + 1)(2L₃ + 1)-dimensional input stimulus, we need N ≥ (2L₁ + 1)(2L₂ + 1) neurons, with each neuron in the population producing at least (2L₃ + 1) measurements. If each neuron produces less than (2L₃ + 1) measurements, a larger population N is needed to faithfully encode the video signal.

More generally, if the n-dimensional input stimulus is an element of a (2L₁ + 1)(2L₂ + 1)…(2L_n + 1)-dimensional RKHS (where the last dimension is time), and the neuron is producing at least at least (2L_n + 1) + 1 spikes per test stimulus, a minimum of (2L₁ + 1)(2L₂ + 1)…(2L_{n − 1} + 1) linearly independent stimuli, or trials with linearly independent stimuli, are needed to identify the receptive field. This condition is sufficient and by duality between channel identification and time encoding, complements the previous necessary condition derived for time decoding machines.

4.2. Extensions

In the derivations above we implicitly assumed that the I/O system was noiseless. In practice, noise is introduced either by the channel or the sampler itself (Lazar et al., 2010). In the presence of noise it is not possible to identify the projection h_n loss-free. However, the analysis/methodology presented above can be extended within an appropriate mathematical setting to I/O systems with noisy measurements. For example, we can still identify an optimal estimate of h_n with respect to an appropriately defined cost function, e.g., by using the Tikhonov regularization method. The regularization methodology extensively discussed in Lazar and Slutskiy (2014a) can be adopted with minor modifications to our setting here.

Finally we note that, for convenience and simplicity of notation, the asynchronous sampler used throughout this paper was the IAF neuron. Extensions of our results to neural circuits built with other biophysically-grounded neuron models can be readily obtained by adapting the methodology developed in Lazar and Slutskiy (2012, 2014a).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by AFOSR under grant #FA9550-12-1-0232 and by NIH under grant #R021 DC012440001.

References

Aertsen, A. M., and Johannesma, P. I. (1981). The spectro-temporal receptive field. A functional characteristic of auditory neurons. Biol. Cybern. 42, 133–143. doi: 10.1007/BF00336731

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Berlinet, A., and Thomas-Agnan, C. (2004). Reproducing Kernel Hilbert Spaces in Probability and Statistics. New York, NY: Kluwer Academic Publishers. doi: 10.1007/978-1-4419-9096-9

CrossRef Full Text | Google Scholar

Brown, E., Moehlis, J., and Holmes, P. (2004). On the phase reduction and response dynamics of neural oscillator populations. Neural Comput. 16, 673–715. doi: 10.1162/089976604322860668

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Clopton, B. M., and Backoff, P. M. (1991). Spectraltemporal receptive fields of neurons of cochlear nucleus of guinea pig. Hear. Res. 52, 329–344. doi: 10.1016/0378-5955(91)90023-3

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Dayan, P., and Abbott, L. (2005). Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. Cambridge, MA: The MIT Press.

Google Scholar

DeAngelis, G. C., Ohzawa, I., and Freeman, R. D. (1995). Receptive-field dynamics in the central visual pathways. Trends Neurosci. 18, 451–458. doi: 10.1016/0166-2236(95)94496-R

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Kowalski, N., Depireux, D. A., and Shamma, S. A. (1996). Analysis of dynamic spectra in ferret primary auditory cortex. i. characteristics of single-unit responses to moving ripple spectra. J. Neurophysiol. 76, 3503–3523.

Pubmed Abstract | Pubmed Full Text | Google Scholar

Lazar, A. A. (2004). Time encoding with an integrate-and-fire neuron with a refractory period. Neurocomputing 58–60, 53–58. doi: 10.1016/j.neucom.2004.01.022

CrossRef Full Text | Google Scholar

Lazar, A. A., and Pnevmatikakis, E. A. (2011). Video time encoding machines. IEEE Trans. Neural Netw. 22, 461–473. doi: 10.1109/TNN.2010.2103323

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lazar, A. A., Pnevmatikakis, E. A., and Zhou, Y. (2010). Encoding natural scenes with neural circuits with random thresholds. Vision Res. 50, 2200–2212. doi: 10.1016/j.visres.2010.03.015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lazar, A. A., and Slutskiy, Y. B. (2010). Identifying dendritic processing. Adv. Neural Inf. Process. Syst. 23, 1261–1269.

Google Scholar

Lazar, A. A., and Slutskiy, Y. B. (2012). Channel identification machines. Comput. Intell. Neurosci. 2012:209590. doi: 10.1155/2012/209590

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lazar, A. A., and Slutskiy, Y. B. (2013). Multisensory encoding, decoding, and identification. Adv. Neural Inf. Process. Syst. 26, 3183–3191.

Google Scholar

Lazar, A. A., and Slutskiy, Y. B. (2014a). Functional identification of spike-processing neural circuits. Neural Comput. 26, 264–305. doi: 10.1162/NECO-a-00543

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lazar, A. A., and Slutskiy, Y. B. (2014b). Spiking neural circuits with dendritic stimulus processors: encoding, decoding, and identification in reproducing kernel Hilbert spaces. J. Comput. Neurosci. doi: 10.1007/s10827-014-0522-8. [Epub ahead of print]

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Lazar, A. A., and Tóth, L. T. (2004). Perfect recovery and sensitivity analysis of time encoded bandlimited signals. IEEE Trans. Circ. Syst. I 51, 2060–2073. doi: 10.1109/TCSI.2004.835026

CrossRef Full Text | Google Scholar

Pillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., Chichilnisky, E. J., et al. (2008). Spatio-temporal correlations and visual signaling in a complete neuronal population. Nature 454, 995–999. doi: 10.1038/nature07140

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

Sherrington, C. S. (1906). The Integrative Action of the Nervous System. New York, NY: C Scribner and Sons.

Google Scholar

Waters, J., Schaefer, A., and Sakmann, B. (2005). Backpropagating action potentials in neurones: measurement, mechanisms and potential functions. Prog. Biophys. Mol. Biol. 87, 145–170. doi: 10.1016/j.pbiomolbio.2004.06.009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar