FastLZeroSpikeInference Tutorial

R tutorial

In this tutorial, we demonstrate how the FastLZeroSpikeInference package can be used to deconvolve calcium imaging data from a single neuron. All demos require the FastLZeroSpikeInference package; installation instructions are provided here.

To illustrate the software, we generate a synthetic dataset according to

$y_t = c_t + \epsilon_t, \\ c_t = \gamma c_{t-1} + z_t, \\ z_t \sim \mathrm{Pois}(\theta), \\$

where $\epsilon_t \sim \mathrm{N}(0, \sigma^2)$ and $c_1 = 1$. The function simulate_ar1 generates data from this model.

sim <- simulate_ar1(n = 10000, gam = 0.95, poisMean = 0.01, sd = 0.15, seed = 1)
plot(sim)

plot of chunk unnamed-chunk-1

The estimate_spikes function estimates spikes based on the calcium trace, exponential decay parameter $\gamma$, and a tuning parameter $\lambda$. This function solves

$\underset{c_1,\ldots,c_T, z_2,\ldots,z_T}{\mathrm{minimize}} \frac{1}{2} \sum_{t=1}^T \left( y_t - c_t \right)^2 + \lambda \sum_{t=2}^T 1_{\left( z_t \neq 0 \right) }$

if the constraint parameter is set to false (default), and

$\underset{c_1,\ldots,c_T, z_2,\ldots,z_T}{\mathrm{minimize}} \frac{1}{2} \sum_{t=1}^T \left( y_t - c_t \right)^2 + \lambda \sum_{t=2}^T 1_{\left( z_t \neq 0 \right) } \mbox{ subject to } z_t = c_t - \gamma c_{t-1} \geq 0,$

if the constraint parameter is set to true. Here’s an example of the function applied to the previously generated data.

# AR(1) model
fit <- estimate_spikes(dat = sim$fl, gam = 0.95, lambda = 1)
print(fit)

## 
##  Call: 
## estimate_spikes(dat = sim$fl, gam = 0.95, lambda = 1)
## 
##  Output: 
## Number of spikes 	 89 
## 
##  Settings: 
## Data length 		 10000 
## Model type 		 ar1 
## Gamma 			 0.95 
## Lambda 			 1

By default, and to save computation time, the calcium concentration is not automatically estimated. However, the estimate_calcium function can be used to estimate the calcium concentration based on a prior fit.

# compute fitted values from prev. fit
fit <- estimate_calcium(fit)
plot(fit)

plot of chunk unnamed-chunk-3

Alternatively, both spikes and calcium concentrations can be computed with the estimate_spikes function with estimate_calcium parameter set to true.

fit <- estimate_spikes(dat = sim$fl, gam = 0.95, lambda = 1, estimate_calcium = T)
plot(fit)

plot of chunk unnamed-chunk-4

To solve the constrained problem

we set the constraint parameter to true in the estimate_spikes function. Consider the following example.

# Constrained AR(1) model
fit <- estimate_spikes(dat = sim$fl, gam = 0.95, lambda = 1, constraint = T, estimate_calcium = T)
print(fit)

## 
##  Call: 
## estimate_spikes(dat = sim$fl, gam = 0.95, lambda = 1, constraint = T, 
##     estimate_calcium = T)
## 
##  Output: 
## Number of spikes 	 89 
## 
##  Settings: 
## Data length 		 10000 
## Model type 		 ar1-pos-constrained 
## Gamma 			 0.95 
## Lambda 			 1

plot(fit)

plot of chunk unnamed-chunk-5