PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
 
PLoS One. 2011; 6(6): e20490.
Published online 2011 June 7. doi:  10.1371/journal.pone.0020490
PMCID: PMC3110192

Denoising Two-Photon Calcium Imaging Data

Matjaz Perc, Editor

Abstract

Two-photon calcium imaging is now an important tool for in vivo imaging of biological systems. By enabling neuronal population imaging with subcellular resolution, this modality offers an approach for gaining a fundamental understanding of brain anatomy and physiology. Proper analysis of calcium imaging data requires denoising, that is separating the signal from complex physiological noise. To analyze two-photon brain imaging data, we present a signal plus colored noise model in which the signal is represented as harmonic regression and the correlated noise is represented as an order An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e001.jpg autoregressive process. We provide an efficient cyclic descent algorithm to compute approximate maximum likelihood parameter estimates by combing a weighted least-squares procedure with the Burg algorithm. We use Akaike information criterion to guide selection of the harmonic regression and the autoregressive model orders. Our flexible yet parsimonious modeling approach reliably separates stimulus-evoked fluorescence response from background activity and noise, assesses goodness of fit, and estimates confidence intervals and signal-to-noise ratio. This refined separation leads to appreciably enhanced image contrast for individual cells including clear delineation of subcellular details and network activity. The application of our approach to in vivo imaging data recorded in the ferret primary visual cortex demonstrates that our method yields substantially denoised signal estimates. We also provide a general Volterra series framework for deriving this and other signal plus correlated noise models for imaging. This approach to analyzing two-photon calcium imaging data may be readily adapted to other computational biology problems which apply correlated noise models.

Introduction

Two-photon microscopy is now widely recognized as a valuable tool for real-time in vivo imaging of biological systems [1], [2]. A two-photon microscope excites fluorophores in a volume of biological sample using pulsed lasers to induce the emission of a fluorescence signal [3]. Typically, a focused laser beam scans the tissue in a two-dimensional raster pattern, producing a fluorescence image that typically spans hundreds of cells [4]. The images facilitate highly informative and quantitative analyses with a range of biological applications. Two-photon imaging of calcium-sensitive fluorescent indicators to investigate neural physiology is particularly appealing because the measured fluorescence is closely related to neural activity [5]. This imaging modality enables analysis of a broad spatial scale, ranging from the structure of dendritic spines (microns) [6][8] to the architecture of neuronal networks (millimeters) [9][12], as well as analysis of a broad temporal scale from fast action potentials (milliseconds) [10], [12][14] to slow calcium waves (seconds) [15], [16].

Properly separating signal from noise, often termed denoising, is a crucial signal processing procedure in the analysis of imaging data. While there has been considerable success in the development of two-photon microscopy hardware and experimental techniques, the corresponding signal processing methodology has received less attention. The observed fluorescence response depends upon several factors: 1) the nature of the stimulus and the modulation of neural activity due to the stimulus; 2) movements due to highly structured physiological processes; 3) spontaneous neural activity; and 4) optical and electrical noise. Current approaches to processing two-photon data consist of averaging the measured fluorescence levels over multiple trials followed by kernel-based smoothing or fitting an appropriate curve to these time-series data [11], [17][19]. Averaging, while highly intuitive and easy to perform, requires a large number of trials which is often not possible in two-photon imaging experiments. Principal components analysis (PCA) has also been used for denoising image data by dimension reduction [14], but in general it does not exploit the stimulus-driven modulation of the response. One PCA-based approach preserves only the PCA components that exceed a certain threshold of correlation with the stimulus sequence [20]. Fourier analysis of two-photon data recorded in response to periodic stimulation allows for signal extraction at the excitation frequency [21], [22] and possibly some of its harmonics [23], but does not model activity at other frequencies. A signal plus colored noise model has been used to analyze functional magnetic resonance (fMRI) data [24]. However, the expectation-maximization estimation algorithm used in that method has high computational complexity.

Nonlinear approaches to signal analysis could be candidates for the analysis of the data considered here. These include projective filtering, wavelet methods, local linear approximation, and several other methods [25][28]. We derive our approach by assuming a nonlinear relationship between the measured fluorescence response, the stimulus, and the colored noise. The Volterra series expansion leads to an approximation of this nonlinear relationship that relates the measured fluorescence linearly to the harmonic signal and the correlated noise. Our analysis shows that this approach to constructing models for two-photon imaging also yields the commonly used signal plus noise models for fMRI data analysis.

We propose a novel statistical signal plus correlated noise (SCN) model for the analysis of the two-photon calcium imaging data in which the stimulus induced structure is represented as a harmonic regression and the temporally correlated noise is represented as an autoregressive process. We present a computationally efficient cyclic descent algorithm for maximum likelihood estimation of the model parameters. Our approach differs from current two-photon image analysis techniques in that we use a formal likelihood framework to not only separate signal and noise, but to also select a model, assess model goodness of fit and make inferences about physiological features in the estimated images. The high computational efficiency of the algorithm makes it amenable to automated analysis of large imaging data sets. By analyzing two-photon calcium imaging data recorded from the ferret primary visual cortex, we demonstrate that our approach accurately models the data and provides significantly denoised images.

Results

Visual stimulation and two-photon image acquisition

Time series traces of two-dimensional images (XYT) with a field-of-view of approximately An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e002.jpg were collected at An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e003.jpg from the primary visual cortex of a ferret using a two-photon microscope (see Methods). The stimulation protocol consisted of square-wave gratings with An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e004.jpg contrast which drifted at An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e005.jpg orthogonally to the orientation and rotated by An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e006.jpg every second (each data frame), i.e., the stimulus rotated An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e007.jpg in An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e008.jpg. The time series of the response of a neuron to this stimulus approximated a full orientation tuning curve. This stimulus was repeated three times. Prior to recording the stimulus responses, 10 image frames were acquired in the absence of any visual stimulus and their mean provided the estimate of the baseline level at each pixel. Manually determined boundaries delineate the set of pixels that define each cell, and each of the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e009.jpg cells thus identified (Figure S1a) consists of An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e010.jpg pixels (meanAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e011.jpgs.d.). The data consist of the time series of fluorescence on each image pixel. The relative fluorescence on a given pixel is An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e012.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e013.jpg is the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e014.jpg time-sample of the measured fluorescence intensity; An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e015.jpg; An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e016.jpg is the baseline level; and we have An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e017.jpg samples. Using this orientation stimulus, initial anatomical images of the neuronal population can be obtained by plotting the pixel-wise maximum fluorescence across the image time-series, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e018.jpg (Figure 1a). The relative fluorescence traces from the imaged cells show the diversity of orientation responses (Figure 1b and S1b).

Figure 1
Example of two-photon fluorescence image and individual cell dynamics.

A signal plus correlated noise model

We assume that in each pixel, the measured fluorescence at each time can be described by a signal plus correlated noise (SCN) model defined as

equation image
(1)

where the signal is defined as the order An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e024.jpg harmonic regression

equation image
(2)

and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e026.jpg is the period of the stimulus. We assume that the temporally correlated noise obeys the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e027.jpg order autoregressive model (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e028.jpg) given by

equation image
(3)

where the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e030.jpg are assumed to be independent, identically distributed Gaussian random variables with mean zero and unknown variance An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e031.jpg. We assume that the zeros of the characteristic polynomial, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e032.jpg, are outside the unit circle to insure stationarity of the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e033.jpg model. We model the signal as a harmonic regression because the measured fluorescence shows a strong sinusoidal response at the period of the stimulus. We postulate that this smooth, periodic structure should be well described by the low-order terms of a Fourier series expansion defined by the harmonic regression model. The An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e034.jpg model represents the highly structured physiological and electronic noise components of the fluorescence measurements. This and other signal plus correlated noise models can be derived from the Volterra series framework that we present in Methods.

Efficient approximate maximum likelihood parameter estimation by cyclic descent

To use the SCN model in Eq. 1 to denoise calcium imaging data, we estimate its parameters An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e035.jpg, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e036.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e037.jpg by maximum likelihood using a cyclic descent algorithm. The cyclic descent algorithm provides an efficient approach for solving this nonlinear estimation problem by iterating between computing the solutions to two highly tractable linear estimation problems (see Methods). That is, at iteration An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e038.jpg, given An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e039.jpg the estimate of the inverse of the covariance matrix of An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e040.jpg from iteration An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e041.jpg, the algorithm computes An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e042.jpg, the weighted least-squares estimate of An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e043.jpg. Given An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e044.jpg, the algorithm computes An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e045.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e046.jpg using Burg algorithm and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e047.jpg using Levinson-Durbin recursion (Methods). Because of the properties of the Burg algorithm, we provide AR parameter estimates that yield a stationary process. The Levinson-Durbin algorithm provides an efficient means of computing An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e048.jpg from An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e049.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e050.jpg. This efficiency is significant for large An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e051.jpg since An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e052.jpg is a An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e053.jpg matrix. We use as the stopping criterion the condition that the relative change in the estimate of An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e054.jpg between iterations is smaller than threshold An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e055.jpg. If this stopping criterion is satisfied, the algorithm stops; otherwise, given An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e056.jpg, the algorithm proceeds to iteration An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e057.jpg. With this stopping criterion based on the residual variance, the cyclic descent algorithm applied to our calcium imaging data consistently converges in 3 to 5 iterations. This class of iterative algorithms are known to converge at least linearly [29], and our results show that the cyclic descent algorithm in fact achieves exponential convergence (Figure S2). Although this algorithm is highly efficient, we can further expedite processing, as may be required for real-time implementation, by reducing the number of iterations. This cyclic descent algorithm avoids computing the gradients and Hessians required for Newton's procedure and the multiple iterations characteristic of the expectation maximization algorithm. A theorem due to Corradi [29] suggests that our cyclic descent algorithm finds the global maximum of the likelihood.

Choice of model order and assessment of goodness-of-fit

Separation of the fluorescence data into signal and correlated noise relies on choosing appropriate values of model orders An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e058.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e059.jpg. To make these selections, we use well-established model selection and goodness-of-fit criteria, namely the corrected Akaike information criterion (AICc) and analyses of the correlation structure and spectra of the residuals from the model fits (Methods). These criteria help determine the optimal tradeoff between model parsimony and estimation accuracy. As a representative example, we consider the AICc for various model orders for one cell (Figure 2a). The AR model alone can capture much of the periodicity in the data including that due to the stimulus response, but unlike the SCN model, it does not decompose the data into stimulus-driven and background components. Our approach therefore is to fit a signal-only model first and determine the optimal An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e060.jpg using AICc (Figure 2b). Next, using the chosen An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e061.jpg, we fit the SCN model and determine the optimal An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e062.jpg (Figure 2b). Goodness-of-fit analysis is another important consideration whose purpose is to insure that the residuals, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e063.jpg, are white (Methods), confirming that the systematic variance in the data has been explained by the SCN model. We find that inclusion of an AR component is necessary to obtain white residuals even when An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e064.jpg is large (Figure 2c). The AR model order suggested by AICc is insufficient and instead a higher order is required to guarantee white residuals. The fluorescence data spectrum shows certain dominant periodicities, some of which correspond to the stimulus frequency and its low harmonics, and are captured by An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e065.jpg (Figure 2d). The nonuniform spectrum of background activity, including the significant activity observed at low frequencies, is captured by the AR component, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e066.jpg. The spectrum of residuals, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e067.jpg, is approximately flat. The normalized cumulative periodogram (NCP) of An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e068.jpg, falls outside the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e069.jpg whiteness bounds (Figure 2e). In contrast, the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e070.jpg NCP nearly coincides with white noise NCP as desired. This analysis assists in determining the required harmonic and AR model orders, which may vary from cell to cell. We find that An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e071.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e072.jpg satisfy the above requirements for most of the cells in our data-set (Figure S3 and Table S1). Once the optimal SCN model order has been determined and goodness-of-fit assessment completed, we use the model to make biological inferences.

Figure 2
Model order selection and goodness of fit.

Tuning curve estimation

We use the SCN model to characterize the relative fluorescence response to the stimulus at a single pixel. The close fit between the data and the signal estimate establishes the validity of our model (Figure 3a). The signal component, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e090.jpg, provides a denoised estimate of the response for three trials of stimulus presentation (Figure 3b). The autocorrelation function and quantiles of the residual, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e091.jpg, confirm that it is consistent with an independent Gaussian process (Figure 3c and d). We construct the denoised response tuning curve, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e092.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e093.jpg is a circular random variable that represents the stimulus orientation. We also obtain the approximate An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e094.jpg confidence intervals (Methods) and analyze the response characteristics. This signal estimate (Figure 3e) captures the key features of the neuronal response, such as the location and width of tuning to the stimulus effect. The use of a Gaussian or cosine curve to fit the data, as is common practice in neuroscience, would constrain the response estimate to have a simple, symmetric shape. Our model allows the tuning curve estimate to reflect the complex shape of the cell response observed in the data with minimal computational complexity.

Figure 3
Decomposition of fluorescence data into signal and noise components.

Image denoising

We can reconstruct denoised images using the signal component estimate, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e100.jpg, at each pixel. A comparison of the fluorescence response estimates of pixels around a cell obtained with conventional across-trial averaging and with our SCN model (Figure 4a) demonstrates the enhanced image contrast and clarity provided by our model. Our denoising method delineates the stimulus response within the cell soma and allows improved observation of calcium dynamics around the cell associated with excitation. In a second cell (Figure 4b), the background activity at the bottom of the frame is substantially reduced in intensity by our approach. The increased contrast of the denoised images reveals additional subcellular processes not discernible in the conventional images obtained by averaging (Figure 4c). This opens up the possibility of future work to characterize the source of these signals and their behavior.

Figure 4
Image denoising and visualization of calcium activity.

Inferring neuronal characteristics

By denoising two-photon imaging data with the SCN model, we provide reliable estimates of several quantities that may be used to describe neuronal behavior. For example, we study the orientation preferences of the primary visual cortex neurons in our sample. At each pixel, the preferred orientation is obtained as the orientation at which the denoised tuning curve peak occurs, i.e., An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e104.jpg. As previously reported [9], neighboring cells show a preference for similar orientations with a smooth spatial variation (Figure 5a and c). Among the cells, there are different degrees of deviation from the mean preferred orientation (Figure 5d). This deviation is particularly high for two of the cells possibly due to somatic and dendritic dynamics. We calculate the orientation selectivity from the estimated tuning curve, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e105.jpg, as the half-width at half-height. Our analysis of orientation selectivity at each pixel reveals both spatial trends and intra-cellular variations (Figure 5b). A wide-ranging level of orientation selectivity is apparent (Figure 5e). These examples demonstrate that the SCN model can facilitate a variety of functional analyses with high reliability.

Figure 5
Spatial distribution of preferred orientation and orientation selectivity after denoising with SCN model.

Neuronal signal-to-noise ratio estimation

The ratio of stimulus-evoked response (signal) to background activity (colored noise) provides a natural definition of the neuronal signal-to-noise ratio (SNR) and a way to compare the relative responsiveness of the cells to the stimulus. We calculate the signal power, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e107.jpg from our harmonic model and the noise power, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e108.jpg, from the AR model to obtain

equation image
(4)

The cells in our data set exhibit a wide range of SNRs (Figure 6a). The locations of the cells with high SNR (Figure 6b) agree closely with the anatomical map (Figure 1a), and therefore the pixel-wise SNR maps can be used to identify robustly the locations of cells that respond to the given stimulation.

Figure 6
Neuronal signal-to-noise ratio (SNR) estimated using the SCN model.

Discussion

We have presented a flexible local likelihood framework for analyzing two-photon calcium imaging data. Our framework appreciably enhances image contrast on a pixel-by-pixel basis by using an SCN model to separate the salient stimulus-evoked neural responses from the complex forms of physiological and recording noise common to two-photon imaging experiments. The cyclic descent algorithm provides a computationally efficient approach for fitting the SCN model to the time-series of fluorescence responses. Our framework suggests a straight-forward yet novel way to track with improved subcellular resolution the temporal dynamics of individual neurons (Figure 4) and obtain significantly denoised images of neuronal populations (Figure 5).

We have formulated our analysis as a harmonic regression plus colored noise problem. Cellular calcium responses have a stochastic nature and exhibit oscillations with a colored noise component [30], as demonstrated by calcium recordings from pancreatic acinar cells [31] and airway myocites [32]. Colored noise also appears in many other contexts in computational biology, including functional magnetic resonance imaging (fMRI) [24], [33], neural voltage-sensitive dye imaging [34], circadian rhythms [35], synaptic background activity in cortical neurons [36][38], gene regulatory networks [39], speech signals [40], cell locomotion patterns [41], and many others. The procedures usually applied to these problems, based on expectation maximization or exact maximum-likelihood procedures, are often computationally intensive. Our approach suggests an alternate approximate maximum likelihood procedure that can be applied to a broad range of such problems, and may offer specific advantages for analysis in real-time computation and high-throughput processing.

We developed our analysis framework using a continuous and periodic stimulus applied to visual cortex cells. Implicit in our approach is the treatment of intrinsic imaging, as the signal dynamics can be easily described using our model. Also, with minor modifications, our framework can be easily extended to imaging protocols using other stimuli. For example, in some two-photon imaging experiments the stimulus is applied either in a random noise-like manner to avoid anticipatory responses [42], or by interspersing blank frames with no relevant excitatory or inhibitory stimulus present [9]. To apply our approach to data recorded from any of these experiments, we simply replace the stimulus represented as a harmonic regression in the current SCN model with an appropriate formulation of the stimulus model for the given protocol. The remainder of the analysis paradigm, including model fitting, model selection, goodness-of-fit assessment and inference, then proceeds as described above.

Our analytical framework is general so that a number of current statistical models that describe imaging modalities can be easily derived from it. For example, a commonly used model for fMRI data analysis [43] can be obtained from our Volterra series framework (see Methods) as

equation image
equation image
(5)

where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e115.jpg is a gamma function used to model the hemodynamic response of the body. The second term on the right represents physiological noise.

We assumed a signal plus Gaussian noise model in our analysis. However, in two-photon microscopy and other optical imaging modalities, the measured fluorescence intensity is a function of the discrete number of incident photons, and is therefore fundamentally a counting process and not necessarily Gaussian. The counting process nature of the two-photon experiments becomes more apparent as the acquisition rate increases [44]. The measured fluorescence in some two-photon imaging experiments may also exhibit non-Gaussian behavior due to distortions introduced during acquisition or post-processing. If the Gaussian assumption no longer holds, we can develop alternative likelihood approaches based on appropriately chosen non-Gaussian models. For example, neuronal spike trains can be extracted from two-photon data using template deconvolution [45], [46]. Non-Gaussian likelihoods based on the theory of point processes and implemented using the generalized linear model could be adapted to analyze these two-photon imaging data.

We modeled the time-series of neural responses in each pixel separately and did not consider inter-pixel dependencies. Such dependencies arise because: the activity of a single cell is captured across multiple pixels; retinotopy and network dependencies may lead to similar response in contiguous regions of the image; and data pre-processing procedures such as spatial smoothing introduce correlations. This problem, although challenging, currently confronts all biological imaging modalities and should ideally be studied by formulating appropriate biologically-based spatiotemporal models.

We have illustrated the application of our framework in offline analyses. However, due to its low computational complexity, our current analysis paradigm can be readily adapted to conduct large-volume, high-throughput imaging data analyses in real-time. These and other related aspects will be the topics of future reports.

Methods

Experimental procedures

Two-photon imaging of the fluorescent calcium indicator Oregon Green Bapta (OGB) was performed in the visual cortex of anesthetized ferrets. Neurons were bulk-loaded with OGB by intracortical injection of the AM-ester conjugated form of OGB using standard techniques [9], [18], [47]. Imaging was performed with a custom-made two-photon laser scanning microscope consisting of a modified Olympus Fluoview confocal scan head and a titanium/sapphire laser providing approximately An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e116.jpg fsec pulses at An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e117.jpg MHz pumped by a An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e118.jpg W solid-state source [48]. Fluorescence was detected using photomultiplier tubes in whole-field detection mode and a An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e119.jpg, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e120.jpg NA lens. Image acquisition was accomplished using Fluoview software. The images were taken from cortical layers 2/3, and this area was readily distinguished from layer 1 on the basis of the relative density of astrocytes and neurons. Visual stimuli, generated with Matlab using the PsychoPhysics Toolbox [49], were delivered via a An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e121.jpg LCD display placed An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e122.jpg m away from the eyes of the animal. Neurons with relative fluorescence clearly distinguishable from the neuropil were chosen for subsequent cellular analysis.

A Volterra series framework for signal plus colored noise imaging models

We assume that the measured fluorescence, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e123.jpg, in a two-photon calcium imaging experiment is a function of a time-varying stimulus, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e124.jpg, and noise in the system, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e125.jpg. We further assume that the response, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e126.jpg, of the biological system depends on a nonlinear transformation of the stimulus input to the biological system. We can express the effect of the input stimulus and noise on the measured fluorescence at a pixel as

equation image
(6)

Expanding the right side of Eq. 6 in a Volterra series [50] yields

equation image
equation image
equation image
(7)

Take a discrete approximation to the first two terms on the right of Eq. 7 and assume that the second-order terms are sufficiently small so that they can be approximated as An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e131.jpg, independent Gaussian noise with mean zero and variance An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e132.jpg. Then Eq. 7 becomes

equation image
(8)

where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e134.jpg denotes the time-samples. We can express An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e135.jpg in terms of its Fourier series expansion. If An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e136.jpg is periodic and smooth, its Fourier expansion can be well represented by a finite series. Thus, taking the first An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e137.jpg terms of the series, we can write

equation image
(9)

where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e139.jpg is the period of the input stimulus. In the two-photon imaging experiment, we assume that the effect of the stimulus on the system is instantaneous so that the discrete kernel can be written in terms of the Kronecker delta function as

equation image
(10)

Substituting Eq. 10 into Eq. 9 yields

equation image
(11)

which is our model given in Eq. 1. If our model had not fit the data well, then we could use Eq. 7 to derive a modified model by including one or more of the second order terms.

Burg algorithm

The Burg algorithm for autoregressive (AR) coefficient estimation uses least squares forward-backward prediction error minimization and is constrained to satisfy Levinson-Durbin recursions (LDR) [51], [52]. For the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e142.jpg model in Eq. 3, the Burg algorithm estimates the coefficients An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e143.jpg and innovations variance An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e144.jpg, given the time series An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e145.jpg for An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e146.jpg, as follows [52].

Require: An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e147.jpgAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e148.jpgAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e149.jpgAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e150.jpg

  forAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e151.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e152.jpg doAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e153.jpgAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e154.jpg

  if An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e155.jpg then

   for An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e156.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e157.jpg doAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e158.jpg

   end for

   for An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e159.jpg toAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e160.jpg doAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e161.jpg

   end for

   for An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e162.jpg to An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e163.jpg doAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e164.jpg

   end for

  end for

end forAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e165.jpg

return An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e166.jpg

Cholesky factorization

The An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e167.jpg covariance matrix of the AR process An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e168.jpg can be written in its Cholesky form as An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e169.jpg. The inverse matrix An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e170.jpg is used in our cyclic descent algorithm and can be calculated efficiently using Levinson-Durbin recursions, where [52]

equation image
(12)

and

equation image
(13)

The coefficient and variance estimates of AR models up to order An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e173.jpg are computed by the Burg algorithm during An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e174.jpg model parameter estimation and are therefore already available to populate An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e175.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e176.jpg. Hence this is a highly efficient procedure for computing An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e177.jpg without matrix inversion.

Cyclic descent algorithm

We use the cyclic descent algorithm for joint estimation of autoregressive and harmonic coefficient vectors, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e178.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e179.jpg, from the fluorescence data vector An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e180.jpg. The algorithm proceeds as follows.

equation image
equation image
equation image

repeat

equation image
equation image
equation image

Compute An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e187.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e188.jpg from An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e189.jpg using Burg algorithm.

Compute An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e190.jpg from its Cholesky factors using Levinson-Durbin recursion untillAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e191.jpg

equation image

returnAn external file that holds a picture, illustration, etc.
Object name is pone.0020490.e193.jpg

Statistical methods

As our SCN model consists of two linear components, it is straight-forward to obtain confidence bounds and construct significance tests for the model parameters. For the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e194.jpg harmonic regression coefficient estimate, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e195.jpg, the approximate 95% confidence interval is An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e196.jpg where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e197.jpg; An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e198.jpg; and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e199.jpg [53]. Similarly, for the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e200.jpg coefficient of the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e201.jpg process, we have the confidence interval An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e202.jpg where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e203.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e204.jpg. Here, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e205.jpg contains the time-lagged samples of the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e206.jpg process and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e207.jpg. Based on these confidence intervals, we can design the t-test of significance for the model coefficients. The alternate hypothesis for the significance of the harmonic coefficients is An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e208.jpg, and the alternative hypothesis for AR coefficient significance is defined similarly. The corresponding parameter is significant if the hypothesis is rejected.

We use the corrected Akaike information criterion (AICc) for model order selection. For the order An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e209.jpg harmonic regression and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e210.jpg model, we define An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e211.jpg and

equation image
(14)

We use residual analysis to confirm the whiteness of our model's residual noise, An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e213.jpg. The normalized autocorrelation function (ACF) of the residuals at lag An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e214.jpg is given by An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e215.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e216.jpg. The approximate 95% whiteness bounds are An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e217.jpg and the corresponding Ljung-Box portmanteau test statistic is An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e218.jpg, where conventionally An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e219.jpg ACF taps are considered. The null hypothesis for the whiteness test is An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e220.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e221.jpg denotes the alpha level, taken as 5% in our analysis.

We use circular statistics to analyze circular random variables such as orientation An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e222.jpg. The circular mean is calculated as An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e223.jpg with its 95% confidence interval given by An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e224.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e225.jpg is circular standard error and

equation image
(15)

is circular dispersion [54].

Ethics Statement

All experimental procedures were approved by the Massachusetts Institute of Technology Committee on Animal Care and adhered to the NIH guidelines for the Care and Use of Laboratory Animals (Protocol No. 0608-069-11).

Supporting Information

Figure S1

Two-photon fluorescence image of a cell population. (a) Anatomical image of a population of 15 cells. Brighter gray shades represent higher fluorescence intensity. ROIs and cell indices indicate all of the cells identified manually. (b) Orientation tuning curve of each cell obtained by averaging the measured relative fluorescence across three trials.

(TIF)

Figure S2

Convergence of the parameter estimates obtained with cyclic descent. (a) Iterative estimates of the model parameters, namely the residual variance (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e227.jpg), intercept (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e228.jpg), harmonic coefficients (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e229.jpg and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e230.jpg; An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e231.jpg), autoregressive coefficients (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e232.jpg; An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e233.jpg) and residual variance (An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e234.jpg, for the fluorescence time series in Figure 3a. (b) Percentage difference between the successive estimates of the model parameters in a. For the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e235.jpg iteration, the percentage difference is calculated as An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e236.jpg, where An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e237.jpg is the estimate of parameter An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e238.jpg at the An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e239.jpg iteration and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e240.jpg. The y-axis has a logarithmic scale.

(TIF)

Figure S3

Model order selection. (a) AICc surface for each of the cells in our data set averaged across the pixels for that cell (blue: low; red: high). (b) AICc as a function of the harmonic order when an AR component is not fit to the residual of the harmonic regression. (c) AICc as a function of the AR order when the optimal harmonic order from b is used. (d) Percentage of pixels of each cell that pass the Ljung-Box whiteness test (blue: An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e241.jpg; red: An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e242.jpg).

(TIF)

Table S1

Optimal harmonic and AR model orders. The table shows the optimal model orders for each cell, obtained using AICc and Ljung-Box test. Based on these results, we conclude that for our data set, a good fit to the data is obtained with approximately An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e243.jpg harmonics and An external file that holds a picture, illustration, etc.
Object name is pone.0020490.e244.jpg AR coefficients.

(DOC)

Footnotes

Competing Interests: The authors have declared that no competing interests exist.

Funding: This work was supported by U.S. National Institutes of Health DP1 OD003646-01, R01 EB006385-01, EY07023 and EY017098. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

1. Wilt B, Burns LD, Ho ETW, Ghosh KK, Mukamel EA, et al. Advances in light microscopy for neuroscience. Annu Rev Neurosci. 2009;32:435–506. [PMC free article] [PubMed]
2. Denk W, Svoboda K. Photon upmanship: why multiphoton imaging is more than a gimmick. Neuron. 1997;18:351–357. [PubMed]
3. Homma R, Baker BJ, Jin L, Garaschuk O, Konnerth A, et al. Dynamic Brain Imaging: Multi-Modal Methods and In Vivo Applications, Humana Pres., chapter Wide-field and two-photon imaging of brain activity with voltage- and calcium-sensitive dyes. 2009. [PubMed]
4. Helmchen F, Denk W. Deep tissue two-photon microscopy. Nature Methods. 2005;2:932–940. [PubMed]
5. Kerr JND, Denk W. Imaging in vivo: watching the brain in action. Nature Neuroscience. 2008;9:195–205. [PubMed]
6. Majewska A, Sur M. Motility of dendritic spines in visual cortex in vivo: Changes during the critical period and effects of visual deprivation. Proc Natl Acad Sci. 2003;100:16024–16029. [PubMed]
7. Barretto RPJ, Messerschmidt B, Schnitzer MJ. In vivo uorescence imaging with highresolution microlenses. Nature Methods. 2009;6:511–512. [PMC free article] [PubMed]
8. Holtmaat A, Bonhoeffer T, Chow DK, Chuckowree J, Paola VD, et al. Long-term, highresolution imaging in the mouse neocortex through a chronic cranial window. Nat Protoc. 2009;4:1128–1144. [PMC free article] [PubMed]
9. Ohki K, Chung S, Ch'ng YH, Kara P, Reid CR. Functional imaging with cellular resolution reveals precise micro-architecture in visual cortex. Nature. 2005;433:597–603. [PubMed]
10. Göbel W, Kampa BM, Helmchen F. Imaging cellular network dynamics in three dimensions using fast 3D laser scanning. Nature Methods. 2006;4:73–79. [PubMed]
11. Kerr JND, Greenberg D, Helmchen F. Imaging input and output of neocortical networks in vivo. Proc Natl Acad Sci. 2005;102:14063–14068. [PubMed]
12. Greenberg DS, Houweling AR, Kerr JND. Population imaging of ongoing neuronal activity in the visual cortex of awake rats. Nature Neuroscience. 2008;11:749–751. [PubMed]
13. Yaksi E, Friedrich RW. Reconstruction of firing rate changes across neuronal populations by temporally deconvolved Ca2+ imaging. Nature Methods. 2006;3:377–383. [PubMed]
14. Mukamel EA, Nimmerjahn A, Schnitzer MJ. Automated analysis of cellular signals from large-scale calcium imaging data. Neuron. 2009;63:747–760. [PMC free article] [PubMed]
15. Garaschuk O, Linn J, Eilers J, Konnerth A. Large-scale oscillatory calcium waves in the immature cortex. Nature. 2000;3:452–459. [PubMed]
16. Mironov SL. Metabotropic glutamate receptors activate dendritic calcium waves and TRPM channels which drive rhythmic respiratory patterns in mice. J Physiol. 2008;586:2277–2291. [PubMed]
17. Basole A, White LE, Fitzpatrick D. Mapping multiple features in the population response of visual cortex. Nature. 2003;423:986–990. [PubMed]
18. Schummers J, Yu H, Sur M. Tuned responses of astrocytes and their inuence on hemodynamic signals in the visual cortex. Science. 2008;320:1638–1643. [PubMed]
19. Martinez LM, Alonso JM, Reid CR, Hirsch JA. Laminar processing of stimulus orientation in cat visual cortex. J Physiol. 2002;540:321–333. [PubMed]
20. Gabbay M, Brennan C, Kaplan E, Sirovich L. A principal components-based method for the detection of neuronal activity maps: application to optical imaging. NeuroImage. 2000;1:313–325. [PubMed]
21. Kalatsky VA, Stryker MP. New paradigm for optical imaging: temporally encoded maps of intrinsic signal. Neuron. 2003;38:529–545. [PubMed]
22. Mrsic-Flogel T, Hübener M, Bonhoeffer T. Brain mapping: New wave optical imaging. Curr Biol. 2003;13:R778–R780. [PubMed]
23. Sornborger A, Sailstad C, Kaplan E, Sirovich L. Spatiotemporal analysis of optical imaging data. NeuroImage. 2003;18:610–621. [PubMed]
24. Purdon PL, Solo V, Weisskoff RM, Brown EN. Locally regularized spatiotemporal modeling and model comparison for functional MRI. NeuroImage. 2001;14:912–923. [PubMed]
25. Kantz H, Schreiber T. Cambridge University Press; 2004. Nonlinear Time Series Analysis.
26. Kantz H, Schreiber T. Nonlinear noise reduction: A case study on experimental data. Phys Rev E. 1993;48 [PubMed]
27. Schreiber T, Grassberger A simple noise-reduction method for real data. Phys Lett A. 1991;160:411–418.
28. Effern A, Lehnertz K, Schreiber T, Grunwald T, David P, et al. Nonlinear denoising of transient signals with application to event-related potentials. Physica D. 2000;140:257–266.
29. Corradi C. A note on the computation of maximum likelihood estimates in linear regression models with autocorrelated errors. J Econometrics. 1979;11:303–317.
30. Perc M, Green AK, Marhl M. Establishing the stochastic nature of intracellular calcium oscillations from experimental data. Biophys Chem. 2008;132:33–38. [PubMed]
31. Perc M, Rupnik M, Gosak M, Marhl M. Prevalence of stochasticity in experimentally observed responses of pancreatic acinar cells to acetylcholine. Chaos. 2009;19:037113. [PubMed]
32. Marhl M, Gosak M, Perc M, Roux E. Importance of cell variability for calcium signaling in rat airway myocytes. Biophys Chem. 2010;148:42–50. [PubMed]
33. Bullmore E, Long C, Suckling J, Fadili J, Calvert G, et al. Colored noise and computational inference in neurophysiological (fMRI) time series analysis: Resampling methods in time and wavelet domains. Human Brain Mapping. 2001;12:61–78. [PubMed]
34. Chen Y, Geisler WS, Seidemann E. Optimal temporal decoding of neural population responses in a reaction-time visual detection task. J Neurophysiol. 2008;99:1366–1379. [PMC free article] [PubMed]
35. Brown EN, Choe Y, Luithardt H, Czeisler CA. A statistical model of the human coretemperature circadian rhythm. Am J Physiol Endocrinol Metab. 2000;279:E669–E683. [PubMed]
36. Brunel N, Latham PE. Firing rate of the noisy quadratic integrate-and-fire neuron. Neural Comput. 2003;15:2281–2306. [PubMed]
37. Destexhe A, Rudolph M, Fellous JM, Sejnowski TJ. Fluctuating synaptic conductances recreate in vivo-like activity in neocortical neurons. Neuroscience. 2001;107:13–24. [PMC free article] [PubMed]
38. Sakai Y, Funahashi S, Shinomoto S. Temporally correlated inputs to leaky integrate-and-fire models can reproduce spiking statistics of cortical neurons. Neural Networks. 1999;12:1181–1190. [PubMed]
39. Huang MC, Wu JW, Luo YP, Petrosyan KG. Fluctuations in gene regulatory networks as Gaussian colored noise. J Chem Phys. 2010;132:155101. [PubMed]
40. Gibson JD, Koo B, Gray SD. Filtering of colored noise for speech enhancement and coding. IEEE Trans Sig Proc. 1991;39:1732–1742.
41. Shenderov AD, Sheetz MP. Inversely correlated cycles in speed and turning in an ameba: An oscillatory model of cell locomotion. Biophys J. 1997;72:2382–2389. [PubMed]
42. Everson RM, Prashanth AK, Gabbay M, Knight BW, Sirovich L, et al. Representation of spatial frequency and orientation in the visual cortex. Proc Natl Acad Sci. 1998;95:8334–8338. [PubMed]
43. Friston KJ, Josephs O, Rees G, Turner R. Nonlinear event-related responses in fMRI. Magn Reson Med. 1998;39:41–52. [PubMed]
44. Sjulson L, Miesenböck G. Optical recording of action potentials and other discrete physiological events: a perspective from signal detection theory. Physiology. 2007;22:47–55. [PubMed]
45. Vogelstein JT, Watson BO, Packer AM, Yuste R, Jedynak B, et al. Spike inference from calcium imaging using sequential Monte Carlo methods. Biophys J. 2009;97:636–655. [PubMed]
46. Ozden I, Lee HM, Sullivan MR, Wang SH. Identification and clustering of event patterns from in vivo multiphoton optical recordings of neuronal ensembles. J Neurophysiol. 2008;100:495–503. [PubMed]
47. Stosiek C, Garaschuk O, Holthoff K, Konnerth In vivo two-photon calcium imaging of neuronal networks. Proc Natl Acad Sci. 2003;100:7319–7324. [PubMed]
48. Majewska A, Yiu G, Yuste R. A custom-made two-photon microscope and deconvolution system. Eur J Physiol. 2000;441:398–408. [PubMed]
49. Brainard DH. The Psychophysics Toolbox. Spatial Vision. 1997;10:433–436. [PubMed]
50. Schetzen M. Krieger; 2006. The Volterra and Wiener Theories of Nonlinear Systems.
51. Box GEP, Jenkins GM, Reinsel GC. Wiley; 2008. Time Series Analysis.
52. Kay SM. Prentice Hall; 1999. Modern Spectral Estimation.
53. Myers RH, Montgomery DC, Vining GG. Wiley; 2002. Generalized Linear Models.
54. Fisher NI. Cambridge; 1993. Statistical Analysis of Circular Data.

Articles from PLoS ONE are provided here courtesy of Public Library of Science