Abstract

We address the problem of detecting the presence of a recurring stimulus by monitoring the voltage on a multiunit electrode located in a brain region densely populated by stimulus reactive neurons. Published experimental results suggest that under these conditions, when a stimulus is present, the measurements are gaussian with typical second-order statistics. In this letter we systematically derive a generic, optimal detector for the presence of a stimulus in these conditions and describe its implementation. The optimality of the proposed detector is in the sense that it maximizes the life span (or time to injury) of the subject. In addition, we construct a model for the acquired multiunit signal drawing on basic assumptions regarding the nature of a single neuron, which explains the second-order statistics of the raw electrode voltage measurements that are high-pass-filtered above 300 Hz. The operation of the optimal detector and that of a simpler suboptimal detection scheme is demonstrated by simulations and on real electrophysiological data.

1.  Introduction

The ability to detect the presence of a stimulus by monitoring inner brain electrical activity is of great importance for both biomedical applications and research purposes. For example, biomedical devices that replace damaged neural regions need to detect neural activity earlier in the path and subsequently decide whether to act. In brain research, detection of neural activity following delivery of a stimulus is a common research methodology. This research paradigm of presenting a stimulus and identifying (detecting) the neural response is so dominant because it is believed that the main function of the brain is to identify and react to stimuli of all kinds, and that even sophisticated behavior can be explained in terms of a response to some complex stimulus (see Brown, Kass, & Mitra, 2004). Most biomedical applications, and many research tests, require the detection to be online; that is, detection takes place in real time and as close as possible to the actual appearance of the stimulus. Furthermore, since several stimuli occur one after the other, the detection must also be sequential. The design of a sequential online detector is the main concern of this letter.

Neural signal acquisition techniques vary considerably. However, here we deal exclusively with multiunit extracellular electrical recordings of brain regions that are densely populated by neurons that react to a specific stimulus. In such recordings, the electrode is located outside the cells, and it is exposed to the activity of many cells rather than to a single cell. Apart from the electrical activity of the nearby stimulus-responsive neurons, the electrode also captures the electric field produced by the huge number of cells that are active in the background. This background noise is greatly reduced by high-pass filtering, which blocks frequencies below 300 Hz. To enable digital processing, the electrode voltage is then sampled and digitized at a rate of a few tens of KHz (in most cases 15–30 KHz) using analog-to-digital converters. Multiunit measurements have become popular for two main reasons. First, they reflect the state of a large population rather than the possibly accidental activity of a single cell. Second, these measurements are robust and can be conducted for extensive periods. The main disadvantage of this technique is the degraded signal quality.

Interestingly, multiunit measurements taken from different sensory brain regions that are known to be associated with different physical stimuli result in a very similar characteristic response: In the absence of a stimulus, the acquired signal consists mainly of noisy fluctuations that most of the time retain some baseline level. When the stimulus is presented to the subject, the acquired signal initially exhibits a noticeably higher level of noise and at times may have single spike–like fluctuations. If the stimulus is presented long enough, the noise level decays within a few tens of milliseconds to a level that is somewhat higher than that in the absence of a stimulus. As the noise level decays, the single spike–like fluctuations become less frequent. When the stimulus is removed, the noise returns to its baseline level. Figure 1 depicts an example of such a response taken from the literature.

Figure 1:

Multiunit response to a stimulus pulse reported by Holstein, Buchwald, & Schwafel, 1969a,1969b). (a) Raw data of multiunit activity in response to a stimulus. (b) The corresponding multiunit peristimulus histogram (PSTH). More recent examples can be found in Tommerdahl, Delemos, Whitsel, Favorov, and Metz (1999) and in Radner et al. (2001).

Figure 1:

Multiunit response to a stimulus pulse reported by Holstein, Buchwald, & Schwafel, 1969a,1969b). (a) Raw data of multiunit activity in response to a stimulus. (b) The corresponding multiunit peristimulus histogram (PSTH). More recent examples can be found in Tommerdahl, Delemos, Whitsel, Favorov, and Metz (1999) and in Radner et al. (2001).

The contributions of this letter are twofold. First, we propose a model that explains the acquired multiunit signal by drawing on basic assumptions on the nature of a single neuron. The model successfully predicts several well-known frequently observed empirical results that are either unexplained by other models, or demand explanations with much greater computational complexity (for a review of neuron models, see Burkitt, 2006; Herz, Gollisch, Machens, & Jaeger, 2006). Predominantly, the model predicts the well-documented sensitivity of neurons to stimuli edges and the known exponential response decay that follows the initial response. The second contribution of this letter has engineering value: we describe a generic implementation of a sequential detector for gaussian models that is optimal in many real-life scenarios, and we relate the parameters of the implementation to the requirements of the detection problem and the natural biological parameters of the electrophysiological signal. In the modeling part, we introduce a modified version of the so-called refractory point process model for a single neuron (see Ricciardi & Esposito, 1966; Teich, Matin, & Cantor, 1978; Camproux, Saunier, Chouvet, Thalabard, & Thomas, 1996; Johnson, 1996; Toyoizumi, Rad, & Paninski, 2009; Nossenson & Messer, 2010) that is computationally simpler: it has a small number of parameters, and it produces a closed-form expression connecting the firing rate to the innervating stimulus. Then we consider the implications of this model for situations where a large number of neurons are involved, that is, a situation where the central limit theorem holds. We show that these two features force the covariance matrix of the resulting gaussian signal to have a specific structure. In addition, we show how to present this covariance matrix in state space. This procedure has technical importance for both detection and estimation algorithms. In the detector construction part, we customize ideas from the pioneering work of Schweppe (1965) for sequential detection of multiunit neural measurements. More specifically, we incorporate knowledge regarding the structure of the covariance matrix with prior knowledge regarding stimuli timing, which is suitable for real-life scenarios as well as for laboratory experiments.

Starting with Weber and Buchwald (1965), the detection schemes that are typically employed in conditions of multiunit activity consist of squaring the signal, integrating, and thresholding. This heuristic method and its variants (see, e.g., Basseville & Nikiforov, 1993) are based on good common sense. Nevertheless, the analytic detection scheme proposed in this letter has two key advantages over these methods: (1) the analytic construction guarantees optimal performance, whereas the heuristic solutions are suboptimal, and (2) the parametric construction of the detector introduced here facilitates well-defined analytic procedures for tuning the detector to the measurement conditions under well-defined quality criteria, whereas the heuristic detectors require a long trial-and-error tuning stage that does not ensure optimality. More recent studies focused on the so-called spike sorting problem (see the review by Lewicki, 1998). The underlying assumption of the spike sorting problem is that the stimulus reactive source and the unrelated spiking sources produce different spike shapes, making it possible to distinguish one type of source from the other. However, in the conditions discussed here, it is impossible to isolate and identify single spikes with reasonable probability because of the ambiguity of observations that result from the cumulative action of a large number of spiking sources embedded in thermal noise.

The work here is also closely related to papers dealing with biological neural networks models. In the detector construction section, our methodology relies heavily on the work of Schweppe (1965) whose techniques were recently used in the context of processing neural signals by Roweis and Ghahramani (1999), Barbieri et al., (2004), Paninski, Pillow, and Lewi, (2007), Beck and Pouget (2007), and Deneve (2008). The differences between these works and the work here are mainly related to the choice of neuron model and the problem formulation. Specifically, we consider the detection problem for a scenario in which stimuli are admitted in blocks and are separated by exponentially distributed intervals.

The rest of this letter is organized as follows. In section 2, we mathematically formulate the detection problem. In section 3, we explain the biophysical origin of the multiunit signal and calculate its statistical properties. In section 4, we derive the equivalent state-space model. In sections 5 and 6, we describe in detail the structure of the optimal detector. In section 7, we analyze the detector performance using simulations and also test it on real electrophysiological data. We summarize and conclude in section 8.

2.  Problem Formulation

In this section we formulate the detection problem using a Bayesian approach. We consider a scenario in which a multiunit electrode is implanted in a densely populated brain region that is sensitive to a certain warning stimulus. The detector monitors the multiunit data and must decide within a very small delay, Td, whether a warning signal is present. If the detection is too late, an aversive stimulus hurts the subject.

The durations of the intervals between the onsets of two consecutive warning signals are independent, and identically distributed (i.i.d.) random variables. The lth intertest-interval interval is designated as TITI(l) and obeys the following probability distribution:
formula
2.1
where Tmin is the shortest possible interval between adjacent stimuli and Tavg is the average interval between stimuli. Both parameters are assumed to be known.

The warning signal and the corresponding neural activity continue until the offset of the aversive signal. The total duration of the warning signal Ts is fixed and is known a priori. Note that if the brain region is responsive to the aversive stimulus rather than to the warning signal, then our formulation holds, with Td equal to zero.

While the aversive stimulus is being delivered, the subject is at risk of injury that is proportional to the stimulus intensity. Here we assume that the intensity of the stimulus is known a priori; hence the intensity is already embedded in the momentary risk.

It is possible to take an action that eliminates the risk of that specific injury, but at the same time, this action increases the risk from other threats. For example, if the aversive stimulus is an air blow to the eye, shutting the eyelid can prevent injury. However, a subject with shut eyes is more exposed to other threats. This example illustrated in Figure 2.

Figure 2:

Illustrative warning stimulus, aversive stimulus, and acquired data waveforms, together with false alarm and misdetection costs. In cases where only one stimulus exists, the formulation remains the same, and Td signifies the allowed delay for asserting the detection flag.

Figure 2:

Illustrative warning stimulus, aversive stimulus, and acquired data waveforms, together with false alarm and misdetection costs. In cases where only one stimulus exists, the formulation remains the same, and Td signifies the allowed delay for asserting the detection flag.

The goal of the automated decision system (henceforth referred to as the detector) is to keep the subject unharmed for as long as possible by continuously minimizing the risk of injury, ℜinjury(t), at any given moment.

To reach a decision, the detector is allowed to use only data that have been acquired so far by the multiunit electrode. This data are designated as rt0 (the electrode voltage from start time to current time).

The momentary risk of an injury, ℜinjury(t), depends on the detector's decision whether to declare that the stimulus is on, Ion(t) = 1, and risk a false alarm, or declare that the stimulus is off Ioff(t) = 1 and risk a misdetection:
formula
2.2
The momentary probability for an injury in an infinitesimal time interval, [t, t + Δt], is proportional to the momentary risk times the decision spacing interval:
formula
2.3
It follows that the probability of staying uninjured at time t = N · Δt (i.e., N time intervals from time axis origin) is
formula
2.4
Thus, at the limit Δt → 0, the probability of staying uninjured can be expressed as
formula
2.5
The goal of the detector is to decide whether to risk misdetection Ioff(t) = 1 or a false alarm Ion(t) = 1 so that the total risk is minimized and the expected life span of the subject is maximized:
formula
2.6
where, E{Tuninjured} is the expectancy of the first injury time:
formula
2.7

The considerations of the detector are illustrated in Figure 3.

Figure 3:

Detector role. Automated decision system for reduction of risk.

Figure 3:

Detector role. Automated decision system for reduction of risk.

We now clarify a few subtle points:

  • • 

    injury(t), CMD(t), and CFA are given in units of 1/sec or Hz.

  • • 

    E{} Stands for ensemble averaging of the costs. The ensemble-averaged costs and their corresponding risks, ℜMD(t) and ℜFA(t), are given in units of Hz.

  • • 

    Ion(t) and its complement are indicator functions for the decision of the detector. Ion(t) = 1 means that at time t, the detector preferred the risk of false alarm and therefore asserts the detection flag.

Table 1 summarizes the quantities introduced in this section.

Table 1:
Random Variables, Parameters, and Observables of the Problem.
MarkMeaningUnit
TITI Intertrial interval. A random variable that represents interval duration between consecutive onsets of aversive stimuli. sec 
Tmin Shortest possible interval between consecutive onsets of the aversive stimuli. sec 
Ts Duration of the neural response to the warning signal. sec 
Td Delay of the onset of the aversive stimulus from the onset of the warning signal. sec 
CMD The momentary cost for misdetecting the aversive stimulus while it is on. Hz 
CFA The momentary cost for falsely alarming on the presence of an aversive stimulus. Hz 
rt0 Voltage on the multiunit electrode, r(t) from start time until time t. These are the observables of the problem. Volts 
Δt Sampling interval. sec 
MarkMeaningUnit
TITI Intertrial interval. A random variable that represents interval duration between consecutive onsets of aversive stimuli. sec 
Tmin Shortest possible interval between consecutive onsets of the aversive stimuli. sec 
Ts Duration of the neural response to the warning signal. sec 
Td Delay of the onset of the aversive stimulus from the onset of the warning signal. sec 
CMD The momentary cost for misdetecting the aversive stimulus while it is on. Hz 
CFA The momentary cost for falsely alarming on the presence of an aversive stimulus. Hz 
rt0 Voltage on the multiunit electrode, r(t) from start time until time t. These are the observables of the problem. Volts 
Δt Sampling interval. sec 

3.  Signal Model

In this section we describe the statistical nature of the acquired signal rt0. In first section 3.1, we define the exact measurement conditions for which this detector is intended. In section 3.2, we explain that under these conditions, the acquired signal rt0 is a gaussian process with a specific structure of covariance function. This property is important for detector construction because the optimization process of any detector relies on the statistical nature of the received signal. A summary of the signal model is given in section 3.3.

3.1.  Expected Measurement Conditions.

The measurement conditions that are suited for the proposed detector are as follows:

  • • 

    The measurement is extracellular, and the electrode is implanted in a brain region that is densely populated with stimulus reactive cells. The emphasis is on stimuli reactive cells because any inner brain extracellular measurement is exposed to a large population of nerve cells (Gold, Henze, Koch, & Buzsaki, 2006; Milstein & Koch, 2008). More quantitatively, we require that the surrounding of the electrode, up to a radius of ∼150μ (Gold et al., 2006; Milstein & Koch, 2008) be densely populated with stimulus-responsive cells. Thus, in a region with a cell density of 10,000 cells/mm3 (Weber, Chen, Hubbard, & Kaufman, 2000), the total number of cells within the required proximity is on the order of many tens to a few hundred.

  • • 

    The operating mode of the neurons in the brain region of interest is a refractory Markov neuron mode (see Ricciardi & Esposito, 1966; Teich et al., 1978; Johnson, 1996; Toyoizumi, Rad, & Paninski, 2009; Nossenson & Messer, 2010). The refractory Markov mode is quite common in sensory brain regions, and the decision whether the neurons in the region of interest are indeed in this mode is not difficult, because such neurons produce an unmistakable response pattern, which is given in Figure 5. Examples of such a response can also be found in Creutzfeldt, Hellweg, and Schreiner (1980), Tsumoto, Creutzfeldt, and Legendy (1978), Caird and Klinke (1983), Aston-Jones and Bloom (1981), and Di Lorenzo and Schwartzbaum (1982).

  • • 

    The received signal is bandpass filtered between 300 Hz and 4000 Hz.

3.2.  The Resulting Signal Model.

The measurement conditions that were specified in section 3.1 dictate a very specific signal model. The received signal, r(t), results from the cumulative action of three processes:
formula
3.1
where
  • • 

    x(t) is the electrode voltage that originates from the stimulus reactive spiking neurons.

  • • 

    n(t) is the noise due to the electrical activity of the background neurons unrelated to the stimulus.

  • • 

    v(t) is noise from the electrical measurement equipment.

Because both the desired signal and the background noise originate from large populations, the central limit theorem holds, and both signal and noise are gaussian processes:
formula
3.2
formula
3.3
formula
3.4
As a result, the total electrode voltage is also a gaussian process:
formula
3.5
The gaussian nature of the multiunit measurement can be verified by plotting an amplitude histogram of the process as in Schwartz, Ramos, and John (1976), Heetderks (1978), and Legatt, Arezzo, and Vaughan (1980).

In the case of a gaussian process, the statistical properties of the received signal are completely determined by the mean, μ(t), and the covariance function, Σ(t1, t2). Moreover, under the specified measurement conditions, each of the components that comprise the received signal has a specific structure of covariance function. We specify these in the next subsections.

3.2.1.  The Mean and Covariance of the Desired Signal.

As stated in section 3.1, we require that the stimulus-responsive neurons are refractory Markov neurons. The simplest form of such a neuron is best described by the two-state diagram depicted in Figure 4. Under this model, the neuron is either armed and ready to fire, or it is in a refractory state in which it cannot fire. The transition from the armed state to the refractory state is accompanied by spike firing. When the neuron is armed, its probability of firing a spike in an infinitesimal time segment [t, t + Δt) is
formula
3.6
where y(t) is proportional to the short-term average intensity of the stimulus and the constant R0 is the average spontaneous firing rate of the neuron in the armed state. Note that y(t) corresponds to the amount of excitatory chemicals that are produced by either the various senses or an external source (added artificially). After a spike is fired, the neuron starts a recovery process that has a fixed rate (independent of stimulus intensity), R1.
Figure 4:

Modeling a neuron as a two-state Markov chain.

Figure 4:

Modeling a neuron as a two-state Markov chain.

The equations describing the two-state model are as follows:
formula
3.7
formula
3.8
By substituting P1 = 1 − P0 in equation 3.7 and then using the integration factor method, we find that the probability of staying in an armed state in any moment as a function of stimulus intensity through time is
formula
3.9
where t0 is the start time and P0(t0) is the probability of being in an armed state at the start time.
The corresponding firing rate of this model is
formula
3.10

Figure 5 depicts the shape of the function Rfire(t) when the stimulus y(t) is a rectangular pulse.

Figure 5:

(a) The expected firing rate of a refractory Markov neuron in response to a rectangular stimulus according to equations 3.9 and 3.10. (b) Single-cell, intracellular measurements taken in the auditory cortex and by Wang, Lu, Snider, and Liang (2005). Other examples can be found in Taberner and Liberman (2005), Ingham and McAlpine (2005), and Reches and Gutfreund, (2008).

Figure 5:

(a) The expected firing rate of a refractory Markov neuron in response to a rectangular stimulus according to equations 3.9 and 3.10. (b) Single-cell, intracellular measurements taken in the auditory cortex and by Wang, Lu, Snider, and Liang (2005). Other examples can be found in Taberner and Liberman (2005), Ingham and McAlpine (2005), and Reches and Gutfreund, (2008).

When an ensemble of refractory Markov neurons is fed with the same stimulus, the spike firing probability of each neuron in the ensemble is altered according to equation 3.10. The extracellular electrode then captures the spatial sum of the spikes arriving from all sources. It is known from the empirical measurements of Logothetis, Kayser, and Oeltermann (2007) and also from the analytical work of Gold et al. (2006) that the spikes arriving at the electrode are deformed versions of an intracellularly measured spike and that the deformation depends on the position of the spiking source relative to the electrode. Furthermore, these studies show that the deformation is limited to a very specific subfamily of deformations: the spikes are attenuated with the distance, and they can be phase-shifted. As we explain, this type of deformation has important implications for the mean response and the variance of the response.

The mean response captured by the electrode, μx(t), is the sum of the expected responses from all individual spiking sources. However, because it is difficult to know a priori the locations of the somas relative to the electrode and because each location leads to a different phase shift, the mean response of an ensemble of refractory Markov neurons, μx(t), may appear in several forms, as shown in Figure 6b. For this reason, the mean response to a stimulus, μx(t), should be found empirically using a grand average method. In other words, the empirical mean is calculated by averaging several responses to a warning stimuli.

Figure 6:

(a) The variance of the response to a rectangular stimulus as a function of time has a characteristic shape that resembles the firing rate curve. (b) The mean response to stimulus μx(t) may come in several forms depending on electrode location, and it should be measured empirically. In the top figure, the mean is clearly seen at the edges, whereas in the lower figure, the mean is zero.

Figure 6:

(a) The variance of the response to a rectangular stimulus as a function of time has a characteristic shape that resembles the firing rate curve. (b) The mean response to stimulus μx(t) may come in several forms depending on electrode location, and it should be measured empirically. In the top figure, the mean is clearly seen at the edges, whereas in the lower figure, the mean is zero.

Unlike the mean response, the covariance function of an ensemble of refractory Markov neurons, Σx(t1, t2), will always have a shape that resembles the firing rate regardless of the electrode location,
formula
3.11
where:
  • • 

    gx(t) is a function that has an amplitude frequency response of a spike that is high-pass-filtered above 300 Hz and an arbitrary phase response. An example of such a spike shape is given in Figure 7. In practice, the exact spectral shape of the spike can also be measured empirically by calculating the periodograms before the response and during the steady-state part of the response. The spectrum increase is due to the stimuli-responsive neurons.

  • • 
    The operation ∫ψ(t1, t2, τ)○dτ is defined as
    formula
    3.12
    This operation consists of a regular integral and an additional correction term. The correction term originates from the fact that the observed random processes are all causal (each spike contributes only after it has started), and they are band-limited either naturally (the spikes) or because of the external low-pass filtering (which limits the thermal noise). The causality and the fact that these processes are band-limited and, hence, smooth introduce a small, yet accumulating correction term to the autocorrelation relative to the simple expression that is obtained without considering the effect of these two properties. An important property of this operation is that the resulting expression for the covariance is phase insensitive, so that any spike shape that is phase-shifted will yield that same covariance function. Recall that the studies by Gold et al. (2006) and Logothetis et al. (2007) indicate that the signal captured by the electrode is a sum of phase-shifted versions of a regular spike. Because the covariance is phase insensitive, it is expected to have the same form regardless of the phase shift of the spiking sources. As shown later, the phase invariance also relaxes the rules for hardware implementation.
Figure 7:

(a) A Hodgkin-Huxley spike (taken from Milstein & Koch, 2008) after high-pass filtering at 300 Hz. (b) Magnitude of the Fourier transform of the curve in a.

Figure 7:

(a) A Hodgkin-Huxley spike (taken from Milstein & Koch, 2008) after high-pass filtering at 300 Hz. (b) Magnitude of the Fourier transform of the curve in a.

The shape of the variance, Σx(t, t), in the case of a rectangular stimulus is shown in Figure 6a.1 In addition, it is possible to further validate equation 3.10 empirically by plotting a three-dimensional graph of the empirical covariance . Equation 3.10 predicts that the covariance is concentrated in the vicinity of the diagonal t1 = t2.

3.2.2.  Covariance of the Background Signal.

The background noise, n(t), is the result of the huge number of neurons that are electrically active regardless of stimulus presence. It is empirically known (Watkins, Santhanam, Shenoy, & Harrison, 2004; Martinez, Pedreira, Ison, & Quian Quiroga, 2009; Novak et al., 2009; Hill, Mehta, & Kleinfeld, 2011) that when this background noise is high-pass-filtered above 300 Hz, the resulting noise is gaussian and stationary, and it has a characteristic and reproducible spectral shape of a filtered spike. Thus, the expression for the covariance stays the same; only the stimulus intensity is set to zero. We are left with the spontaneous firing rate of the irrelevant neurons:
formula
3.13

Note that due to stationarity, the exact spectral shape can also be measured empirically by calculating the periodogram.

3.2.3.  Additive White Noise.

The received signal also contains a weak white noise component caused by the electrical equipment and the residuals of the filtering. The covariance function of this noise component is
formula
3.14
where σ2v is the value of the noise floor of the empirical spectrum at the right end of the frequency axis.

3.3.  Summary of the Model.

The multiunit electrode is exposed to many neurons. Some of the neurons react to the stimulus and are considered the signal. Other neurons are unrelated to the stimulus and fire spontaneously. We refer to these neurons as noise. Due to the large number of cells involved, both noise and signal obey the central limit theorem and have gaussian distributions. Moreover, these distributions have a specific structure of covariance function that is described by equations 3.11 to 3.14. The covariance intensity depends on the stimuli intensity through the firing rate, which is given by equations 3.9 and 3.10. The signal mean also depends on the firing rate and the electrode neuron channel, but unlike the covariance function, it is phase sensitive and can appear in several forms depending on the electrode location relative to the cells. Table 2 summarizes the model.

Table 2:
Summary of Model.
r(t) = x(t) + n(t) + v(t
x(t) ∼ Nx(t), Σx(t1,t2)) 
n(t) ∼ N(0, Σn(t1,t2)) 
v(t) ∼ N(0, Σv(t1,t2)) 
 
 
Σv(t1,t2) = σ2v○δ(t2t1
The mean μx(t) depends on brain-electrode geometry and has to be measured empirically. 
A coarse estimate of the firing rate Rfire(τ) can be extracted empirically by performing a short-lasting single-cell PSTH (peristimulu time histogram) measurement or by using parametric methods, that is, by estimating the parameters R0, R1, and then calculating Rfire(τ) using equations 3.9 and 3.10
The parameters σ2n, σ2v are found empirically. 
gx(t) is a function that has the same magnitude of the Fourier transform of a high-pass-filtered spike produced by the stimulus-responsive cell and an arbitrary phase response. An example of such a function is given in Figure 7, and a possible digital implementation of this function is described in Table 4
gn(t) is a function that has the same magnitude of Fourier transform of a high-pass-filtered spike produced by the unrelated cells and an arbitrary phase response. 
r(t) = x(t) + n(t) + v(t
x(t) ∼ Nx(t), Σx(t1,t2)) 
n(t) ∼ N(0, Σn(t1,t2)) 
v(t) ∼ N(0, Σv(t1,t2)) 
 
 
Σv(t1,t2) = σ2v○δ(t2t1
The mean μx(t) depends on brain-electrode geometry and has to be measured empirically. 
A coarse estimate of the firing rate Rfire(τ) can be extracted empirically by performing a short-lasting single-cell PSTH (peristimulu time histogram) measurement or by using parametric methods, that is, by estimating the parameters R0, R1, and then calculating Rfire(τ) using equations 3.9 and 3.10
The parameters σ2n, σ2v are found empirically. 
gx(t) is a function that has the same magnitude of the Fourier transform of a high-pass-filtered spike produced by the stimulus-responsive cell and an arbitrary phase response. An example of such a function is given in Figure 7, and a possible digital implementation of this function is described in Table 4
gn(t) is a function that has the same magnitude of Fourier transform of a high-pass-filtered spike produced by the unrelated cells and an arbitrary phase response. 

4.  State-Space Representation of Signal Model

In this section, we show that the signal model summarized in Table 2 has an equivalent state-space representation. The equivalence is in the sense that both the state-space model and the model depicted in Table 2 have the same mean and covariance functions, and since both models are gaussian, they will have identical statistical properties. We present the state-space model because this formulation is necessary to construct a sequential online detector.

We start the construction of the state-space model by repeating equation 3.11:
formula
4.1
Now, ignore for a moment equation 4.1 and consider the following linear system, which is driven by a varying-amplitude white noise:
formula
4.2
The covariance of x(t) is
formula
4.3
Observe that equations 4.1 and 4.3 would be identical if we set f(t) as the square root of the firing rate of a refractory Markov neuron:
formula
4.4

At this point, we note that the integral in equation 4.2 and its corresponding autocorrelation function, equation 4.3, are known as the Stratonovich integral (Stratonovich, 1966).

Next, we present the spike response, gx(t), using a recursive infinite impulse response (IIR) filter:
formula
4.5
Figure 8 compares the impulse response of an IIR filter with Qx,a = 6, Qx,b = 4 to an unfiltered Hodgkin and Huxley spike (H&H) taken from Gold et al. (2006). This comparison is for illustration alone, as here we discuss only the case of a filtered spike. Table 4 in appendix A lists the coefficients of a filtered version of this spike, which is shown in Figure 7.
Figure 8:

Comparison of spike shape versus state-space approximation with Qx,a = 6 and Qx,b = 4 in the time domain (a) and in the Fourier frequency domain (b). The spike shape was taken from Milstein and Koch (2008).

Figure 8:

Comparison of spike shape versus state-space approximation with Qx,a = 6 and Qx,b = 4 in the time domain (a) and in the Fourier frequency domain (b). The spike shape was taken from Milstein and Koch (2008).

Recall that equations 3.11 and 3.12 imply that the phase response of Gx(s) is immaterial for the expression of the covariance function, so there exists some degree of freedom for choosing the coefficients bx,m.

We have shown that the following autoregressive moving average (ARMA) equation yields the same covariance matrix of the process x(jΔt):
formula
4.6

Note that by adding the nonzero mean, μx(t), the ARMA equivalent of x(t) also has the same mean.

Similarly, the covariance Σn(t1, t2) associated with the background spontaneous spikes is obtained if n(t) obeys the following stochastic difference equation:
formula
4.7
Thus, r(kΔt), which is the sum of the uncorrelated random processes x(jΔt), n(jΔt) and v(jΔt), can be expressed using a sum of two ARMA processes that rely on data from the past and additional white noise components from the most recent observation (innovations):
formula
4.8
where ur[jΔt] is the process mean and dw[jl], dζ[jl], and dξj are independent standard white noise processes.
The final step in the construction of the state space is to convert the two difference equations for x[j] and n[j], which are on the order of {Qx,a, Qx,b} and {Qn,a, Qn,b}, respectively, to a set of (Qx,max + Qn,max) first-order difference equations where
formula
4.9
formula
4.10
This type of representation is known as a state-space representation. Note that there are many possible state-space representations, and all are acceptable as long as the equation for the observations, 4.8, either stays the same or has the same covariance function. We chose the following common set of equations:
formula
4.11
formula
4.12
where
  1. is a vector of internal state variables of size (Qx,max + Qn,max) × 1. The first entry of this vector, , corresponds to x(kΔt) minus the latest innovation of x(kΔt). The following (Qx,a − 1) entries are used for the construction of . Entry number Qx,a + 1 of this vector, , corresponds to n(kΔt) minus the latest innovation of n(kΔt). The following (Qn,a − 1) entries are used for the construction of .

  2. The matrices and are given in appendix A. They generate and .

  3. The vector sums the state variables and and the vector adds the latest innovations that were missing (see step 1). The values of and are also given in appendix A.

5.  Optimal Detection of Presence of a Stimulus

In this section, we develop the optimal detector for the problem formulated in section 2.

Recall that the role of the detector is to minimize the risk of an injury given the observations
formula
5.1
The risk, ℜinjury(t), is the average cost of associating the current moment with the wrong group, that is, declaring on the current moment as a moment with the stimulus is present in a case where there is no stimulus or vice versa. Since our cost functions are piecewise constants, the average risk can be expressed as follows:
formula
5.2
where PFA(tTFA | rt0) is the probability that the current time falls inside the false alarm zone and PMD(tTMD | rt0) is the probability that the current time is within the detection window.
From equation 5.2, it is evident that in order to minimize the risk of an injury, ℜinjury(t), the following decision rule should be applied:
formula
5.3
Thus, we raise the detection flag if the risk of misdetection is larger than the risk of a false alarm. To reach a decision, the detector must calculate the conditional probabilities PFA(tTFA | rt0) and PMD(tTMD | rt0) and set the output of the detector according to equation 6.1. Since the duration of the stimulus is a priori known, asking whether the stimulus is active more than Td seconds (this is the legal delay from onset) is equivalent to asking whether the elapsed time from the most recent stimulus onset is greater than Td and smaller than Ts (stimulus duration). Phrasing the detection problem this way is more convenient because two different hypotheses regarding the location of stimulus onset are mutually exclusive events. This makes it possible to calculate the total probability of an onset presence in one of several possible locations by simply summing the corresponding probabilities of these different onset locations. Defining
formula
5.4
Then the probabilities PMD(tTMD | rt0) and PFA(tTFA | rt0) can be expressed as the sum of the onset location probabilities:
formula
5.5
formula
5.6
where t is the current time and is the moment of last onset.
To avoid using an infinite sum of onset locations, we define
formula
5.7
Using this definition, equation 5.6 can be written as
formula
5.8
Furthermore, based on principles from the pioneering work of Schweppe (1965), we show in appendix B that the posterior probabilities, and , of the specific problem that was phrased in section 2 can be calculated in a recursive manner; that is, the values of these probabilities at time t are an update of their previous set of values at time t − Δt using the most recent observation.2

In section 6, we describe the implementation of the detector that sequentially calculates the probabilities for each possible delay from the most recent onset. The detector raises the detection flag if the sum of the probabilities corresponding to onset delays greater than Td and smaller than Ts is larger than the sum of probabilities of the remaining alternatives.

6.  Detector Structure

In this section, we describe in detail the implementation of the detector. The design is a customization of the work of Schweppe (1965) for sequential detection of multiunit neural signals. The applicability of the general ideas of Schweppe (1965) for processing neural signals has been discussed in Roweis and Ghahramani (1999), Barbieri et al. (2004), Paninski et al. (2007), Beck and Pouget (2007), and Deneve (2008), but the exact equations here are different.

The basic idea of the detector is to sequentially calculate normalized versions of the posterior probabilities, and and set the detection flag according to the following equivalent version of decision rule 5.3:
formula
6.1
The detector is implemented digitally, and thus the timescale is discrete:
formula
6.2
where Δt is the duration of the clock cycle:
formula
6.3
The possible interval lengths from the last onset are also discretized, and thus the detector tests only the following discrete and finite set of events:
formula
6.4
and also the event:
formula
6.5
Thus, Mmin is set as
formula
6.6
Note that according to the problem formulation, the event (onset is imminent) is contained in the event .

As shown in Figure 9, the detector consists of three main substructures:

  • • 

    The register array which stores the normalized posterior probabilities

  • • 

    The Kalman (Kalman, 1960) estimators array

  • • 

    The detection flag assertion unit

Figure 9:

Detector structure. The Kalman array of estimators, the register array, and the decision unit.

Figure 9:

Detector structure. The Kalman array of estimators, the register array, and the decision unit.

We describe these three subunits in more detail in the following subsections.

6.1.  The Register Array.

The register array stores logarithmic versions of Mmin normalized posterior probabilities. Each of the posterior probabilities has a corresponding register designated as dk[jΔt]. The first register, d1[jΔt], corresponds to the posterior probability:
formula
The second register, d2[nΔt], corresponds to the posterior probability,
formula
and so on. The last register, , corresponds to the posterior probability:
formula
Rather than storing the the posterior probability itself, the value that is stored in each register at time t = jΔt is the following function of the posterior probability:
formula
6.7
As a result of this normalization, register #Mmin is constantly kept at zero:
formula
6.8
The registers are updated every clock cycle (every Δt seconds) following the arrival of new sampled data from the electrode. The update rules of the registers are designed such that the new value of the register will indeed reflect the logarithm of the normalized posterior probability. Here we only summarize the update rules that follow from direct calculation of the posterior probabilities, which is given in appendix B. To phrase the update rules in a compact manner, we use the following short notations:
formula
6.9
and
formula
6.10
The parameter α is the reciprocal of the mean number of clock cycles added to the minimal intertrial intervals (Tmin), and it is a known a priori. The value of qk[nΔt] is calculated by the Kalman estimators array (see section 6.2) using the method in Schweppe (1965). Here, assume that qk[nΔt] is known.

Using the above notations, the register update rules can be stated as follows:

  1. Rule 1. The update rule for the first register, d1[k], is:
    formula
  2. Rule 2. The update procedure of the other (Mmin − 2) registers, is
    formula
  3. Rule 3. Lmin[j] is defined as follows:
    formula

Proof. See appendix B.

6.2.  The Kalman Estimators Array.

The array of Kalman estimators consists of Mmin units. Each unit outputs the quantity qk(jΔt), which is the logarithm of the probability of the current observation conditioned by past observations and the hypothesis that the stimulus occurred kΔt seconds ago, that is,
formula
6.14
where r[jΔt] is the most recently received electrode sample (current time is t = jΔt):
  • • 

    The first Kalman estimator, unit 1, is designed under the assumption that stimulus onset occurred Δt seconds ago: .

  • • 

    The second Kalman estimator, unit 2, is designed under the assumption that stimulus onset occurred 2Δt seconds ago: .

  • • 
    The last Kalman estimator, unit #Mmin, is designed under the assumption that no stimulus has occurred in the last Tmin seconds:
    formula
  • • 

    All the Kalman estimators receive a mean free version of the signal under the corresponding hypothesis. In other words, the mean of the corresponding hypothesis, μx[k], is removed (see Figure 6b, which shows examples for signal mean).

A schematic description of the different hypotheses associated with each estimator appears in Figure 10.

Figure 10:

The different hypotheses of the Kalman estimators. Each Kalman estimator performs the Kalman recursion using a different hypothesis regarding the current variance. The first Kalman estimator always assumes that stimulus onset occurs now, and therefore it is designed for increased variance. The second Kalman filter assumes that onset occurred on the previous sample, and so on. The last filter assumes no stimulus. See also Figure 5.

Figure 10:

The different hypotheses of the Kalman estimators. Each Kalman estimator performs the Kalman recursion using a different hypothesis regarding the current variance. The first Kalman estimator always assumes that stimulus onset occurs now, and therefore it is designed for increased variance. The second Kalman filter assumes that onset occurred on the previous sample, and so on. The last filter assumes no stimulus. See also Figure 5.

The Kalman estimators are chained (see Figure 9). The purpose of the chaining is to transfer the internal state variable of one estimator to the next. In this way, the kth Kalman estimator always maintains the same hypothesis regarding offset from stimulus although time progresses (similar to a moving assembly line). In addition to the state vector, the filters also pass on the expected estimation error covariance matrix . The latter, however, could be computed at the offline parameter estimation stage and set as a constant after the initial transient period of online detection is over.

Every clock cycle (every Δt [sec]), each one of the Mmin estimators receives three new inputs and produces three new outputs. The inputs of the kth Kalman filter at time t = jΔt are the mean free version of the current electrode voltage, r[jΔt] − ux[k], and two additional inputs that originate from estimator #(k-1): and . Each filter executes once the recursion that appears in equations C.1 to C.5 of appendix C, and passes on new values of the qk[j] to the register array, and new values of to the next estimator, which is estimator #(k+1). Note that the recursions executed by different estimators are different because they incorporate different values of and (see the state space section, equation A.3). Estimator #k uses and at all times.

The first Kalman estimator (#1) and the last Kalman estimator (#Mmin) are an exception. The first estimator receives the inputs from Kalman estimator #Mmin. The last Kalman estimator (#Mmin) produces two different sets of predictions:

  1. . This prediction is sent to Kalman estimator #1. It corresponds to the hypothesis that the stimulus will occur in the next sample.

  2. . This prediction is sent back to the Kalman estimator #Mmin. It corresponds to the hypothesis that no stimulus will occur in the next sample.

At the initialization of online stage (t = 0), the values of are set differently:
formula
6.15
The procedure for calculating is given in appendix C.

6.3.  Assertion of the Detection Flag.

Recall from section 5 that the detection flag is asserted according to the following decision rule:
formula
6.16
Using equation 5.8, we express the numerator of equation 6.16 as follows:
formula
6.17
Using equation 5.5, we express the denominator of equation 6.16 as follows:
formula
6.18
which yields this decision rule:
  1. Rule 4. The decision flag is raised according to the following condition:
    formula
    where
    formula
    is the legal delay from the onset of the warning signal expressed as number of clock cycles and
formula
6.21

is the number of clock cycles in which the warning signal persists.

7.  Results

7.1.  Simulation Study.

In this section we study the operation of the detector by means of simulations. The ultimate goal is to quantify the probability of misdetection () and the probability of false alarm (PFA[k]), both as functions of time, and consequently estimate the expected life span of the subject (or expected time to injury; see Figure 3) using the following relations:
formula
7.1
where
formula
7.2
and
formula
7.3

Since analytic evaluation of PFA(t) and PMD(t) is a formidable task, we evaluated PMD(t) and PFA(t) using simulations.3

In our simulations, we generated a gaussian signal with a time-dependent autocorrelation function as described in sections 3 and 4. The complete list of the parameters used appears in Table 3. After signal generation, we ran the detector and marked the times of false alarm events and misdetection events. Two examples for such simulation runs, one in low SNR and the other in high SNR, are given in Figures 11a and 11c, respectively.

Table 3:
Simulation Parameters.
ParameterValueDescription
Δt 8.33e5 sec Simulation time step 
R0 40 Hz Spontaneous firing rate of a neuron in the armed state 
R1 380 Hz Recovery rate of a neuron 
y 380 Hz Stimulus 
Rss 36.1 Hz Rss = R0R1/(R0 + R1
σ2x/Rss 0.0065 V2t Signal power due to spontaneous activity of unrelated cells 
σ2n (two values: in low and high SNR)  Signal power due to spontaneous activity of unrelated cells 
σ2v 0.001 V2t Power of additive white noise 
Mmin 1152 Minimal allowed number of simulation steps between adjacent stimuli 
Ms 256 Stimulus duration 
Md Allowed detection delay 
α 6.51eProbability of stimulus onset to occur after minimal duration is over 
CFA/CMD  Ratio of cost of false alarm to cost of misdetection; several values were used 
ParameterValueDescription
Δt 8.33e5 sec Simulation time step 
R0 40 Hz Spontaneous firing rate of a neuron in the armed state 
R1 380 Hz Recovery rate of a neuron 
y 380 Hz Stimulus 
Rss 36.1 Hz Rss = R0R1/(R0 + R1
σ2x/Rss 0.0065 V2t Signal power due to spontaneous activity of unrelated cells 
σ2n (two values: in low and high SNR)  Signal power due to spontaneous activity of unrelated cells 
σ2v 0.001 V2t Power of additive white noise 
Mmin 1152 Minimal allowed number of simulation steps between adjacent stimuli 
Ms 256 Stimulus duration 
Md Allowed detection delay 
α 6.51eProbability of stimulus onset to occur after minimal duration is over 
CFA/CMD  Ratio of cost of false alarm to cost of misdetection; several values were used 
Figure 11:

Two simulations of a multiunit response in low SNR (top), and high SNR (bottom). (a,c) In each of the left figures, the top signal is the generated multiunit signal, the detector output flag is the third signal from the top (solid black line), and the false alarm and misdetection costs appear above and below the detector output flag, respectively. (b,d) The right figures are the performance curves (〈PMD〉 as a function of 〈PFA〉) in two different SNR conditions. Note that the ad hoc detector is closer to the optimum in low-SNR conditions.

Figure 11:

Two simulations of a multiunit response in low SNR (top), and high SNR (bottom). (a,c) In each of the left figures, the top signal is the generated multiunit signal, the detector output flag is the third signal from the top (solid black line), and the false alarm and misdetection costs appear above and below the detector output flag, respectively. (b,d) The right figures are the performance curves (〈PMD〉 as a function of 〈PFA〉) in two different SNR conditions. Note that the ad hoc detector is closer to the optimum in low-SNR conditions.

Figure 12 shows the evolution of the errors probabilities PMD[k] and PFA[k] along time in two different SNRs. The probability of a misdetection error (see Figure 12a) is at its peak on stimulus onset. In high-SNR conditions, the probability of misdetection errors declines sharply shortly after stimulus onset; that is, the detection in high SNR is better and quicker. The high concentration of misdetection errors near stimulus onset also means that allowing a greater detection delay (Md) can also reduce these errors significantly. The distribution of the false alarm errors is higher just after stimulus fall, and the remaining errors are spread approximately uniformly between the end of the minimal intertrial interval (Tmin) to stimulus onset.

Figure 12:

PFA(jΔt) and PMD(jΔt) along time in high and low SNR conditions. In both SNRs, threshold was chosen so that 〈PFA〉 = 0.0035 and the number of ensemble averages was Navg = 32. The stimuli are marked by the dark gray rectangular puls, and the errors in low and high SNRs are marked in light gray and dotted black, respectively. (a) The distribution of the misdetection errors is at its peak on stimulus onset. In high SNR, the probability of a misdetection is lower than in low SNR, and it declines sharply shortly after stimulus onset; the detection is better and quicker. (b) The distribution of a false alarm error is at its peak at stimulus fall, and the remaining errors are distributed approximately uniformly between the end of the minimal silent interval (Tmin) to stimulus onset. The total false alarm probability in both SNRs is about the same because we specifically chose the thresholds to meet this condition.

Figure 12:

PFA(jΔt) and PMD(jΔt) along time in high and low SNR conditions. In both SNRs, threshold was chosen so that 〈PFA〉 = 0.0035 and the number of ensemble averages was Navg = 32. The stimuli are marked by the dark gray rectangular puls, and the errors in low and high SNRs are marked in light gray and dotted black, respectively. (a) The distribution of the misdetection errors is at its peak on stimulus onset. In high SNR, the probability of a misdetection is lower than in low SNR, and it declines sharply shortly after stimulus onset; the detection is better and quicker. (b) The distribution of a false alarm error is at its peak at stimulus fall, and the remaining errors are distributed approximately uniformly between the end of the minimal silent interval (Tmin) to stimulus onset. The total false alarm probability in both SNRs is about the same because we specifically chose the thresholds to meet this condition.

As a first-order approximation for the risk, we considered the time-averaged error from each type:
formula
7.4
where
formula
7.5
formula
7.6
and Tsim is the simulation run time.
The approximation of the error probabilities using their constant time averaged values leads to the relation
formula
7.7

Thus, the risk average over time is an approximation of the reciprocal of the expected time to injury. It is important to note that 〈PMD〉 and 〈PFA〉 depend on stimulus duration and also on the spacing between stimuli, as shown in Figure 12a. Therefore, one can fairly compare the 〈PMD〉 of two tests only if the stimuli duration and spacing are the same. Figure 11b shows the trade-off between the time-averaged misdetection error and the false alarm error, which is achieved by changing the ratio between the misdetection cost to false alarm cost. For a Bayesian detector, the ratio between the costs is the detection threshold. This graph is useful because the exact risk of a false alarm (CFA) is a subjective quantity that can vary from one case to another.

Given the nontrivial structure of the detector, it is also reasonable to inquire whether there are simpler solutions that achieve similar performance. Here we compared the detector against a digital realization of the Weber and Buchwald (1965) detector, which basically consists of a squaring stage, an integrator, and a threshold mechanism. The integration length of this detector is a free parameter that needs to be tuned. We used a moving-average digital integrator and found that an integration length of Ns/4 yielded the best performance. Figures 11b and 11d show a performance comparison between the optimal detector and the digital realization of the Weber and Buchwald (1965) detector. The optimal detector has a higher probability of true detection for every false alarm probability. The performance of the square-and-integrate detector is closer to the optimum in low-SNR conditions than in high-SNR conditions. This feature can be explained by showing that in low-SNR conditions, the optimal decision function, equation 6.19, could also be expressed as a weighted sum of squares of filtered (whitened) data.

7.2.  Experimental Results.

In addition to the simulations, we have also tested the detector on real electrophysiological data recorded from the Pontine nucleus of a rat as part of the ReNaChip project (www.renachip.org).4 The experimental procedure for obtaining the data is described in Taub and Mintz (2010), and a comprehensive analysis of more of these experiments is in preparation (Nossenson, Taub, & Messer, 2011).

Figure 13a shows an example of raw multiunit data together with the stimuli indicator and the detection output flag. The resemblance between the real data (see Figure 13a) and the simulated data (see Figure 11a and 11c) is evident. Figure 13b shows a performance comparison between the proposed detector and the square-and-integrate detector (Weber & Buchwald, 1965; Basseville & Nikiforov, 1993). The figure shows that the performance of the optimal detector is also superior when testing it on real electrophysiological multiunit data. The real data performance curve resembles the simulated data performance curve in Figure 11d, and in both figures, the performance gap in favor of the optimal detector is particularly noticeable in low false alarm probabilities.

Figure 13:

Example for the detection of auditory stimuli from real electrophysiological data taken from the Pontine nucleus of a rat. (a) The top signal is the raw data. The solid black line (second from the bottom) is the detector output flag. The actual locations of the stimuli appear at the bottom, and the misdetection cost is in light gray. (b) The performance of the optimal detector for these data (black line) and the performance of the ad hoc detector (gray with rectangular marks). Data supplied courtesy of Aryeh H. Taub.

Figure 13:

Example for the detection of auditory stimuli from real electrophysiological data taken from the Pontine nucleus of a rat. (a) The top signal is the raw data. The solid black line (second from the bottom) is the detector output flag. The actual locations of the stimuli appear at the bottom, and the misdetection cost is in light gray. (b) The performance of the optimal detector for these data (black line) and the performance of the ad hoc detector (gray with rectangular marks). Data supplied courtesy of Aryeh H. Taub.

8.  Summary, Discussion, and Future Work

8.1.  Summary.

We have discussed the problem of optimal detection of stimuli by observing multiunit recordings. More specifically, we have focused on recordings from brain regions that are densely populated by refractory Markov stimulus-responsive neurons.

We started by defining the optimality criteria to be the life span of the subject in situations where aversive stimuli are admitted in blocks that are spaced by exponentially distributed intervals. Then we showed that under the required measurement conditions, the observed signal is gaussian with a specific structure of the covariance function. We used the optimality criteria as well as the statistical nature of the observed signal to derive an optimal decision rule for declaring a stimulus presence. We tested the optimal detector by means of simulations and also using real electrophysiological data, and we showed that its performance is indeed superior. The detector proposed here is particularly useful in cases where good detection probabilities are required in very low false alarm probabilities. It is also useful as a reference for evaluating the performances of suboptimal detectors.

8.2.  Discussion: Relation to Other Detection Schemes.

The analytic detector introduced here has interesting ties to heuristic detection schemes and the analysis of neural networks. Specifically, the artificial detector introduced here continuously compares the incoming data to many delayed versions of a specific hypothesis it is programmed to test. The comparisons with the hypotheses are technically done by eliminating the expected mean response and then scaling, whitening, and squaring the resulting error. The resulting square errors from these comparisons are continuously accumulated and are propagated along a path that in our case consists of electronic registers. Simultaneously, as the square errors accumulate and propagate, the artificial detector sends the total square error from each register to another unit, which is in charge of raising the detection flag. The detection flag is raised following a calculation that involves summing exponential functions of the errors.

The heuristic detector also sums scaled square errors; however, several differences exist. First, the errors calculated by the heuristic method are not whitened before they are summed. Second, the error scaling method and the length of the accumulation window are different. Third, in the decision step of the heuristic approaches, the square errors are summed linearly and are then compared to a threshold, whereas in the detector introduced here, the exponents of the errors are summed and then thresholded.

Interestingly enough, as part of our signal analysis, we considered a neuron model that is capable of propagating positive quantities along a path, and can also exhibit an exponential response to the stimulus. This raises interesting associations that deserve exploration.

8.3.  Future Work.

The parametric construction of the detector enables a well-defined procedure for tuning it to measurement conditions. However, we only briefly described these parameter estimation procedures. The parameters to be estimated include the spike spectral shape, the background noise spectral shape, and the firing rate curve. We are currently working on an experimental paper (Nossenson et al., 2011) in which we discuss the multiunit response to auditory stimuli in the Pontine nucleus. In the same paper we also describe the parameter estimation procedures in more detail.

Appendix A:  State-Space Matrices

  1. is a matrix of size (Qx,max + Qn,max) × (Qx,max + Qn,max), which consists of two submatrices:
    formula
    formula
    See Table 4 for possible values of ak.
  2. is a time-dependent matrix of size (Qx,max + Qn,max) × 2, which consists of two column vectors:
    formula
    formula
    formula
    See Table 4 for possible values of bk.
  3. is a constant vector of size 1 × (Qx,max + Qn,max):
    formula
  4. is the following 1 × 3 vector:
    formula
Table 4:
Digital Implementation of g(t).
m0123456
am  −4.1519  −6.7732 3.6776 −1.1235  
bm −0.1329  −0.2523 −0.0530 0.1161 −0.0300 −0.0006 
m0123456
am  −4.1519  −6.7732 3.6776 −1.1235  
bm −0.1329  −0.2523 −0.0530 0.1161 −0.0300 −0.0006 

, , Qa = 6, Qb = 6.

Magnitude of the Fourier transform of a spike shape taken from Milstein and Koch (2008) and then high-pass-filtered above 300 Hz.

Appendix B:  Calculation of the Posterior Probabilities by the Detector

In this section we explain the recursive calculation of the probabilities and . Throughout the section, we use the following brief notations:
formula
B.1
formula
B.2
and
formula
B.3
The parameter α is the reciprocal of the mean number of clock cycles that are added to the minimal intertrial intervals (Tmin). The parameter Mmin is the number of clock cycles within the period Tmin. The quantity qk(t) stands for the logarithm of the probability that the current observation fits the hypothesis that a stimulus occurred kΔt seconds ago. The method for calculating qk(t) was described by Schweppe (1965) and is summarized in appendix C. In this section, assume that qk(t) is known.

We now show that the probabilities and can be expressed using their values at the previous sampling moment and using qk(t), Mmin and α:

Claim 1. The probability that an onset occurred in the most recent observation is
formula
B.4
Proof. The update rule for Pr(ttonset = Δt| rt0) originates from the fact that an onset can take place now only if the last onset occurred exactly ago or if the last onset occurred more than ago. Thus, the probability that an onset is currently observed can be written as
formula
B.5
Claim 2. The probability that an onset occurred more than Tmin ago is
formula
B.6
Proof: The update rule for originates from the fact that an onset occurred ago, or before, if no new onset is taking place now and the last onset occurred either exactly ago or before (more than ago). Thus, the probability can be written as the sum of two mutually exclusive events:
formula
B.7
Claim 3. The probability that an onset occurred kΔt seconds ago with Δt < kΔt < Tmin is
formula
B.8
Proof. The probability that an onset occurred ago where , depends on how well the observations that have been received up to now fit these hypotheses. These probabilities are calculated as follows:
formula
B.9

Remark 1. Since all the probabilities above are multiplied by the same factor, , the latter is immaterial for the decision rule.

Remark 2. By dividing the expressions in equations B.5 and B.8 by B.6, which is the expression for , and then taking the logarithm on both sides, we obtain the register update rules.

Appendix C:  Calculation of the Conditional Probabilities, qk(t)

In this section we summarize the method developed by Schweppe (1965) for calculating the logarithm of the conditional probability:
formula
According to Schweppe (1965), the formula for calculating qk(t) is
formula
C.1
where r[jΔt] is the received electrode voltage at time t = jΔt and the quantities Σr,k|k−1[jΔt] and are obtained using the Kalman algorithm:
  1. Kalman predict step. In this step we assume that the estimated state, and the covariance of the state estimation error, are known, and we use them to generate a prediction for the next state as follows:
    formula
    formula
    formula
    formula
  2. Kalman correct step. In this step, we correct our state prediction using the newly received sample, r(t). We also update the covariance matrix of the state estimation error:
    formula
    formula
    where is defined as
    formula

At the initialization of online operation (t = 0), the values of are set differently:

  1. All state vectors are set to zero:
    formula
  2. The value of is determined by solving the following equation:
    formula
    This equation results from taking the expectation of the square of equation 4.11. The other matrices, with k = 1, 2, …, Mmin−1 are calculated in a recursive manner using equation C.4 where equation C.10 is used as the initial value of the recursion.

Note that the matrices are deterministic quantities that settle to stable values after some transient period. The steady-state value of the matrix is found by solving the discrete algebraic Riccati equation:
formula
C.11
The other steady-state values of the matrices, with k = 1, 2, …, Mmin−1 are calculated in a recursive manner using equation C.4 where equation C.11 is used as the initial value of the recursion.

Acknowledgments

The research leading to these results has received partial funding from the European Community's 7th Framework Programme (FP7/2007-2011) under grant agreement no. 216809. We thank all our partners in this program and, in particular, Arie Taub for contributing the electrophysiological multiunit data. We thank Mati Mintz for his explanations on brain research, Yosi Shacham for the physical explanations, and Ytai Ben-Tsvi, Majd Zreik, Roni Hogri, Robert Prückl, Ivan Herreros, Andrea Giovannucci, and Paolo Del Giudice for the helpful discussions. We are grateful for the efforts of Angela Silmon, Mira Marcus-Kalish, Christoph Guger, and Paul F.M.J. Verschure for making the ReNaChip project happen. We thank the anonymous reviewers for their constructive comments.

Notes

1

The variance is a special case of the covariance where t1 = t2 = t.

2

We chose the work of Schweppe (1965) over later works by Duncan (1968) and Kailath (1969) since it fits the Stratonovich formulation, whereas the other two works fit the Ito formulation.

3

Note that the error probabilities PFA(t) and PMD(t) are no longer conditioned by the observations rt0. Part of the difficultly originates from the need to average out all possible observations at all possible scenarios of stimuli delivery.

4

These data were supplied courtesy of Aryeh H. Taub.

References

Aston-Jones
,
G.
, &
Bloom
,
F.
(
1981
).
Nonrepinephrine-containing locus coeruleus neurons in behaving rats exhibit pronounced responses to non-noxious environmental stimuli
.
Journal of Neuroscience
,
1
(
8
),
887
900
.
Barbieri
,
R.
,
Frank
,
L.
,
Nguyen
,
D.
,
Quirk
,
M.
,
Solo
,
V.
,
Wilson
,
M.
, et al
(
2004
).
Dynamic analyses of information encoding in neural ensembles
.
Neural Computation
,
16
(
2
),
277
307
.
Basseville
,
M.
, &
Nikiforov
,
I.
(
1993
).
Detection of abrupt changes: Theory and application
.
Upper Saddle River, NJ
:
Prentice Hall
.
Beck
,
J.
, &
Pouget
,
A.
(
2007
).
Exact inferences in a neural implementation of a hidden Markov model
.
Neural Computation
,
19
(
5
),
1344
1361
.
Brown
,
E.
,
Kass
,
R.
, &
Mitra
,
P.
(
2004
).
Multiple neural spike train data analysis: State-of-the-art and future challenges
.
Nature Neuroscience
,
7
(
5
),
456
461
.
Burkitt
,
A.
(
2006
).
A review of the integrate-and-fire neuron model: I. Homogeneous synaptic input
.
Biological Cybernetics
,
95
(
1
),
1
19
.
Caird
,
D.
, &
Klinke
,
R.
(
1983
).
Processing of binaural stimuli by cat superior olivary complex neurons
.
Experimental Brain Research
,
52
(
3
),
385
399
.
Camproux
,
A.
,
Saunier
,
F.
,
Chouvet
,
G.
,
Thalabard
,
J.
, &
Thomas
,
G.
(
1996
).
A hidden Markov model approach to neuron firing patterns
.
Biophysical Journal
,
71
(
5
),
2404
2412
.
Creutzfeldt
,
O.
,
Hellweg
,
F.
, &
Schreiner
,
C.
(
1980
).
Thalamocortical transformation of responses to complex auditory stimuli
.
Experimental Brain Research
,
39
(
1
),
87
104
.
Deneve
,
S.
(
2008
).
Bayesian spiking neurons I: Inference
.
Neural Computation
,
20
(
1
),
91
117
.
Di Lorenzo
,
P.
, &
Schwartzbaum
,
J.
(
1982
).
Coding of gustatory information in the pontine parabrachial nuclei of the rabbit: Temporal patterns of neural response
.
Brain Research
,
251
(
2
),
245
257
.
Duncan
,
T.
(
1968
).
Evaluation of likelihood functions
.
Information and Control
,
13
(
1
),
62
74
.
Gold
,
C.
,
Henze
,
D.
,
Koch
,
C.
, &
Buzsaki
,
G.
(
2006
).
On the origin of the extracellular action potential waveform: A modeling study
.
Journal of Neurophysiology
,
95
(
5
),
3113
3128
.
Heetderks
,
W.
(
1978
).
Criteria for evaluating multiunit spike separation techniques
.
Biological Cybernetics
,
29
(
4
),
215
220
.
Herz
,
A.
,
Gollisch
,
T.
,
Machens
,
C.
, &
Jaeger
,
D.
(
2006
).
Modeling single-neuron dynamics and computations: A balance of detail and abstraction
.
Science
,
314
(
5796
),
80
85
.
Hill
,
D.
,
Mehta
,
S.
, &
Kleinfeld
,
D.
(
2011
).
Quality metrics to accompany spike sorting of extracellular signals
.
Journal of Neuroscience
,
31
(
24
),
8699
8705
.
Holstein
,
S.
,
Buchwald
,
J.
, &
Schwafel
,
J.
(
1969a
).
Progressive changes in auditory response patterns to repeated tone during normal wakefulness and paralysis
.
Brain Research
,
16
(
1
),
133
148
.
Holstein
,
S.
,
Buchwald
,
J.
, &
Schwafel
,
J.
(
1969b
).
Tone response patterns of the auditory nuclei during normal wakefulness, paralysis, and anesthesia
.
Brain Research
,
15
(
2
),
483
499
.
Ingham
,
N.
, &
McAlpine
,
D.
(
2005
).
GABAergic inhibition controls neural gain in inferior colliculus neurons sensitive to interaural time differences
.
Journal of Neuroscience
,
25
(
26
),
632
645
.
Johnson
,
D.
(
1996
).
Point process models of single-neuron discharges
.
Journal of Computational Neuroscience
,
3
(
4
),
275
299
.
Kailath
,
T.
(
1969
).
A general likelihood-ratio formula for random signals in gaussian noise
.
IEEE Transactions on Information Theory
,
15
(
3
),
350
361
.
Kalman
,
R.
(
1960
).
A new approach to linear filtering and prediction problems. Trans
.
ASME, J. Basic Eng.
,
82D
,
35
45
.
Legatt
,
A.
,
Arezzo
,
J.
, &
Vaughan
Jr.,
H.
(
1980
).
Averaged multiple unit activity as an estimate of phasic changes in local neuronal activity: Effects of volume-conducted potentials
.
Journal of Neuroscience Methods
,
2
(
2
),
203
217
.
Lewicki
,
M.
(
1998
).
A review of methods for spike sorting: The detection and classification of neural action potentials
.
Network: Computation in Neural Systems
,
9
(
4
),
53
78
.
Logothetis
,
N.
,
Kayser
,
C.
, &
Oeltermann
,
A.
(
2007
).
In vivo measurement of cortical impedance spectrum in monkeys: Implications for signal propagation
.
Neuron
,
55
(
5
),
809
823
.
Martinez
,
J.
,
Pedreira
,
C.
,
Ison
,
M.
, &
Quian Quiroga
,
R.
(
2009
).
Realistic simulation of extracellular recordings
.
Journal of Neuroscience Methods
,
184
(
2
),
285
293
.
Milstein
,
J.
, &
Koch
,
C.
(
2008
).
Dynamic moment analysis of the extracellular electric field of a biologically realistic spiking neuron
.
Neural Computation
,
20
(
8
),
2070
2084
.
Nossenson
,
N.
, &
Messer
,
H.
(
2010
).
Modeling neuron firing pattern using a two state Markov chain
. In
Sensor Array and Multichannel Signal Processing Workshop (SAM), 2010 IEEE
(pp.
41
44
).
Piscataway, NJ
:
IEEE
.
Nossenson
,
N.
,
Taub
,
A. H.
, &
Messer
,
H.
(
2010
).
The pontine nucleus response to auditory stimuli and its detection
.
Unpublished manuscript
.
Novak
,
P.
,
Klemp
,
J.
,
Ridings
,
L.
,
Lyons
,
K.
,
Pahwa
,
R.
, &
Nazzaro
,
J.
(
2009
).
Effect of deep brain stimulation of the subthalamic nucleus upon the contralateral subthalamic nucleus in Parkinson disease
.
Neuroscience Letters
,
463
(
1
),
12
16
.
Paninski
,
L.
,
Pillow
,
J.
, &
Lewi
,
J.
(
2007
).
Statistical models for neural encoding, decoding, and optimal stimulus design
.
Progress in Brain Research
,
165
,
493
507
.
Radner
,
W.
,
Sadda
,
S. R.
,
Humayun
,
M. S.
,
Suzuki
,
S.
,
Melia
,
M.
,
Weiland
,
J.
, et al
(
2001
).
Light-driven retinal ganglion cell responses in blind RD mice after neural retinal transplantation
.
Investigative Ophthalmology and Visual Science
,
42
(
5
),
1057
1065
.
Reches
,
A.
, &
Gutfreund
,
Y.
(
2008
).
Stimulus-specific adaptations in the gaze control system of the barn owl
.
Journal of Neuroscience
,
28
(
6
),
1523
1533
.
Ricciardi
,
L.
, &
Esposito
,
F.
(
1966
).
On some distribution functions for non-linear switching elements with finite dead time
.
Biological Cybernetics
,
3
(
3
),
148
152
.
Roweis
,
S.
, &
Ghahramani
,
Z.
(
1999
).
A unifying review of linear gaussian models
.
Neural Computation
,
11
(
2
),
305
345
.
Schwartz
,
E.
,
Ramos
,
A.
, &
John
,
E.
(
1976
).
Single cell activity in chronic unit recording: A quantitative study of the unit amplitude spectrum
.
Brain Research Bulletin
,
1
(
1
),
57
68
.
Schweppe
,
F.
(
1965
).
Evaluation of likelihood functions for gaussian signals
.
IEEE Transactions on Information Theory
,
11
(
1
),
61
70
.
Stratonovich
,
R.
(
1966
).
A new representation for stochastic integrals and equations
.
Siam J. Control
,
4
,
362
371
.
Taberner
,
A.
, &
Liberman
,
M.
(
2005
).
Response properties of single auditory nerve fibers in the mouse
.
Journal of Neurophysiology
,
93
(
1
),
557
569
.
Taub
,
A.
, &
Mintz
,
M.
(
2010
).
Amygdala conditioning modulates sensory input to the cerebellum
.
Neurobiology of Learning and Memory
,
94
,
521
529
.
Teich
,
M.
,
Matin
,
L.
, &
Cantor
,
B.
(
1978
).
Refractoriness in the maintained discharge of the cat's retinal ganglion cell
.
J. Opt. Soc. Am
,
68
(
3
),
386
402
.
Tommerdahl
,
M.
,
Delemos
,
K.
,
Whitsel
,
B.
,
Favorov
,
O.
, &
Metz
,
C.
(
1999
).
Response of anterior parietal cortex to cutaneous flutter versus vibration
.
Journal of Neurophysiology
,
82
(
1
),
16–33
.
Toyoizumi
,
T.
,
Rad
,
K.
, &
Paninski
,
L.
(
2009
).
Mean-field approximations for coupled populations of generalized linear model spiking neurons with Markov refractoriness
.
Neural Computation
,
21
(
5
),
1203
1243
.
Tsumoto
,
T.
,
Creutzfeldt
,
O.
, &
Legendy
,
C.
(
1978
).
Functional organization of the corticofugal system from visual cortex to lateral geniculate nucleus in the cat
.
Experimental Brain Research
,
32
(
3
),
345
364
.
Wang
,
X.
,
Lu
,
T.
,
Snider
,
R.
, &
Liang
,
L.
(
2005
).
Sustained firing in auditory cortex evoked by preferred stimuli
.
Nature
,
435
(
7040
),
341
346
.
Watkins
,
P.
,
Santhanam
,
G.
,
Shenoy
,
K.
, &
Harrison
,
R.
(
2004
).
Validation of adaptive threshold spike detector for neural recording
. In
Proceedings of Engineering in Medicine and Biology Society, 26th Annual International Conference of the IEEE
(Vol. 2, pp.
4079
4082
).
Piscataway, NJ
:
IEEE
.
Weber
,
A.
,
Chen
,
H.
,
Hubbard
,
W.
, &
Kaufman
,
P.
(
2000
).
Experimental glaucoma and cell size, density, and number in the primate lateral geniculate nucleus
.
Investigative Ophthalmology and Visual Science
,
41
(
6
),
1370
1379
.
Weber
,
D.
, &
Buchwald
,
J.
(
1965
).
A technique for recording and integrating multiple unit activity simultaneously with the EEG in chronic animals
.
Electroencephalography and Clinical Neurophysiology
,
19
(
2
),
190
192
.