## 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.

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, *T _{d}*, whether a warning signal is present. If the detection is too late, an aversive stimulus hurts the subject.

*l*th intertest-interval interval is designated as

*T*(

_{ITI}*l*) and obeys the following probability distribution: where

*T*

_{min}is the shortest possible interval between adjacent stimuli and

*T*

_{avg}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 *T _{s}* 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

*T*equal to zero.

_{d}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.

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 *r ^{t}*

_{0}(the electrode voltage from start time to current time).

*t*,

*t*+ Δ

_{t}], is proportional to the momentary risk times the decision spacing interval: It follows that the probability of staying uninjured at time

*t*=

*N*· Δ

_{t}(i.e.,

*N*time intervals from time axis origin) is Thus, at the limit Δ

_{t}→ 0, the probability of staying uninjured can be expressed as The goal of the detector is to decide whether to risk misdetection

*I*(

_{off}*t*) = 1 or a false alarm

*I*(

_{on}*t*) = 1 so that the total risk is minimized and the expected life span of the subject is maximized: where,

*E*{

*T*} is the expectancy of the first injury time:

_{uninjured}The considerations of the detector are illustrated in Figure 3.

We now clarify a few subtle points:

- •
ℜ

_{injury}(*t*),*C*(_{MD}*t*), and*C*are given in units of 1/sec or Hz._{FA} - •
*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. - •
*I*(_{on}*t*) and its complement are indicator functions for the decision of the detector.*I*(_{on}*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.

Mark . | Meaning . | Unit . |
---|---|---|

T _{ITI} | Intertrial interval. A random variable that represents interval duration between consecutive onsets of aversive stimuli. | sec |

T_{min} | Shortest possible interval between consecutive onsets of the aversive stimuli. | sec |

T _{s} | Duration of the neural response to the warning signal. | sec |

T _{d} | Delay of the onset of the aversive stimulus from the onset of the warning signal. | sec |

C _{MD} | The momentary cost for misdetecting the aversive stimulus while it is on. | Hz |

C _{FA} | The momentary cost for falsely alarming on the presence of an aversive stimulus. | Hz |

r^{t}_{0} | 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 |

Mark . | Meaning . | Unit . |
---|---|---|

T _{ITI} | Intertrial interval. A random variable that represents interval duration between consecutive onsets of aversive stimuli. | sec |

T_{min} | Shortest possible interval between consecutive onsets of the aversive stimuli. | sec |

T _{s} | Duration of the neural response to the warning signal. | sec |

T _{d} | Delay of the onset of the aversive stimulus from the onset of the warning signal. | sec |

C _{MD} | The momentary cost for misdetecting the aversive stimulus while it is on. | Hz |

C _{FA} | The momentary cost for falsely alarming on the presence of an aversive stimulus. | Hz |

r^{t}_{0} | 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 *r ^{t}*

_{0}. 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

*r*

^{t}_{0}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/mm

^{3}(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.

*r*(

*t*), results from the cumulative action of three processes: 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.

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, Σ(*t*_{1}, *t*_{2}). 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.

*t*,

*t*+ Δ

_{t}) is where

*y*(

*t*) is proportional to the short-term average intensity of the stimulus and the constant

*R*

_{0}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),

*R*

_{1}.

*P*

_{1}= 1 −

*P*

_{0}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 where

*t*

_{0}is the start time and

*P*

_{0}(

*t*

_{0}) is the probability of being in an armed state at the start time.

Figure 5 depicts the shape of the function *R _{fire}*(

*t*) when the stimulus

*y*(

*t*) is a rectangular pulse.

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.

_{x}(

*t*

_{1},

*t*

_{2}), will always have a shape that resembles the firing rate regardless of the electrode location, where:

- •
*g*(_{x}*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 ∫ψ(
*t*_{1},*t*_{2}, τ)○*d*τ is defined as 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.

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 *t*_{1} = *t*_{2}.

#### 3.2.2. Covariance of the Background Signal.

*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:

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

#### 3.2.3. Additive White Noise.

^{2}

_{v}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.

r(t) = x(t) + n(t) + v(t) |

x(t) ∼ N(μ_{x}(t), Σ_{x}(t_{1,}t_{2})) |

n(t) ∼ N(0, Σ_{n}(t_{1,}t_{2})) |

v(t) ∼ N(0, Σ_{v}(t_{1,}t_{2})) |

Σ_{v}(t_{1,}t_{2}) = σ^{2}_{v}○δ(t_{2} − t_{1}) |

The mean μ_{x}(t) depends on brain-electrode geometry and has to be measured empirically. |

A coarse estimate of the firing rate R(τ) 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 _{fire}R_{0}, R_{1}, and then calculating R(τ) using equations 3.9 and 3.10. _{fire} |

The parameters σ^{2}_{n}, σ^{2}_{v} are found empirically. |

g(_{x}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. |

g(_{n}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) ∼ N(μ_{x}(t), Σ_{x}(t_{1,}t_{2})) |

n(t) ∼ N(0, Σ_{n}(t_{1,}t_{2})) |

v(t) ∼ N(0, Σ_{v}(t_{1,}t_{2})) |

Σ_{v}(t_{1,}t_{2}) = σ^{2}_{v}○δ(t_{2} − t_{1}) |

The mean μ_{x}(t) depends on brain-electrode geometry and has to be measured empirically. |

A coarse estimate of the firing rate R(τ) 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 _{fire}R_{0}, R_{1}, and then calculating R(τ) using equations 3.9 and 3.10. _{fire} |

The parameters σ^{2}_{n}, σ^{2}_{v} are found empirically. |

g(_{x}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. |

g(_{n}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.

*x*(

*t*) is

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).

*g*(

_{x}*t*), using a recursive infinite impulse response (IIR) filter: Figure 8 compares the impulse response of an IIR filter with

*Q*

_{x,a}= 6,

*Q*

_{x,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.

Recall that equations 3.11 and 3.12 imply that the phase response of *G _{x}*(

*s*) is immaterial for the expression of the covariance function, so there exists some degree of freedom for choosing the coefficients

*b*

_{x,m}.

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

_{n}(

*t*

_{1},

*t*

_{2}) associated with the background spontaneous spikes is obtained if

*n*(

*t*) obeys the following stochastic difference equation: 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): where

*u*[

_{r}*j*Δ

_{t}] is the process mean and

*dw*

_{[j−l]},

*d*ζ

_{[j−l]}, and

*d*ξ

_{j}are independent standard white noise processes.

is a vector of internal state variables of size (

*Q*_{x,max}+*Q*_{n,max}) × 1. The first entry of this vector, , corresponds to*x*(*k*Δ_{t}) minus the latest innovation of*x*(*k*Δ_{t}). The following (*Q*_{x,a}− 1) entries are used for the construction of . Entry number*Q*_{x,a}+ 1 of this vector, , corresponds to*n*(*k*Δ_{t}) minus the latest innovation of*n*(*k*Δ_{t}). The following (*Q*_{n,a}− 1) entries are used for the construction of .The matrices and are given in appendix A. They generate and .

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.

_{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: where

*P*(

_{FA}*t*∈

*T*|

_{FA}**r**

^{t}

_{0}) is the probability that the current time falls inside the false alarm zone and

*P*(

_{MD}*t*∈

*T*|

_{MD}**r**

^{t}

_{0}) is the probability that the current time is within the detection window.

_{injury}(

*t*), the following decision rule should be applied: 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

*P*(

_{FA}*t*∈

*T*|

_{FA}**r**

^{t}

_{0}) and

*P*(

_{MD}*t*∈

*T*|

_{MD}**r**

^{t}

_{0}) 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

*T*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

_{d}*T*and smaller than

_{d}*T*(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

_{s}*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 *T _{d}* and smaller than

*T*is larger than the sum of probabilities of the remaining alternatives.

_{s}## 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.

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

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

### 6.1. The Register Array.

*M*

_{min}normalized posterior probabilities. Each of the posterior probabilities has a corresponding register designated as

*d*[

_{k}*j*Δ

_{t}]. The first register,

*d*

_{1}[

*j*Δ

_{t}], corresponds to the posterior probability: The second register,

*d*

_{2}[

*n*Δ

_{t}], corresponds to the posterior probability, and so on. The last register, , corresponds to the posterior probability: 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: As a result of this normalization, register #

*M*

_{min}is constantly kept at zero:

_{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: and The parameter α is the reciprocal of the mean number of clock cycles added to the minimal intertrial intervals (

*T*

_{min}), and it is a known a priori. The value of

*q*[

_{k}*n*Δ

_{t}] is calculated by the Kalman estimators array (see section 6.2) using the method in Schweppe (1965). Here, assume that

*q*[

_{k}*n*Δ

_{t}] is known.

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

**Proof.** See appendix B.

### 6.2. The Kalman Estimators Array.

*M*

_{min}units. Each unit outputs the quantity

*q*(

_{k}*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, 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: . - •
- •
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.

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 *k*th 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 *M*_{min} estimators receives three new inputs and produces three new outputs. The inputs of the *k*th Kalman filter at time *t* = *j*Δ_{t} are the mean free version of the current electrode voltage, *r*[*j*Δ_{t}] − *u _{x}*[

*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

*q*[

_{k}*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 (#*M*_{min}) are an exception. The first estimator receives the inputs from Kalman estimator #*M*_{min}. The last Kalman estimator (#*M*_{min}) produces two different sets of predictions:

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

. This prediction is sent back to the Kalman estimator #

*M*_{min}. It corresponds to the hypothesis that no stimulus will occur in the next sample.

*t*= 0), the values of are set differently: The procedure for calculating is given in appendix C.

### 6.3. Assertion of the Detection Flag.

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

## 7. Results

### 7.1. Simulation Study.

*P*[

_{FA}*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:

Since analytic evaluation of *P _{FA}*(

*t*) and

*P*(

_{MD}*t*) is a formidable task, we evaluated

*P*(

_{MD}*t*) and

*P*(

_{FA}*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.

Parameter . | Value . | Description . |
---|---|---|

Δ_{t} | 8.33e5 sec | Simulation time step |

R_{0} | 40 Hz | Spontaneous firing rate of a neuron in the armed state |

R_{1} | 380 Hz | Recovery rate of a neuron |

y | 380 Hz | Stimulus |

R _{ss} | 36.1 Hz | R = _{ss}R_{0}R_{1}/(R_{0} + R_{1}) |

σ^{2}_{x}/R _{ss} | 0.0065 V^{2}/Δ_{t} | Signal power due to spontaneous activity of unrelated cells |

σ^{2}_{n} (two values: in low and high SNR) | Signal power due to spontaneous activity of unrelated cells | |

σ^{2}_{v} | 0.001 V^{2}/Δ_{t} | Power of additive white noise |

M_{min} | 1152 | Minimal allowed number of simulation steps between adjacent stimuli |

M _{s} | 256 | Stimulus duration |

M _{d} | 1 | Allowed detection delay |

α | 6.51e4 | Probability of stimulus onset to occur after minimal duration is over |

C/_{FA}C _{MD} | Ratio of cost of false alarm to cost of misdetection; several values were used |

Parameter . | Value . | Description . |
---|---|---|

Δ_{t} | 8.33e5 sec | Simulation time step |

R_{0} | 40 Hz | Spontaneous firing rate of a neuron in the armed state |

R_{1} | 380 Hz | Recovery rate of a neuron |

y | 380 Hz | Stimulus |

R _{ss} | 36.1 Hz | R = _{ss}R_{0}R_{1}/(R_{0} + R_{1}) |

σ^{2}_{x}/R _{ss} | 0.0065 V^{2}/Δ_{t} | Signal power due to spontaneous activity of unrelated cells |

σ^{2}_{n} (two values: in low and high SNR) | Signal power due to spontaneous activity of unrelated cells | |

σ^{2}_{v} | 0.001 V^{2}/Δ_{t} | Power of additive white noise |

M_{min} | 1152 | Minimal allowed number of simulation steps between adjacent stimuli |

M _{s} | 256 | Stimulus duration |

M _{d} | 1 | Allowed detection delay |

α | 6.51e4 | Probability of stimulus onset to occur after minimal duration is over |

C/_{FA}C _{MD} | Ratio of cost of false alarm to cost of misdetection; several values were used |

Figure 12 shows the evolution of the errors probabilities *P _{MD}*[

*k*] and

*P*[

_{FA}*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 (

*M*) 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 (

_{d}*T*

_{min}) to stimulus onset.

Thus, the risk average over time is an approximation of the reciprocal of the expected time to injury. It is important to note that 〈*P _{MD}*〉 and 〈

*P*〉 depend on stimulus duration and also on the spacing between stimuli, as shown in Figure 12a. Therefore, one can fairly compare the 〈

_{FA}*P*〉 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 (

_{MD}*C*) is a subjective quantity that can vary from one case to another.

_{FA}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 *N _{s}*/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.

## 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

m . | 0 . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|---|---|

a _{m} | −4.1519 | −6.7732 | 3.6776 | −1.1235 | |||

b _{m} | −0.1329 | −0.2523 | −0.0530 | 0.1161 | −0.0300 | −0.0006 |

m . | 0 . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . |
---|---|---|---|---|---|---|---|

a _{m} | −4.1519 | −6.7732 | 3.6776 | −1.1235 | |||

b _{m} | −0.1329 | −0.2523 | −0.0530 | 0.1161 | −0.0300 | −0.0006 |

, , *Q _{a}* = 6,

*Q*= 6.

_{b}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

*T*

_{min}). The parameter

*M*

_{min}is the number of clock cycles within the period

*T*

_{min}. The quantity

*q*(

_{k}*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

*q*(

_{k}*t*) was described by Schweppe (1965) and is summarized in appendix C. In this section, assume that

*q*(

_{k}*t*) is known.

We now show that the probabilities and can be expressed using their values at the previous sampling moment and using *q _{k}*(

*t*),

*M*

_{min}and α:

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

## Appendix C: Calculation of the Conditional Probabilities, *q*_{k}(t)

*q*

_{k}(t)*q*(

_{k}*t*) is 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:

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

- The value of is determined by solving the following equation: This equation results from taking the expectation of the square of equation 4.11. The other matrices, with
*k*= 1, 2, …,*M*_{min−1}are calculated in a recursive manner using equation C.4 where equation C.10 is used as the initial value of the recursion.

*k*= 1, 2, …,

*M*

_{min−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 *t*_{1} = *t*_{2} = *t*.

^{3}

Note that the error probabilities *P _{FA}*(

*t*) and

*P*(

_{MD}*t*) are no longer conditioned by the observations

*r*

^{t}_{0}. 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.