Progress in modern neuroscience critically depends on our ability to observe the activity of large neuronal populations with cellular spatial and high temporal resolution. However, two bottlenecks constrain efforts towards fast imaging of large populations. First, the resulting large video data is challenging to analyze. Second, there is an explicit tradeoff between imaging speed, signal-to-noise, and field of view: with current recording technology we cannot image very large neuronal populations with simultaneously high spatial and temporal resolution. Here we describe multi-scale approaches for alleviating both of these bottlenecks. First, we show that spatial and temporal decimation techniques based on simple local averaging provide order-of-magnitude speedups in spatiotemporally demixing calcium video data into estimates of single-cell neural activity. Second, once the shapes of individual neurons have been identified at fine scale (e. g. , after an initial phase of conventional imaging with standard temporal and spatial resolution), we find that the spatial/temporal resolution tradeoff shifts dramatically: after demixing we can accurately recover denoised fluorescence traces and deconvolved neural activity of each individual neuron from coarse scale data that has been spatially decimated by an order of magnitude. This offers a cheap method for compressing this large video data, and also implies that it is possible to either speed up imaging significantly, or to “zoom out” by a corresponding factor to image order-of-magnitude larger neuronal populations with minimal loss in accuracy or temporal resolution.
A major goal of neuroscience is to understand interactions within large populations of neurons, including their network dynamics and emergent behavior. This ideally requires the observation of neural activity over large volumes. Recently, light-sheet microscopy and genetically encoded indicators have enabled unprecedented whole-brain imaging of tens of thousands of neurons at cellular resolution [1]. However, light-sheet microscopy generally suffers from slow volumetric speeds (e. g. [2], but see also [3,4]) and is usually applied to small and transparent brains. In scattering brains, current technologies with single-neuron resolution are usually based on slow, serially-scanned two-photon (2P) imaging methods that can only sample from O (102 − 103) neurons simultaneously with adequate temporal resolution [5]. Recent advances have enabled faster light-sheet imaging in cortex [6] and fast volumetric 2P imaging [7], but we must still contend with critical trade-offs between temporal and spatial resolution—and the need for even faster imaging of even larger neural populations. Another critical challenge is the sheer amount of data generated by these large-scale imaging methods. A crucial step for further neural analysis involves a transition from voxel-space to neuron-source space: i. e. , we must detect the neurons and extract and demix each neuron’s temporal activity from the video. Simple methods such as averaging voxels over distinct regions of interest (ROIs) are fast, but more statistically-principled methods based on constrained non-negative matrix factorization (CNMF) better conserve information, yield higher signal-to-noise ratio, recover more neurons, and enable the demixing of spatially overlapping neurons [8]. The methods described in [8] were not optimized for very large datasets, but NMF is a key machine learning primitive that has enjoyed more than a decade of intensive algorithmic optimization [9–12] that we can exploit here to scale the CNMF approach. We find that a very simple idea leads to order-of-magnitude speedups: by decimating the data (i. e. , decreasing the resolution of the data by simple local averaging [13]), we can obtain much faster algorithms with minimal loss of accuracy. Decimation ideas do not just lead to faster computational image processing, but also offer prescriptions for faster image acquisition over larger fields of view (FOV), and for observing larger neural populations. Specifically, we propose the following two-phase combined image acquisition/analysis approach. In the first phase, we use conventional imaging methods to obtain estimates of the visible neuronal locations and shapes. After this cell-identification phase is complete we switch to low-spatial-resolution imaging, which in the case of camera-based imaging simply corresponds to “zooming out” on the image, i. e. , expanding the spatial size of each voxel. This has the benefit of projecting a larger FOV onto the same number of voxels; alternatively, if the number of voxels recorded per second is a limiting factor, then recording fewer (larger) voxels per frame implies that we can image at higher frame-rates. We are thus effectively trading off spatial resolution for temporal resolution; if we cut the spatial resolution too much we may no longer be able to clearly identify or resolve single cells by eye in the obtained images. However, we show that, given the high-spatial-resolution information obtained in the first imaging phase, the demixing stage of CNMF can recover the temporal signals of interest even from images that have undergone radical spatial decimation (an order of magnitude or more). In other words, CNMF significantly shifts the tradeoff between spatial and temporal resolution, enabling us to image larger neuronal populations at higher temporal resolution. The rest of this paper is organized as follows. We first describe how temporal and spatial decimation (along with several other improvements) can be used within the CNMF algorithm to gain order-of-magnitude speed-ups in calcium imaging video processing. Next we investigate how decimation can enable faster imaging of larger populations for light-sheet and 2P imaging. We show the importance of the initial cell identification phase, quantitatively illustrate how CNMF changes the tradeoff between spatial and temporal resolution, and discuss how spatially decimated imaging followed by demixing can be interpreted as a simple compression and decoding scheme. We show that good estimates of the neural shapes can be obtained on a small batch of standard-resolution data, corresponding to a short cell-identification imaging phase. Finally we demonstrate that interleaved imaging that translates the pixels by subpixel shifts on each frame further improves the fidelity of the recovered neural time series.
Constrained non-negative matrix factorization (see Methods) relies on the observation that the spatiotemporal fluorescence activity (represented as a space-by-time matrix) can be expressed in terms of a product of two matrices: a spatial matrix A that encodes the location and shape of each neuron and a temporal matrix C that characterizes the calcium concentration within each neuron over time. Placing constraints on the spatial footprint of each neuron (e. g. , enforcing sparsity and locality of each neural shape) and on the temporal activity (modeling the observed calcium in terms of a filtered version of sparse, non-negative neural activity) significantly improves the estimation of these components compared to vanilla NMF [8]. Below (cf. Fig 1), we describe a number of algorithmic improvements on the basic approach described in [8]: an iterative block-coordinate descent algorithm in which we optimize for components of A with C held fixed, then for C with A held fixed. We begin by considering imaging data obtained at low temporal resolution, specifically a whole-brain light-sheet imaging recording acquired at a rate of 2 Hz using nuclear localized GCaMP6f in zebrafish. We restricted our analysis to a representative patch shown in Fig 2A, extracted from a medial z-layer of the telencephalon (pallium). (Similar analyses were also performed on patches from midbrain and hindbrain, with similar conclusions.) The neural centers were detected automatically using the greedy method from [8]. To ensure that the spatial components in A are localized, we constrained them to lie within spatial sub-patches (dashed squares in Fig 2A; see also Methods). The first algorithmic improvement follows from the realization that some of the constraints applied in CNMF are unnecessary, at least during early iterations of the algorithm, when only crude estimates for A and C are available. Specifically, [8] imposed temporal constraints on C in each iteration: namely, C was modeled as a filtered version of a nonnegative neural activity signal S—i. e. , CG = S, for an invertible matrix G—and therefore CG is constrained to be non-negative. We found that enforcing a simpler non-negativity constraint on C instead of CG (and then switching to impose the constraint on CG only once the estimates of A and C were closer to convergence) led to a simpler algorithm enabling faster early iterations with no loss in accuracy. Next we found that significant additional speed-ups in this simplified problem could be obtained by simply changing the order in which the variables in this simplified block-coordinate descent scheme are updated [12]. Instead of updating the temporal activity and spatial shape of one neuron at a time (Fig 1A, black line) as in [8], which is known as hierarchical alternating least squares (HALS, [9]) or rank-one residue iteration (RRI, [14]), it turned out to be beneficial to update the activities of all neurons while keeping their shapes fixed, and then updating all shapes while keeping their activities fixed (Fig 1A, vermilion line). The ensuing method is a constrained version of the fast hierarchical alternating least squares (fast HALS, [15]) for NMF; one major advantage of this update ordering is that in each iteration we operate on smaller matrices obtained as dot-products of the data matrix Y with A or C, and there is no need to compute the large residual matrix Y − AC (which is of the same size as the original video) [10,11]. (In the comparisons below we computed the residual to quantify performance, but excluded the substantial time spent on its computation from the reported wall time values.) Next we reasoned that to obtain a good preliminary estimate of the spatial shape matrix A, it is likely unnecessary to use the original data at full temporal resolution [13]. Thus we experimented with the following approach: downsample temporally by a factor of k, then run constrained fast HALS (as described above) for 30 iterations, and then finally return to the original (non-downsampled) data and run a few more iterations of fast HALS until convergence. We experimented with three different downsampling methods: 1) selection of the k-th frame (this could be considered a kind of stochastic update rule, since we are forming updates based only on a subset of the data); 2) forming a median over the data in each block of k frames (applying the median over each pixel independently); and 3) forming a mean over each block of k frames. The mean approach (3) led to significantly more accurate and stable results than did the subsampling approach (1), consistent with the results of [12], and was about an order of magnitude faster than the median approach (2) with similar accuracy, so we restrict our attention to the mean approach (3) for the remainder of this work. (A further advantage of approach (3) relative to (1) is that (1) can miss fast activity transients.) Fig 1B shows the results obtained for a varying number of decimation factors k; we conclude that temporal decimation provides another significant speedup over the results shown in Fig 1A. The starting point of each line is the MSE for keeping the neural shapes obtained on decimated data fixed and solving for the time series on the full data. Refining the shapes on the full data further decreases the MSE, however by less then 1%, hence good shapes are obtained even using merely decimated data. Besides downsampling methods, we also considered dimensionality reduction via structured random projections [16] or singular value decomposition (SVD, [17], see Methods). We compressed the data by the same factor k = 30, but found that both of these dimensionality-reduction methods were less efficient than simple decimation (Fig 1C). We further evaluated whether we can gain improvements by further compressing the decimated data via SVD or random projections, such that the reduced dimension is just slightly larger than the number of neurons. However, we did not obtain any improvements beyond plain decimation. We found this result to hold also for smaller patches that contained fewer neurons. Further speed gains were obtained when applying spatial decimation (computing a mean within l × l pixel blocks) in addition to temporal decimation over the 30 preliminary fast HALS iterations (Fig 1C); see Algorithm 1 for full details. Strikingly, spatial decimation led not only to faster but also to better solutions of the biconvex factorization problem (where solution quality is measured by the residual sum of square errors, ||Y − AC||2, RSS), apparently because the spatially-decimated solutions are near better local optima in the squared-error objective function than are the non-decimated solutions. In summary, by simplifying the early iterations of the CNMF algorithm (by removing the temporal deconvolution constraints to use fast HALS iterations on temporally and spatially subsampled data), we obtained remarkable speed-ups without compromising the accuracy of the obtained solution, at least in terms of the sum-of-squares objective function. But how do these modifications affect the extracted neural shapes and activity traces? We ran the algorithm without decimation until convergence and with decimation for 1 and 10 s respectively. Fig 2 shows the results for three neurons with overlapping patches. Both shapes and activity traces agree well even if the decimated algorithm is run for merely 1 s (Fig 2B and 2C) and are nearly identical if run longer (Fig 2D); hence, decimation does not impair the final obtained accuracy. Our focus has been on speeding up CNMF, one computational bottleneck of the entire processing pipeline. For completeness, we report the times spent on each step of the pipeline in Table 1 and compare to the previous CNMF version of [8]. After loading, the data was decimated temporally by a factor of 30 to speed up the detection of the neural centers using the greedy initialization method from [8]. We further decimated spatially and ran fast HALS for 30 iterations before finally returning to the whole data and performing five final fast HALS iterations. Each trace was normalized by the fluorescence at the resting state (known as ΔF/F) to account for baseline drift using a running percentile filter. Finally, the fluorescence traces were denoised via sparse non-negative deconvolution, using the recently developed fast method of [18], which eliminated another computational bottleneck present in the original CNMF implementation (last row of Table 1). We have shown that decimation leads to much faster computational processing of calcium video data. More importantly, these results inspired us next to propose a method for faster image acquisition or for imaging larger neural populations. The basic idea is quite simple: if we can estimate the quantities of interest (A and C) well given decimated data, then why collect data at the full resolution at all? Since spatial decimation by a factor of l conceptually reduces the number of pixels recorded over a given FOV by a factor of l2 (though of course this situation is slightly more complex in the case of scanning two-photon imaging; we will come back to this issue below), we should be able to use our newly-expanded pixel budget to image more cells, or image the same population of cells faster. As we will see below, this basic idea can be improved upon significantly: if we have a good estimate for the spatial neural shape matrix A at the original spatial resolution, then we can decimate more drastically (thus increasing this l2 factor) with minimal loss in accuracy of the estimated activity C. This, finally, leads to the major proposal of this paper: first perform imaging with standard spatial resolution via conventional imaging protocols. Next perform the ROI detection and CNMF described above to obtain a good estimate of A. Then begin acquiring spatially l-decimated images and use the C-estimation step of CNMF to extract and demix the imaged activity. As we will see below, this two-phase imaging approach can potentially enable the accurate extraction of demixed neural activity even given quite large decimation factors l, with a correspondingly large increase in the resulting “imaging budget. ”
The basic message of this paper is that standard approaches for imaging calcium responses in large neuronal population—which have historically been optimized so that humans can clearly see cells blink in the resulting video—lead to highly redundant data, and we can exploit this redundancy in several ways. In the first part of the paper, we saw that we can decimate standard calcium imaging video data drastically, to obtain order-of-magnitude speedups in processing time with no loss (and in some cases even some gain) in accuracy of the recovered signals. In the second part of the paper, we saw that, once the cell shapes and locations are identified, we can drastically reduce the spatial resolution of the recording (losing the ability to cleanly identify cells by eye in the resulting heavily-pixelated movies) but still faithfully recover the neural activity of interest. This in turn leads naturally to a proposed two-phase imaging approach (first, identify cell shapes and locations at standard resolution; then image at much lower spatial resolution) that can be seen as an effort to reduce the redundancy of the resulting video data. We anticipate a number of applications of the results presented here. Regarding the first part of the paper: faster computational processing times are always welcome, of course, but more fundamentally, the faster algorithms developed here open the door towards guided experimental design, in which experimenters can obtain images, process the data quickly, and immediately use this to guide the next experiment. With more effort this closed-loop approach can potentially be implemented in real-time, whether for improving optical brain-machine interfaces [26], or enabling closed-loop optogenetic control of neuronal population dynamics [27,28]. Highly redundant data streams are by definition highly compressible. The results shown in S1 and S2 Videos illustrate clearly that spatially-decimated image acquisition (the second phase of our two-phase imaging approach) can be seen as a computationally trivial low-loss compression scheme. Again, regarding applications of this compression viewpoint: reductions in memory usage are always welcome—but more fundamentally, this type of compression could for example help enable wireless applications in which bandwidth and power-budget limitations are currently a significant bottleneck [21,29–31]. Regarding applications of the proposed two-phase imaging approach: we can potentially use this approach to image either more cells, or image cells faster, or some combination of both. In most of the paper we have emphasized the first case, in which we ‘zoom out’ to image larger populations at standard temporal resolution. However, a number of applications require higher temporal resolution. One exciting example is the larval zebrafish, where it is already possible to image the whole brain, but light-sheet whole-brain volumetric imaging rates are low [1] and current efforts are focused on faster acquisition [3,4, 32]. Higher temporal resolution is also needed for circuit connectivity inference [33,34] or the real-time closed-loop applications discussed above, where we need to detect changes in activity as quickly as possible. Finally, genetically encoded voltage indicators [35] may soon enable imaging of neuronal populations with single-cell, millisecond-scale resolution; these indicators are still undergoing intense development [36–39] but when more mature the resulting signals will be much faster than currently-employed calcium indicators and significantly higher temporal resolution will be required to capture these signals. A number of previous papers can be interpreted in terms of reducing the redundancy of the output image data. Our work can be seen as one example of the general theme of increasing the ratio N/D, with N denoting the number of imaged neurons and D the number of observations per timestep, with demixing algorithms used post hoc to separate the overlapping contributions of each cell to each observed pixel. In a compressed sensing framework, [40] proposed to image randomized projections of the spatial calcium concentration at each timestep, instead of measuring the concentration at individual locations. In [41], [42], and [43], information is integrated primarily across depth, either by creating multiple foci, axially extended point spread functions (PSFs), or both, respectively. In contrast to these methods, [7] instead scanned an enlarged near-isotropic PSF, generated with temporal focusing, to quickly interrogate cells in a single plane at low spatial resolution. This approach is closest in spirit to the one-phase spatially decimated imaging approach analyzed in Figs 6–8, and could potentially be combined with our two-phase approach to achieve further speed/accuracy gains. We expect that different strategies for increasing N/D will have different advantages in different situations. One advantage of the approach developed here is its apparent simplicity—at least at a conceptual level, we just need to ‘zoom out’ without the need for radically new imaging hardware. Throughout this work we have remained deliberately agnostic regarding the physical implementation of the spatial decimation; all of the decimation results presented here were based on software decimation after acquisition of standard-resolution images. Thus to close we turn now to a discussion of potential experimental caveats. One critical assumption in our simulations is that the total recorded photon flux per frame is the same for each decimation level l. This is a reasonable assumption for light-sheet imaging (assuming we are not limited by laser power or by the peak or average light power on the sample): in this case, increasing the effective pixel size could be achieved easily, either isotropically with a telescope, or anisotropically, with a cylindrical lens or anamorphic prism pair. However, faster whole-brain light-sheet imaging requires faster shifts of the light sheet and imaged focal plane. This challenge can be solved by extended depth-of-field (EDoF) pupil encoding [3,4, 32], remote focusing [44], or with an electrically tunable lens [45]. Higher light-sheet imaging rates can also be obtained with swept confocally-aligned planar excitation (SCAPE) microscopy [6]. In short, we believe our proposed two-phase imaging approach fits well with a variety of proven light sheet methods; for similar reasons, the two-phase approach would also fit well with light-field imaging methods [46–48]. In traditional two-photon imaging the situation is more complicated. The image is created by serially sweeping a small, diffraction limited point across the sample. Along the “fast” axis, the beam moves continuously, and the integrated signal across a line is constant, regardless of detection pixelation—the signal is simply partitioned into more or fewer bins. Along the “slow” axis, however, the galvonometers are moved in discrete steps, and low pixel numbers generally mean that portions of the image are not scanned, increasing frame speed, but concomitantly these ‘missed’ areas generate no signal. This consequently reduces the total number of photons collected. Thus to achieve the same photon flux over the larger (lower spatially sampled) pixels, while maintaining the same SNR, we require an enlarged PSF, which maps a larger sampled volume to each pixel. This approach was recently demonstrated to be effective in [7]; alternative strategies for enlarging the PSF could involve fixed diffractive optical elements [49] or spatial light modulator (SLM) systems [50]. Programmable phase-only SLMs offer the additional benefit of being able to dynamically change the size and shape of the excitation PSF, even between frames, which may help disambiguate closely spaced sources, and effectively control the recorded source sparsity. In any instantiation, maximal imaging speed will be limited by the time required to collect enough photons for adequate SNR, which in turn is limited by photophysics and the light tolerance of the sample. In future work we plan to pursue both light-sheet and 2P implementations of the proposed two-phase imaging approach, to quantify the gains in speed and FOV size that can be realized in practice. We also expect techniques for denoising, demixing, and deconvolution of calcium imaging video to continue to improve in the near future, as more accurate nonlinear, non-Gaussian models for calcium signals and noise are developed; as new demixing methods become available, we can easily swap these methods in in place of the CNMF approach used here. We expect that the basic points about temporal and spatial decimation discussed in this paper will remain valid even as newer and better demixing algorithms become available.
Light-sheet imaging of zebrafish was conducted according to protocols approved by the Institutional Animal Care and Use Committee of the Howard Hughes Medical Institute, Janelia Research Campus. Two-photon imaging of mouse was carried out in accordance with animal protocols approved by the Columbia University Institutional Animal Care and Use Committee. The calcium fluorescence of the whole brain of a larval zebrafish was recorded using light-sheet imaging. It was a transgenic (GCaMP6f) zebrafish embedded in agarose but with the agarose around the tail removed. The fish was in a fictive swimming virtual environment as described in [51]. The closed loop setting, characterized by visual feedback being aligned with the recorded motor activity, was periodically interrupted by open loop phases. Whole-brain activity was recorded for 1,500 seconds with a rate of 2 volumes per second. In vivo two-photon imaging was performed in a transgenic (GCaMP6s) mouse through a cranial window in visual cortex. The mouse was anesthetized (isoflurane) and head-fixed on a Bruker Ultima in vivo microscope with resonant scanners, and spontaneous activity was recorded. The field of view extended over 350 μm × 350 μm and was recorded for 100 seconds with a resolution of 512×512 pixels at 20 frames per second. In the case of 2P imaging, the field of view contained N = 187 ROIs. It was observed for a total number of T = 2,000 timesteps and had a total number of D = 512×512 pixels. We restricted our analysis of the zebrafish data to a representative patch of size D = 96×96 pixels containing N = 46 ROIs, extracted from a medial z-layer of the whole-brain light-sheet imaging recording of T = 3,000 frames. The observations at any point in time can be vectorized in a single column vector of length D; thus all the observations can be described by a D × T matrix Y. Following [8], we model Y as Y = A C + b f ⊤ + E (1) where A ∈ R + D × N is a spatial matrix that encodes the location and shape of each neuron, C ∈ R + N × T is a temporal matrix that characterizes the calcium concentration of each neuron over time, b ∈ R + D, f ∈ R T are nonnegative vectors encoding the background spatial structure and global intensity, respectively, and E is additive Gaussian noise with mean zero and diagonal covariance. Our model assumes a linear relationship between fluorescence and calcium concentration as well as Gaussian noise. As emphasized in [8], more elaborate models for C can be incorporated in the alternating updates, but we did not pursue this generalization here. For the zebrafish data we ensured that the spatial components are localized, by constraining them to lie within spatial patches (which are not large compared to the size of the cell body) around the neuron centers, thus imposing sparsity on A by construction. Because of the low temporal resolution of these recordings, the inferred neural activity vectors are not expected to be particularly sparse, and therefore we do not impose sparsity in the temporal domain. This leads to the optimization problem minimize A, b, C, f ∥ Y - A C - b f ⊤ ∥ subject to: A, b, C ≥ 0, A (d, n) = 0 ∀ d ∉ P n (2) where Pn denotes the n-th fixed spatial patch. This problem is biconvex, i. e. solving for C and f with fixed A and b is a convex problem; likewise solving for A and b with fixed C and f is convex. As discussed in the Results section, we solve this problem by block-coordinate descent, first applied to much smaller decimated data and then using this solution as a warm start for the optimization on the full data. In the resulting Algorithm 1 we appended for concision b and f as an additional column or row to A and C respectively. The matrix products A⊤Y and CY⊤ in Algorithm 1 are computationally expensive for the full data. These matrix products can also be performed on GPU instead CPU; whereas for the comparatively small 96×96 patches we did not obtain any speed-ups using a GPU, we verified on patches of size 256×256 that some modest overall speedups (a factor of 1. 5–2) can be obtained by porting this step to a GPU. For the decimated data the matrix products are cheap enough to iterate just once over all neurons and instead alternate more often between updating shapes and activities (instead of performing many iterations within HALSactivity or HALSshape in Algorithm 1). In early iterations our estimates of A and C are changing significantly and it is better to perform just one block-coordinate descent step for each neuron to update A (and similarly for C); for later iterations, and on the full data where it is more expensive to compute A⊤Y and CY⊤, we increase the inner iterations in HALSactivity or HALSshape. Algorithm 1 is a constrained version of fast HALS. To further improve on fast HALS, [52] suggested to replace cyclic variable updates with a greedy selection scheme focusing on nonzero elements. This was unnecessary here because most nonzero elements are prespecified by the patches Pn; i. e. , we are already focusing on the nonzero elements. Because for the 2P data the observed imaging rate is much higher than the decay rate of the calcium indicator, we constrain the temporal traces C to obey the calcium indicator dynamics, to enable further denoising and deconvolution of the data. As in [8], we approximate the calcium concentration dynamics using an autoregressive process of order 2 (AR (2) ), C (n, t) = γ 1 C (n, t - 1) + γ 2 C (n, t - 2) + S (n, t), (3) where S (n, t) is the number of spikes that neuron n fired at timestep t. This equation can be conveniently expressed in matrix form as S = CG for a suitable sparse matrix G. We estimate the noise level of each pixel σd by averaging the power spectral density (PSD) over a range of high frequencies, and estimate the coefficients of the AR (2) process for each cell following [22]. Then we solve for A, b, C, f using the following iterative matrix updates: minimize A, b ∥ A ∥ 1, subject to: A, b ≥ 0, ∥ Y (d, :) - A (d, :) C - b (d) f ⊤ ∥ ≤ σ d T ∀ d ∈ { 1,2, . . . , D } (4) minimize C, f ∥ C G ∥ 1, subject to: C G ≥ 0, ∥ Y (d, :) - A (d, :) C - b (d) f ⊤ ∥ ≤ σ d T ∀ d ∈ { 1,2, . . . , D }. (5) These updates are initialized with the results from constrained fast HALS (Alg 1). They impose sparsity on the spatial as well as temporal components, using the estimates of the noise variance as hard constraints to derive a parameter-free convex program. Following the approach in [53] the spike signal S is relaxed from nonnegative integers to arbitrary nonnegative values. The basis pursuit denoising problems in Eqs (4 and 5) can be solved with one of the methods described in [8]. However, a faster update of the temporal matrix C is achieved by using OASIS [18]. Every spatial component in A was normalized to have unit ℓ2-norm, with the corresponding temporal component scaled accordingly. Following [8] we then sort the components according to the product of the maximum value that the temporal component attains and the ℓ4-norm of the corresponding spatial footprints, to penalize overly broad and/or noisy spatial shapes. In order to calculate the ΔF/F values we divided the fluorescence trace C (n, :) of each neuron by its baseline fluorescence that was obtained by projecting the rank-1 background bf⊤ onto the shape A (: , n) of the neuron, (Δ F F) t ≔ A (: , n) ⊤A (: , n) C (n, t) A (: , n) ⊤ b f (t). While we normalized A and C such that A (: , n) ⊤A (: , n) = 1, the ΔF/F values do not depend on this specific normalization because the scaling factor cancels out. To compress the data using truncated SVD in Fig 1C, we followed [17] and computed the eigenvectors V ∈ R M × T belonging to the M largest eigenvalues of the time by time covariance matrix Y⊤Y, which were then used to obtain the compressed data YV⊤. This method was faster than the randomized method due to [54]. The spatial background b for the compressed data was again initialized as 20% percentile of the original data and the temporal background as f ˜ = V 1 T. Because the compressed data and temporal traces can be negative we did not enforce the non-negativity constraint of C ˜ and f ˜, which was crucial as enforcing it indeed yielded worse results. Random compression was performed as described in [16]. We applied the structured random compression algorithm [54] to Y and Y⊤ to obtain L ∈ R D × M and R ∈ R M × T. Specifically, we drew a Gaussian random matrix Ω ∈ R T × M and performed the QR decomposition of YΩ to obtain an orthonormal basis L. Analogously we obtained R for Y⊤. The iterated alternating fast HALS updates were C ← HALSactivity (L⊤Y, L⊤A, C) and A ← HALSshape (YR⊤, A, CR⊤), with L⊤Y and YR⊤ computed once initially. Applying the code of [55] to the raw data we identified the neural shape matrix A1 and spatial background b1. We use the convention that the presence of a lower index l signifies an estimate and its value the decimation factor, i. e. index l = 1 denotes an estimate inferred without decimation. Further, we also obtained the denoised and deconvolved traces C1, S1 as well as f1. To emulate imaging with lower spatial resolution, spatial decimation was performed by converting A1 back into a 512 × 512 × N tensor (Y into a 512 × 512 × T video tensor) and calculating the average of non-overlapping patches of size l × l or l × 1 pixels for each of the N neural shapes (T timesteps). Converting the tensors back to matrices yielded the decimated neural shapes Al (data Yl). We proceeded analogously for the spatial background to obtain bl. The corresponding temporal traces were estimated by solving Eq (5) (with Yl replacing Y, Cl replacing C, etc.), initializing Cl and fl with the result of plain NMF that does not impose temporal constraints, i. e. solving minimizeC l, f l ∥ Y l - A l C l - b l f l ⊤ ∥ subject to Cl ≥ 0. In order to obtain the results for 1-phase imaging without previous shape identification we solved Eq (4) for the decimated data Yl, initializing Al, bl by decimating A1, b1 and setting the temporal components to C1, f1. With increasing decimation factors an increasing number of shapes got purged and absorbed in the background, reflecting the fact that it would have been difficult to detect all ROIs on low resolution data in the first place. Using the obtained remaining shapes we again solved Eq (5) as above. The correlation values for purged neurons were set to zeros for the mean values reported in Fig 6. To obtain some form of ground truth (Figs 6D–6F, 7C and 8C) we generated a simulated dataset Ys by taking the inferred quantities as actual ground truth: As ≔ A1, bs ≔ b1, Cs ≔ C1, fs ≔ f1. We calculated the residual Y − AsCs − bsfs⊤ and reshuffled it randomly but signal dependent for each pixel in time. We partitioned the residual for each pixel into 200 strata according to signal size and reshuffled it within each strata, thus retaining any potential link between noise variance and signal mean. The simulated dataset Ys was obtained by adding the reshuffled residual to AsCs + bsfs⊤ and the same analysis as for the original data was performed. We performed additional control simulations that also took the inferred quantities as actual ground truth, Y *: = A 1 C 1 + b 1 f 1 ⊤, but did not rely on a reshuffling procedure. Instead, we either added Gaussian noise, Y d, t N ∼ N (Y d, t *, σ 2), or Poisson noise ν + Y d, t P ∼ P (ν + Y d, t *), where ν is the photon count that is not due to calcium fluorescence (but rather dark counts and background light). Whereas the Gaussian noise had fixed variance σ2, the Poisson model results in heteroscedastic noise because its variance grows with the mean. The variance of the Gaussian noise was chosen to be equal to the average variance of the Poisson noise. The results shown in S2 Fig agree with those obtained by reshuffling the residual and are similar for Gaussian and Poisson noise, at least on average, though there is some spread if individual traces are considered (S2C Fig). Whereas the model (Eq 1) assumes Gaussian noise, it nevertheless performs well under Poisson noise, consistent with the results of [53]. Another control simulation merely took A1, b1 and f1 as ground truth. However, instead of taking the denoised fluorescence traces C1, which by construction followed the autoregressive model, the fluorescence traces C D were obtained from two datasets that combined electrophysiological recording and calcium imaging with GCaMP6f (11 cells) or GCaMP6s (9 cells) [23]. The calcium response kernel k ^ for each recorded cell was determined by solving the linear regression problem k ^ = arg min k ∥ y - s * k ∥ 2 where y is the noisy fluorescence data of a neuron, s its spike counts per bin, and * denotes convolution. Because only few cells were recorded, but for a longer duration than the 100 s of our two-photon dataset, we assigned to each of our ROIs the ground truth fluorescence trace by randomly selecting a recorded cell, taking a 100 s interval of its spike train that included at least three spikes, and convolving it with the cell’s kernel k. The trace was scaled such that its maximum value was equal to the maximal value of the inferred trace C1. Poisson noise was added, ν + Y D ∼ P (ν + A 1 C D + b 1 f 1 ⊤), as described above. Results are shown in S3 Fig. Projecting the noise of each pixel onto the neural shapes yields the noise of each neural time series. In practice the latter is estimated based on the noisy trace obtained by projecting the fluorescence data onto the shapes. For interleaved imaging (Fig 8) the shape of each neuron differs between alternating frames due to the varying pixelization. Therefore, instead of using one noise level for all timesteps, we estimated two noise levels σodd and σeven based on the PSD for all odd and even frames respectively. The residuals in the noise constraint of the non-negative deconvolution were weighted accordingly by the inverse of the noise level. where y is the noisy fluorescence data of a neuron (cell index suppressed) obtained by subtracting the contribution of all other neurons as well as the background from the spatio-temporal raw pixel fluorescence data and projecting the odd and even frames of the result onto the considered neuron’s shapes aodd and aeven respectively. yodd and yeven denote the vectors obtained by taking only every second component of y starting with the first/second respectively. The denoised fluorescence c is denoted analogously. For simplicity we estimated the coefficients of the AR (2) process based on all frames without separating by noise level. All analyses were performed on a MacBook Pro with Intel Core i5-5257U 2. 7 GHz CPU and 16 GB RAM. We wrote custom Python scripts that called the Python implementation [55] of CNMF [8]. Our scientific Python installation included Intel Math Kernel Library (MKL) Optimizations for improved performance of vectorized math routines.
The voxel rate of imaging systems ultimately sets the limit on the speed of data acquisition. These limits often mean that only a small fraction of the activity of large neuronal populations can be observed at high spatio-temporal resolution. For imaging of very large populations with single cell resolution, temporal resolution is typically sacrificed. Here we propose a multi-scale approach to achieve single cell precision using fast imaging at reduced spatial resolution. In the first phase the spatial location and shape of each neuron is obtained at standard spatial resolution; in the second phase imaging is performed at much lower spatial resolution. We show that we can apply a demixing algorithm to accurately recover each neuron’s activity from the low-resolution data by exploiting the high-resolution cellular maps estimated in the first imaging phase. Thus by decreasing the spatial resolution in the second phase, we can compress the video data significantly, and potentially acquire images over an order-of-magnitude larger area, or image at significantly higher temporal resolution, with minimal loss in accuracy of the recovered neuronal activity. We evaluate this approach on real data from light-sheet and 2-photon calcium imaging.
fluorescence imaging
statistical noise
applied mathematics
neuroscience
simulation and modeling
algorithms
mathematics
statistics (mathematics)
computational neuroscience
neuroimaging
research and analysis methods
gaussian noise
imaging techniques
animal cells
calcium imaging
cellular neuroscience
cell biology
single neuron function
neurons
biology and life sciences
cellular types
physical sciences
computational biology
Multi-scale approaches for high-speed imaging and analysis of large neural populations