Calcium Deconvolution from Fluorescence Signals

We assume we observe the fluorescence signal for T timesteps, and denote by s_t the number of spikes that the neuron fired at the t-th timestep, t = 1, …, T, cf. Fig 1. Following [5 (link), 13 (link)], we approximate the calcium concentration dynamics c using a stable autoregressive process of order p (AR(p)) where p is a small positive integer, usually p = 1 or 2,
$c_t=\sum _{i=1}^p{\gamma }_i{c}_{t-i}+s_t.$
The observed fluorescence $\mathit{y}\in {ℝ}^T$ is related to the calcium concentration as [5 (link)–7 ]:
$y_t=a{c}_t+b+{ϵ}_t,{ϵ}_t\sim \mathcal{N}(0,{\sigma }^2)$
where a is a non-negative scalar, b is a scalar offset parameter, and the noise is assumed to be i.i.d. zero mean Gaussian with variance σ². For the remainder we assume units such that a = 1 without loss of generality. We begin by assuming b = 0 for simplicity, but we will relax this assumption later. (We also assume throughout that all parameters in sight are fixed; in case of e.g. drifting baselines b we could generalize the algorithms discussed here to operate over shorter temporal windows, but we do not pursue this here.) The parameters γ_i and σ can be estimated from the autocovariance function and the power spectral density (PSD) of y respectively [13 (link)]. The autocovariance approach assumes that the spiking signal s comes from a homogeneous Poisson process and in practice often gives a crude estimate of γ_i. We will improve on this below by fitting the AR coefficients directly, which leads to better estimates, particularly when the spikes have some significant autocorrelation.
The goal of calcium deconvolution is to extract an estimate $\hat{\mathit{s}}$ of the neural activity s from the vector of observations y. As discussed in [5 (link), 13 (link)], this leads to the following non-negative LASSO problem for estimating the calcium concentration:
$\underset{\widehat{\mathit{c}},\widehat{\mathit{s}}}{\text{minimize}}\frac{1}{2}\parallel \widehat{\mathit{c}}{-\mathit{y}\parallel }^2+\lambda {\parallel \widehat{\mathit{s}}\parallel }_1\text{subject to}\widehat{\mathit{s}}=G\widehat{\mathit{c}}\ge 0$
where the ℓ₁ penalty on $\hat{\mathit{s}}$ enforces sparsity of the neural activity and the lower triangular matrix G is defined as:
$G=\left(\begin{array}{ccccccc}1& 0& \cdots & \cdots & \cdots & \cdots & 0\\ -{\gamma }_1& 1& \ddots & \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ -{\gamma }_p& \cdots & -{\gamma }_1& 1& 0& \cdots & 0\\ 0& -{\gamma }_p& \cdots & -{\gamma }_1& 1& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& -{\gamma }_p& \cdots & -{\gamma }_1& 1\end{array}\right).$
The deconvolution matrix G is banded with bandwidth p for an AR(p) process. Equivalently, s = c * g with g a finite impulse response filter of order p (p + 1 filter taps) and * denoting convolution. To produce calcium trace c, spike train s is filtered with the inverse filter of g, an infinite impulse response h, c = s * h. (Although our main focus is on the autoregressive model, we will discuss more general convolutional observation models below as well, and touch on nonlinear effects such as saturation in the Appendix.) Following the approach in [5 (link)], note that the spike signal $\hat{\mathit{s}}$ is relaxed from non-negative integers to arbitrary non-negative values.