Measuring the signal-to-noise ratio of a neuron

Significance Neurons represent both signal and noise in binary electrical discharges termed action potentials. Hence, the standard signal-to-noise ratio (SNR) definition of signal amplitude squared and divided by the noise variance does not apply. We show that the SNR estimates a ratio of expected prediction errors. Using point process generalized linear models, we extend the standard definition to one appropriate for single neurons. In analyses of four neural systems, we show that single neuron SNRs range from −29 dB to −3 dB and that spiking history is often a more informative predictor of spiking propensity than the signal or stimulus activating the neuron. By generalizing the standard SNR metric, we make explicit the well-known fact that individual neurons are highly noisy information transmitters.

The signal-to-noise ratio (SNR), a commonly used measure of fidelity in physical systems, is defined as the ratio of the squared amplitude or variance of a signal relative to the variance of the noise. This definition is not appropriate for neural systems in which spiking activity is more accurately represented as point processes. We show that the SNR estimates a ratio of expected prediction errors and extend the standard definition to one appropriate for single neurons by representing neural spiking activity using point process generalized linear models (PP-GLM). We estimate the prediction errors using the residual deviances from the PP-GLM fits. Because the deviance is an approximate χ 2 random variable, we compute a bias-corrected SNR estimate appropriate for single-neuron analysis and use the bootstrap to assess its uncertainty. In the analyses of four systems neuroscience experiments, we show that the SNRs are −10 dB to −3 dB for guinea pig auditory cortex neurons, −18 dB to −7 dB for rat thalamic neurons, −28 dB to −14 dB for monkey hippocampal neurons, and −29 dB to −20 dB for human subthalamic neurons. The new SNR definition makes explicit in the measure commonly used for physical systems the often-quoted observation that single neurons have low SNRs. The neuron's spiking history is frequently a more informative covariate for predicting spiking propensity than the applied stimulus. Our new SNR definition extends to any GLM system in which the factors modulating the response can be expressed as separate components of a likelihood function. SNR | signal-to-noise ratio | neuron | simulation | point processes T he signal-to-noise ratio (SNR), defined as the amplitude squared of a signal or the signal variance divided by the variance of the system noise, is a widely applied measure for quantifying system fidelity and for comparing performance among different systems (1)(2)(3)(4). Commonly expressed in decibels as 10log 10 (SNR), the higher the SNR, the stronger the signal or information in the signal relative to the noise or distortion. Use of the SNR is most appropriate for systems defined as deterministic or stochastic signals plus Gaussian noise (2,4). For the latter, the SNR can be computed in the time or frequency domain.
Use of the SNR to characterize the fidelity of neural systems is appealing because information transmission by neurons is a noisy stochastic process. However, the standard concept of SNR cannot be applied in neuronal analyses because neurons transmit both signal and noise primarily in their action potentials, which are binary electrical discharges also known as spikes (5)(6)(7)(8). Defining what is the signal and what is the noise in neural spiking activity is a challenge because the putative signals or stimuli for neurons differ appreciably among brain regions and experiments. For example, neurons in the visual cortex and in the auditory cortex respond respectively to features of light (9) and sound stimuli (10) while neurons in the somatosensory thalamus respond to tactile stimuli (11). In contrast, neurons in the rodent hippocampus respond robustly to the animal's position in its environment (11,12), whereas monkey hippocampal neurons respond to the process of task learning (13). As part of responding to a putative stimulus, a neuron's spiking activity is also modulated by biophysical factors such as its absolute and relative refractory periods, its bursting propensity, and local network and rhythm dynamics (14,15). Hence, the definition of SNR must account for the extent to which a neuron's spiking responses are due to the applied stimulus or signal and to these intrinsic biophysical properties.
Formulations of the SNR for neural systems have been studied. Rieke et al. (16) adapted information theory measures to define Gaussian upper bounds on the SNR for individual neurons. Coefficients of variation and Fano factors based on spike counts (17)(18)(19) have been used as measures of SNR. Similarly, Gaussian approximations have been used to derive upper bounds on neuronal SNR (16). These approaches do not consider the point process nature of neural spiking activity. Moreover, these measures and the Gaussian approximations are less accurate for neurons with low spike rates or when information is contained in precise spike times.
Lyamzin et al. (20) developed an SNR measure for neural systems using time-dependent Bernoulli processes to model the neural spiking activity. Their SNR estimates, based on variance formulae, do not consider the biophysical properties of the neuron and are more appropriate for Gaussian systems (16,21,22). The Poisson regression model used widely in statistics to relate count observations to covariates provides a framework for studying the SNR for non-Gaussian systems because it provides an analog of the square of the multiple correlation coefficient (R 2 ) used to measure goodness of fit in linear regression analyses Significance Neurons represent both signal and noise in binary electrical discharges termed action potentials. Hence, the standard signal-to-noise ratio (SNR) definition of signal amplitude squared and divided by the noise variance does not apply. We show that the SNR estimates a ratio of expected prediction errors. Using point process generalized linear models, we extend the standard definition to one appropriate for single neurons. In analyses of four neural systems, we show that single neuron SNRs range from −29 dB to −3 dB and that spiking history is often a more informative predictor of spiking propensity than the signal or stimulus activating the neuron. By generalizing the standard SNR metric, we make explicit the well-known fact that individual neurons are highly noisy information transmitters. (23). The SNR can be expressed in terms of the R 2 for linear and Poisson regression models. However, this relationship has not been exploited to construct an SNR estimate for neural systems or point process models. Finally, the SNR is a commonly computed statistic in science and engineering. Extending this concept to non-Gaussian systems would be greatly aided by a precise statement of the theoretical quantity that this statistic estimates (24,25).
We show that the SNR estimates a ratio of expected prediction errors (EPEs). Using point process generalized linear models (PP-GLM), we extend the standard definition to one appropriate for single neurons recorded in stimulus−response experiments. In analyses of four neural systems, we show that single-neuron SNRs range from −29 dB to −3 dB and that spiking history is often a more informative predictor of spiking propensity than the signal being represented. Our new SNR definition generalizes to any problem in which the modulatory components of a system's output can be expressed as separate components of a GLM.

Theory
A standard way to define the SNR is as the ratio where σ 2 signal is structure in the data induced by the signal and σ 2 noise is the variability due to the noise. To adapt this definition to the analysis of neural spike train recordings from a single neuron, we have: to (i) define precisely what the SNR estimates; (ii) extend the definition and its estimate to account for covariates that, along with the applied stimulus or signal input, also affect the neural response; and (iii) extend the SNR definition and its estimate so that it applies to point process models of neural spiking activity.
By analyzing the linear Gaussian signal plus noise model (Supporting Information), we show that standard SNR computations (Eq. S5) provide an estimator of a ratio of EPEs (Eq. S4). For the linear Gaussian model with covariates, this ratio of EPEs is also well defined (Eq. S6) and can be estimated as a ratio of sum of squares of residuals (Eq. S7). The SNR definition further extends to the GLM with covariates (Eq. S8). To estimate the SNR for the GLM, we replace the sums of squares by the residual deviances, their extensions in the GLM framework Eqs. S9 and S10. The residual deviance is a constant multiple of the Kullback−Leibler (KL) divergence between the data and the model. Due to the Pythagorean property of the KL divergence of GLM models with canonical link functions (26-28) evaluated at the maximum likelihood estimates, the SNR estimator can be conveniently interpreted as the ratio of the explained KL divergence of the signal relative to the noise. We propose an approximate bias correction for the GLM SNR estimate with covariates (Eq. S11), which gives the estimator better performance in low signal-to-noise problems such as single-neuron recordings. The GLM framework formulated with point process models has been used to analyze neural spiking activity (5)(6)(7)29). Therefore, we derive a point process GLM (PP-GLM) SNR estimate for single-neuron spiking activity recorded in stimulus− response experiments.
A Volterra Series Expansion of the Conditional Intensity Function of a Spiking Neuron. Volterra series are widely used to model biological systems (30), including neural spiking activity (16). We develop a Volterra series expansion of the log of the conditional intensity function to define the PP-GLM for single-neuron spiking activity (31). We then apply the GLM framework outlined in Supporting Information to derive the SNR estimate.
We assume that on an observation interval ð0, T, we record spikes at times 0 < u 1 < u 2 < ..... < u J < T. If we model the spike events as a point process, then the conditional intensity function of the spike train is defined by (5) where NðtÞ is the number of spikes in the interval ð0, t for t ∈ ð0, T and H t is the relevant history at t. It follows that for Δ small, We assume that the neuron receives a stimulus or signal input and that its spiking activity depends on this input and its biophysical properties. The biophysical properties may include absolute and relative refractory periods, bursting propensity, and network dynamics. We assume that we can express log λðtjH t Þ in a Volterra series expansion as a function of the signal and the biophysical properties (31). The first-order and second-order terms in the expansion are where sðtÞ is the signal at time t, dNðtÞ is the increment in the counting process, β S ðuÞ is the one-dimensional signal kernel, β H ðtÞ is the one-dimensional temporal or spike history kernel, h 1 ðu, vÞ is the 2D signal kernel, h 2 ðu, vÞ is the 2D temporal kernel, and h 3 ðu, vÞ is the 2D signal−temporal kernel. Eq. 4 shows that up to first order, the stimulus effect on the spiking activity and the effect of the biophysical properties of the neuron, defined in terms of the neuron's spiking history, can be expressed as separate components of the conditional intensity function. Assuming that the second-order effects are not strong, then the approximate separation of these two components makes it possible to define the SNR for the signal, also taking account of the effect of the biophysical properties as an additional covariate and vice versa. We expand the log of the conditional intensity function in the Volterra series instead of the conditional intensity function itself in the Volterra series to ensure that the conditional intensity function is positive. In addition, using the log of the conditional intensity function simplifies the GLM formulation by using the canonical link function for the local Poisson model.
Likelihood Analysis Using a PP-GLM. We define the likelihood model for the spike train using the PP-GLM framework (5). We assume the stimulus−response experiment consists of R independent trials, which we index as r = 1, ..., R. We discretize time within a trial by choosing L large and defining the L subintervals Δ = T −1 L. We choose L large so that each subinterval contains at most one spike. We index the subintervals ℓ = 1, ...L and define n r,ℓ to be 1 if, on trial r, there is a spike in the subinterval ((ℓ−1)Δ,ℓΔ) and it is 0 otherwise. We let n r = ðn r,1 , ...n r,L Þ be the set of spikes recorded on trial r in ð0, T. Let H r,ℓ = fn r,ℓ−J , ..., n r,ℓ−1 g be the relevant history of the spiking activity at time ℓΔ. We define the discrete form of the Volterra expansion by using the first two terms of Eq. 4 to obtain β H,j n r,ℓ−j , [5] where β = ðβ 0 , β S , β H Þ′, β S = ðβ S,0 ,..., β S,K Þ′, and β H = ðβ H,1 ,..., β H,J Þ′, and hence the dependence on the stimulus goes back a period of KΔ, whereas the dependence on spiking history goes back a period of JΔ. Exponentiating both sides of Eq. 5 yields .

[6]
The first and third terms on the right side of Eq. 6 measure the intrinsic spiking propensity of the neuron, whereas the second term measures the effect of the stimulus or signal on the neuron's spiking propensity. The likelihood function for β given the recorded spike train is (5) Lðn, βÞ = exp .

[7]
Likelihood formulations with between-trial dependence (32) are also possible but are not considered here.
The maximum likelihood estimate of β can be computed by maximizing Eq. 7 or, equivalently, by minimizing the residual deviance defined as where n = ðn 1 , ..., n R Þ and Lðn, nÞ is the saturated model or the highest possible value of the maximized log likelihood (26). Maximizing log Lðn, βÞ to compute the maximum likelihood estimate of β is equivalent to minimizing the deviance, because Lðn, nÞ is a constant. The deviance is the generalization to the GLM of the sum of squares from the linear Gaussian model (33). As in the standard GLM framework, these computations are carried out efficiently using iteratively reweighted least squares. In our PP-GLM likelihood analyses, we use Akaike's Information Criterion (AIC) to help choose the order of the discrete kernels β H and β S (34). We use the time-rescaling theorem and Kolmogorov−Smirnov (KS) plots (35) along with analyses of the Gaussian transformed interspike intervals to assess model goodness of fit (36). We perform the AIC and time-rescaling goodness-of-fit analyses using cross-validation to fit the model to half of the trials in the experiments (training data set) and then evaluating AIC, the KS plots on the second half the trials (test data set). The model selection and goodness-of-fit assessments are crucial parts of the SNR analyses. They allow us to evaluate whether our key assumption is valid, that is, that the conditional intensity function can be represented as a finite-order Volterra series whose second-order terms can be neglected. Significant lack of fit could suggest that this assumption did not hold and would thereby weaken, if not invalidate, any subsequent inferences and analyses.
SNR Estimates for a Single Neuron. Applying Eq. S11, we have that for a single neuron, the SNR estimate for the signal given the spike history (biophysical properties) with the approximate bias corrections iŝ and that for a single neuron, the SNR estimates of the spiking propensity given the signal isŜ where dimðβÞ is the dimension or the number of parameters inβ. Application of the stimulus activates the biophysical properties of the neuron. Therefore, to measure the effect of the stimulus, we fit the GLM with and without the stimulus and use the difference between the deviances to estimate theŜNR S (Eq. 9). Similarly, to measure the effect of the spiking history, we fit the GLM with and without the spike history and use the difference between the deviances to estimate theŜNR H (Eq. 10).

STATISTICS
In experiment 1 (Fig. 1A), neural spike trains were recorded from 12 neurons in the primary auditory cortex of anesthetized guinea pigs in response to a 200 μs/phase biphasic electrical pulse at 44.7-μA applied in the inferior colliculus (10). Note that the neural recordings were generally multi-unit responses recorded on 12 sites but we refer to them as neurons in this paper. The stimulus was applied at time 0, and spiking activity was recorded from 10 ms before the stimulus to 50 ms after the stimulus during 40 trials. In experiment 2, neural spiking activity was recorded in 12 neurons from the ventral posteromedial (VPm) nucleus of the thalamus (VPm thalamus) in rats in response to whisker stimulation (Fig. 1B) (11). The stimulus was deflection of the whisker at a velocity of 50 mm/s at a repetition rate of eight deflections per second. Each deflection was 1 mm in amplitude and began from the whiskers' neutral position as the trough of a single sine wave and ended smoothly at the same neutral position. Neural spiking activity was recorded for 3,000 ms across 51 trials.
In experiment 3 (Fig. 1C), neural spiking activity was recorded in 13 neurons in the hippocampus of a monkey executing a location scene association task (13). During the experiment, two to four novel scenes were presented along with two to four welllearned scenes in an interleaved random order. Each scene was presented for between 25 and 60 trials. In experiment 4, the data were recorded from 10 neurons in the subthalamic nucleus of human Parkinson's disease patients (Fig. 1D) executing a directed movement task (15). The four movement directions were up, down, left, and right. The neural spike trains were recorded in 10 trials per direction beginning 200 ms before the movement cue and continuing to 200 ms after the cue.
The PP-GLM was fit to the spike trains of each neuron using likelihood analyses as described above. Examples of the model goodness of fit for a neuron from each system is shown in Supporting Information. Examples of the model estimates of the stimulus and history effects for a neuron from each system are shown in Fig. 2.  (Fig. 3, black bars). The higher SNRs (from Eq. 11) in experiments 1 and 2 ( Fig. 3 A and B) are consistent with the fact that the stimuli are explicit, i.e., an electrical current and mechanical displacement of the whisker, respectively, and that the recording sites are only two synapses away from the stimulus. It is also understandable that SNRs are smaller for the hippocampus and thalamic systems in which the stimuli are implicit, i.e., behavioral tasks (Fig. 3 C and D).
A Simulation Study of Single-Neuron SNR Estimation. To analyze the performance of our SNR estimation paradigm, we studied simulated spiking responses of monkey hippocampal neurons with specified stimulus and history dynamics. We assumed four known SNRs of −8.3 dB, −17.4 dB, −28.7 dB, and -∞ dB corresponding, respectively, to stimulus effects on spike rates ranges of 500, 60, 10, and 0 spikes per second (Fig. 4, row 1). For each of the stimulus SNRs, we assumed spike history dependence (Fig. 4, row 2) to be similar to that of the neuron in Fig. 1C. For each of four stimulus effects, we simulated 300 experiments, each consisting of 25 trials (Fig. 4, row 3). To each of the 300 simulated data sets at each SNR level, we applied our SNR estimation paradigm: model  S , the SNR estimates due to the stimulus correcting for the spiking history. The black bars are the 95% bootstrap confidence intervals for SNR dB S . The gray dots areŜNR dB H , the SNR estimates due to the intrinsic biophysics of the neuron correcting for the stimulus. The gray bars are the 95% bootstrap confidence intervals for SNR dB H . The encircled points are the SNR and 95% confidence intervals for the neural spike train raster plots in Fig. 1. fitting, model order selection, goodness-of-fit assessment, and estimation ofŜNR dB S (Fig. 4, row 4) andŜNR dB H (Fig. 4, row 5). The bias-corrected SNR estimates show symmetric spread around their true SNRs, suggesting that the approximate bias correction performed as postulated (Fig. 4, rows 4 and 5). The exception is the case in which the true SNR was −∞ and our paradigm estimatesŜNR dB S as large negative numbers (Fig. 4D,  row 4). TheŜNR dB S are of similar magnitude as the SNR estimates in actual neurons (see SNR = −18.1 dB in the third neuron in Fig.  3C versus −17.4 dB in the simulated neuron (Fig. 4B).

A Simulation Study of SNR Estimation for Single Neurons with No
History Effect. We repeated the simulation study with no spike history dependence for the true SNR values of −1.5 dB, −16.9 dB, −27.9 dB, and -∞ dB, with 25 trials per experiment and 300 realizations per experiment (Fig. 5). Removing the history dependence makes the simulated data within and between trials independent realizations from an inhomogeneous Poisson process. The spike counts across trials within a 1-ms bin obey a binomial model with n = 25 and the probability of a spike defined by the values of the true conditional intensity function times 1 ms. Hence, it is possible to compute analytically the SNR and the bias in the estimates. We used our paradigm to computeŜNR dB S . For comparison, we also computed the variance-based SNR proposed by Lyamzin et al. (20) BothŜNR dB S and the variance-based estimates were computed from the parameters obtained from the same GLM fits (see Eq. S16). For each simulation in Fig. 5, the true SNR value based on our paradigm is shown (vertical lines).
The histograms ofŜNR dB S (Fig. 5, row 3) are spread symmetrically about the true expected SNR. The variance-based SNR estimate overestimates the true SNR in Fig. 5A and underestimates the true SNR in Fig. 5 B and C. These simulations illustrate that the variance-based SNR is a less refined measure of uncertainty, as it is based on only the first two moments of the spiking data, whereas our estimate is based on the likelihood that uses information from all of the moments. At best, the variance-based SNR estimate can provide a lower bound for the information content in the non-Gaussian systems (16). Variance-based SNR estimators can be improved by using information from higher-order moments (37), which is, effectively, what our likelihood-based SNR estimators do.

Discussion
Measuring the SNR of Single Neurons. Characterizing the reliability with which neurons represent and transmit information is an important question in computational neuroscience. Using the PP-GLM framework, we have developed a paradigm for estimating the SNR of single neurons recorded in stimulus response experiments. To formulate the GLM, we expanded the log of the conditional intensity function in a Volterra series (Eq. 4) to represent, simultaneously, background spiking activity, the stimulus or signal effect, and the intrinsic dynamics of the neuron. In the application of the methods to four neural systems, we found that the SNRs of neuronal responses (Eq. 11) to putative stimulisignals-ranged from −29 dB to −3 dB (Fig. 1). In addition, we showed that the SNR of the intrinsic dynamics of the neuron (Eq. 12) was frequently higher than the SNR of the stimulus (Eq. 11). These results are consistent with the well-known observation that, in general, neurons respond weakly to putative stimuli (16,20). Our approach derives a definition of the SNR appropriate for neural spiking activity modeled as a point process. Therefore, it offers important improvements over previous work in which the SNR estimates have been defined as upper bounds derived from Gaussian approximations or using Fano factors and coefficients of variation applied to spike counts. Our SNR estimates are straightforward to compute using the PP-GLM framework (5) and public domain software that is readily available (38). Therefore, they can be computed as part of standard PP-GLM analyses.
The simulation study (Fig. 5) showed that our SNR methods provide a more accurate SNR estimate than recently reported variance-based SNR estimate derived from a local Bernoulli model (20). In making the comparison between the two SNR estimates, we derived the exact prediction error ratios analytically, and we used the same GLM fit to the simulated data to construct the SNR estimates. As a consequence, the differences are only due to differences in the definitions of the SNR. The more accurate  performance of our SNR estimate is attributable to the fact that it is based on the likelihood, whereas the variance-based SNR estimate uses only the first two sample moments of the data. This improvement is no surprise, as it is well known that likelihood-based estimates offer the best information summary in a sample given an accurate or approximately statistical model (34). We showed that for each of the four neural systems, the PP-GLM accurately described the spike train data in terms of goodness-of-fit assessments.
A General Paradigm for SNR Estimation. Our SNR estimation paradigm generalizes the approach commonly used to analyze SNRs in linear Gaussian systems. We derived the generalization by showing that the commonly computed SNR statistic estimates a ratio of EPEs (Supporting Information): the expected prediction of the error of the signal representing the data corrected for the nonsignal covariates relative to the EPE of the system noise. With this insight, we used the work of ref. 26 to extend the SNR definition to systems that can be modeled using the GLM framework in which the signal and relevant covariates can be expressed as separate components of the likelihood function. The linear Gaussian model is a special case of a GLM. In the GLM paradigm, the sum of squares from the standard linear Gaussian model is replaced by the residual deviance (Eq. S10). The residual deviance may be viewed as an estimated KL divergence between data and model (26). To improve the accuracy of our SNR estimator, particularly given the low SNRs of single neurons, we devised an approximate bias correction, which adjusts separately the numerator and the denominator (Eqs. 9 and 10). The bias-corrected estimator performed well in the limited simulation study we reported (Figs. 4 and 5). In future work, we will replace the separate bias corrections for the numerator and denominator with a single bias correction for the ratio, and extend our paradigm to characterize the SNR of neuronal ensembles and those of other non-Gaussian systems.
In Supporting Information, we describe the relationship between our SNR estimate and several commonly used quantities in statistics, namely the R 2 , coefficient of determination, the F statistic, the likelihood ratio (LR) test statistic and f 2 , Cohen's effect size. Our SNR analysis offers an interpretation of the F statistic that is not, to our knowledge, commonly stated. The F statistic may be viewed as a scaled estimate of the SNR for the linear Gaussian model, where the scale factor is the ratio of the degrees of freedom (Eq. S21). The numerator of our GLM SNR estimate (Eq. S9) is a LR test statistics for assessing the strength of the association between data Y and covariates X 2 . The generalized SNR estimator can be seen as generalized effect size. This observation is especially important because it can be further developed for planning neurophysiological experiments, and thus may offer a way to enhance experimental reproducibility in systems neuroscience research (39).
In summary, our analysis provides a straightforward way of assessing the SNR of single neurons. By generalizing the standard SNR metric, we make explicit the well-known fact that individual neurons are noisy transmitters of information.