Fast online deconvolution of calcium imaging data

Fluorescent calcium indicators are a popular means for observing the spiking activity of large neuronal populations, but extracting the activity of each neuron from raw fluorescence calcium imaging data is a nontrivial problem. We present a fast online active set method to solve this sparse non-negative deconvolution problem. Importantly, the algorithm 3progresses through each time series sequentially from beginning to end, thus enabling real-time online estimation of neural activity during the imaging session. Our algorithm is a generalization of the pool adjacent violators algorithm (PAVA) for isotonic regression and inherits its linear-time computational complexity. We gain remarkable increases in processing speed: more than one order of magnitude compared to currently employed state of the art convex solvers relying on interior point methods. Unlike these approaches, our method can exploit warm starts; therefore optimizing model hyperparameters only requires a handful of passes through the data. A minor modification can further improve the quality of activity inference by imposing a constraint on the minimum spike size. The algorithm enables real-time simultaneous deconvolution of O(105) traces of whole-brain larval zebrafish imaging data on a laptop.

S1 Appendix. Technical appendix 1 Algorithm for isotonic regression without pooling 2 For ease of exposition Alg A shows the pseudocode of the isotonic regression algorithm 3 used to convey the core idea. However, this naïve implementation lacks pooling, 4 rendering it inefficient. It repeatedly updates all values x t , ..., x τ during backtracking 5 and calculates the updated value using Eq (7) without exploiting that part of the sum 6 in the numerator has already been computed as an earlier result. It is thus merely 7 O(T 2 ), whereas introducing pools addresses both issues and yields an O(T ) algorithm. while t > 1 and x t < x t −1 do track back 5: For sake of generality we consider the case of weighted regression with weights θ.
This generalization is not only of theoretical interest. These weights could be used to 11 give lower weight to time points with higher variance for heteroscedastic data, for where the nonlinear function f can include saturation effects. This is often taken to be 17 the Hill equation, i.e., f (c) = ac n c n +k d + b, with Hill coefficient n, dissociation constant k d , 18 scaling factor a and baseline b [1]. Applying Newton's algorithm to optimize forŝ (or 19 equivalentlyĉ) results for each Newton step in a weighted constrained regression 20 problem as in Eq (S1), which can be solved efficiently with OASIS. Hence, incorporating 21 OASIS into Newton's algorithm enables the algorithm to handle nonlinear and 22 non-Gaussian measurements.

23
For an AR(1) process introducing weights changes Eq (10) to and its solution is a modification of Eq (11) by adding the weights We merely need to initialize each pool as (v t , w t , t t , l t ) = (y t − µt θt , θ t , t, 1) for each time 26 step t and the updates according to Eqs (12-14) guarantee that Eq (S4) holds for all 27 values v i = c ti as we prove in the next section.

28
For an AR(p) process introducing weights changes Eq (30) to and the same modified initialization holds.
where we used the contingency of the pools, t i+1 = t i + l i . Thus after the update Eq (S6) holds for the merged pool too. It remains to show this also for the values: Initial calcium fluorescence as well known in the AR / linear systems literature [2]. The first pool is merged with 54 the second one whenever the constraint v 2 ≥ d l1 v 1 is violated.

55
Explicit expressions of the hyperparameter updates for AR(2) 56 We solve the noise constrained problem by increasing λ in the dual formulation until the 57 noise constraint is tight. We start with some small λ, e.g. λ = 0, to obtain a first 58 partitioning into pools P. 59 We denote all except the differing last two components of µ by µ = λ(1 − γ 1 − γ 2 ) 60 (Eq 27) and express the components of µ as µ t = µ κ t with else.
Given some µ(λ), the value of the first pool used to model the initial calcium The other values in the pool are implicitly defined by 67 c ti+t = γ 1ĉti+t−1 + γ 2ĉti+t−2 for t = 1, ..., l i − 1. (S13)