Denoising Two-Photon Calcium Imaging Data

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 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 twophoton 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 twodimensional 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][7][8] to the architecture of neuronal networks (millimeters) [9][10][11][12], as well as analysis of a broad temporal scale from fast action potentials (milliseconds) [10,[12][13][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][18][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][26][27][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.

Visual stimulation and two-photon image acquisition
Time series traces of two-dimensional images (XYT) with a field-of-view of approximately 250|250mm were collected at 1Hz from the primary visual cortex of a ferret using a two-photon microscope (see Methods). The stimulation protocol consisted of square-wave gratings with 100% contrast which drifted at 3Hz orthogonally to the orientation and rotated by 10 0 every second (each data frame), i.e., the stimulus rotated 360 0 in 36sec. 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 n~15 cells thus identified ( Figure S1a) consists of 79+25 pixels (mean+s.d.). The data consist of the time series of fluorescence on each image pixel. The relative fluorescence on a given pixel is Df k =f~f k {f 0 ð Þ=f 0 , where f k is the k th time-sample of the measured fluorescence intensity; k~1, . . . ,K; f 0 is the baseline level; and we have K~108 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, max k Df k f g (Figure 1a). The relative fluorescence traces from the imaged cells show the diversity of orientation responses (Figure 1b and S1b).

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 where the signal is defined as the order h harmonic regression and t is the period of the stimulus. We assume that the temporally correlated noise obeys the p th order autoregressive model (AR(p)) given by where the e k are assumed to be independent, identically distributed Gaussian random variables with mean zero and unknown variance s 2 . We assume that the zeros of the characteristic polynomial, 1{ X p j~1 a j z {j , are outside the unit circle to insure stationarity of the AR(p) 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 AR(p) 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 b~m, a 1 , b 1 , . . . ,a h , b h ð Þ , a~a 1 , . . . , a p À Á and s 2 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 n, fromâ a (n) andŝ s 2(n) . This efficiency is significant for large K since W is a K|K matrix. We use as the stopping criterion the condition that the relative change in the estimate ofŝ s 2 between iterations is smaller than threshold g. If this stopping criterion is satisfied, the algorithm stops; otherwise, givenŴ {1(n) , the algorithm proceeds to iteration nz1. 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 h and p. 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 h using AICc ( Figure 2b). Next, using the chosen h, we fit the SCN model and determine the optimal p ( Figure 2b). Goodness-of-fit analysis is another important consideration whose purpose is to insure that the residuals,ê e k , 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 h 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ŝ s k (Figure 2d). The nonuniform spectrum of background activity, including the significant activity observed at low frequencies, is captured by the AR component,v v k . The spectrum of residuals,ê e k , is approximately flat. The normalized cumulative periodogram (NCP) ofv v k , falls outside the 95% whiteness bounds (Figure 2e). In contrast, theê e k 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 h~4 and p~10 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.

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,ŝ s k , provides a denoised estimate of the response for three trials of stimulus presentation ( Figure 3b). The autocorrelation function and quantiles of the residual,ê e k , confirm that it is consistent with an independent Gaussian process ( Figure 3c and d). We construct the denoised response tuning curve,û u(w), where w is a circular random variable that represents the stimulus orientation. We also obtain the approximate 95% 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.

Image denoising
We can reconstruct denoised images using the signal component estimate,ŝ s k , 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.

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., w w~arg max w fû u(w)g. 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,û u(w), 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.

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, P s from our harmonic model and the noise power, P v , from the AR model to obtain 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.

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][37][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. Relative fluorescence data f k (blue) in three consecutive trials and estimate (red) of the signal componentŝ s k (i.e., the stimulus-evoked activity). (c) Autocorrelation function (red) of residual noise,ê e k , lies within the 95% whiteness bounds (blue). (d) The quantile-quantile plot of the residuals confirms Gaussianity. The results in (c) and (d) prove that the residuals are independently and identically distributed Gaussian, and the systematic variance in the data has been explained by the harmonic regression and autoregressive terms. (e) Orientation tuning curve obtained from the denoised signal estimate in (b). The SCN model provides a smooth fit to the across-trials mean of the data. Point-wise approximate 95% confidence intervals are also shown. The SCN model preserves the complex, asymmetric shape of the response tuning curve. doi:10.1371/journal.pone.0020490.g003 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 where K g (t)~t a e bt 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 postprocessing. 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 largevolume, high-throughput imaging data analyses in real-time. These and other related aspects will be the topics of future reports.

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 100 fsec pulses at 80 MHz pumped by a 10 W solid-state source [48]. Fluorescence was detected using photomultiplier tubes in whole-field detection mode and a 20|, 0:95 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 17 00 LCD display placed 0:15 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, f (t), in a twophoton calcium imaging experiment is a function of a time-varying stimulus, g(t), and noise in the system, v(t). We further assume that the response, s(t), 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 Expanding the right side of Eq. 6 in a Volterra series [50] yields 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 e k , independent Gaussian noise with mean zero and variance s 2 . Then Eq. 7 becomes where k~1, . . . , K denotes the time-samples. We can express s(g k ) in terms of its Fourier series expansion. If s(g k ) is periodic and smooth, its Fourier expansion can be well represented by a finite series. Thus, taking the first h terms of the series, we can write where t 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 Substituting Eq. 10 into Eq. 9 yields 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 AR(p) model in Eq. 3, the Burg algorithm estimates the coefficients fa 1 , . . . ,a p g and innovations variance s 2 , given the time series v k for k~1, . . . ,K, as follows [52].
Require: E v k f g~0 end for for k~nz1 toK do

Cyclic descent algorithm
We use the cyclic descent algorithm for joint estimation of autoregressive and harmonic coefficient vectors,â andb, from the fluorescence data vector f~f k ½ k~1,...,K . The algorithm proceeds as follows.
Computeâ a (n) andŝ s 2(n) fromv (n) using Burg algorithm.  We use the corrected Akaike information criterion (AICc) for model order selection. For the order h harmonic regression and AR(p) model, we define AIC~K log s 2 z2(2hzpz1) and We use residual analysis to confirm the whiteness of our model's residual noise, e k . The normalized autocorrelation function (ACF) of the residuals at lag t is given by r t (e)~d t =d 0 , where d t~K , where conventionally T~20 ACF taps are considered. The null hypothesis for the whiteness test is H 0~Q *x 2 a,T{p , where a denotes the alpha level, taken as 5% in our analysis.
We use circular statistics to analyze circular random variables such as orientation w. The circular mean is calculated as w w~tan {1 X N i~1 sin w i = 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).