calcium imaging

Processing Volumetric Data

After collecting volumetric data with the Optotune you can now process your data as follows.

First, use sbxsplit() to generate separate data files for each “slice” of the optotune.

For example, if I have a data file gn7_008_001.sbx collected some data with an optotune waveform having a period of 5 the command

>> sbxsplit('gn7_008_001')

will generate a set of 5 files named gn7_008_001_ot_NNN, where NNN.sbx will range from 000 to 004, corresponding to each separate optical plane.

Second, the resulting files can then be aligned and segmented by treating each plane individually.  This will generate the corresponding *.signals and *.segment files.

Finally, you can call:

>> sbxmerge('gn7_008_001')

This will generate gn7_008_001_merged.signals and gn7_008_001_merged.segment.

The signals matrix will have as many rows as frames were present during acquisition while interpolating the missing samples for each plane (that is, when the optotune was sampling from other planes).

The interpolated signals are then deconvolved as usual to generate an estimate of spiking in the spks matrix.  From here on you can process the data as if it came from a typical experiment where only one plane was sampled.

The mask variable in the segmented file will have a size of [ny nx plane], where [ny, nx] is the size of each frame and plane is the period of the optotune waveform.  Each cell has a unique ID value corresponding to its column in the signals matrix.


Note that each setting of the optotune waveform is treated independently, even though thay may potentially represent the same plane (as it may happen using sinusoidal or triangular z-scanning).

An extremely simple deconvolution algorithm

A number of sophisticated deconvolution algorithms have been proposed in the literature, most of which involve solving an L_1 optimization problem.

While trying to develop a fast, on-line spike detection algorithm for Scanbox I tried a simple idea that seems to work reasonably well and the resulting algorithm is trivial.

In all likelihood this must have been published somewhere.  If you happen to know a relevant reference, just let me know in the comments below.

Suppose the calcium signal, y_n, is the response of a linear system with an exponential impulse response to a train of impulses u_n.

Then, the response satisfies y_n = \alpha y_{n-1} + u_n + \epsilon_n, where we assume \epsilon_n is Gaussian noise.

One simple way to estimate u_n is to first find the value of \alpha that minimizes the error  E \{ ( y_n - \alpha y_{n-1} - u_n )^2 \}.

A few calculations yield \alpha = ( \mu^2 - m_{12}  ) /   ( \mu^2 - m_{02}  ), where \mu = E\{ y_n \}, m_{12} = E\{y_n y_{n-1}\} and m_{02} = E\{y_n^2\}.

As all these quantities are easily estimated recursively as new data arrive, one can also maintain a recursive estimate of the optimal \alpha.

Then, the estimated input to the system can be estimated by computing \hat{u}_n = y_n -  \alpha y_{n-1}.  So, all we need to do is to keep a recursive estimate of \alpha and the last two samples of the signal.

Below is an example of how this simple algorithm performs.deconv_easyA simulated spike train, u_n (second panel form top) is convolved with an exponential impulse response and corrupted by additive Gaussian white noise to yield a simulated signal y_n (top panel).  The third panel shows the estimated \hat{u}_n. Finally, we threshold this signal to obtain the estimated spike train shown in the bottom panel (thresholding is done by means of Otsu’s algorithm, implemented by Matlab’s gray thresh).

Here is the code if you want to play with it:

% generate ca-trace

u = double(rand(1,10000)<0.1); % original spikes
h = exp(-(0:40)/10);       % exponential impulse response
y = conv(h,u);             % convolve
y = y+0.15*randn(size(y)); % additive gaussian noise 

% estimate alpha, uhat and input spike train
m = mean(y);
m02 = mean(y.^2);
m12 = mean(y(2:end-1) .* y(1:end-2));
a = (m^2-m12) / (m^2-m02);
uhat = [0 y(2:end) - a * y(1:end-1)];    % uhat
shat = uhat>graythresh(uhat);         % estimated spike train

Here is another example with a higher noise level and some bursting:


Here is how it runs on a some of the cases of the Chen et al dataset that have varying degrees of difficulty.  The red dots represent the spikes detected as described in the original paper, the black dots represents the spikes as detected by the deconvolution algorithm above.

chen_010 chen_009 chen_008 chen_007 chen_006 untitled chen_005 chen_004 chen_003 chen_002 chen_001

Given the simplicity of the algorithm, it performs remarkably well and can be run across many ROIs in real time.

But wait… there is more. One may expect the impulse response to be not purely exponential. In such a case, it may be helpful to generalize the above to a higher-order model. You can then estimate \hat{u}_n in just two lines of Matlab as follows:

a = lpc(y,10); % order=10
uhat = y-filter([0 -a(2:end)],1,y);

The above should bring back memories to anyone who studied speech compression and recognition at some point in time.

The graph below shows what happens when the impulse response is a gamma function (left).  Now, the LPC coefficients are non-zero for lags larger than one.


Still, and simple thresholding of \hat{u}_n still recovers the original spike train even under the presence of substantial noise:

gamma_lpc_resultGo ahead, load some of your calcium imaging data and plot both your signal y_n along with \hat{u}_n estimated by the two lines above. Here is are a couple of data segments, again, from the Chen et al dataset.

uhat1 uhat0

Note: There are other pre-processing steps I took (like projecting out the neuropil signal) that I have not described above.

The calculation:

We want to minimize J = E\{( y_n - \alpha y_{n-1} - u_n )^2\}.  Taking the derivative with respect to \alpha and equating to zero leads to

E\{y_n y_{n-1}\} - \alpha E\{y_{n-1}^2\} - E\{u_n y_{n-1}\} = 0


m_{12}- \alpha  m_{02} - \lambda \mu = 0

where E\{ u_n\} = \lambda and E\{ y_n \} = \mu.

But taking expectations on both sides of y_n = \alpha y_{n-1} + u_n leads to

\mu = \alpha \mu + \lambda or \lambda = \mu (1-\alpha).

Substituting \lambda in the equation above and solving for \alpha gives \alpha = (\mu^2 - m_{12})/(\mu^2 - m_{02}).