jovo / smc-oopsi

code to implement our sequential monte carlo (smc) algorithm to perform optimal optical spike inference (oopsi) for individual neurons

This is a repository containing the code to implement the algorithms employed in Vogelstein et al, 2009, and subsequent related advances.
Publication version of manuscript is now available from [Cell Press] 4 with subscription, or [here] 5 for free.

Note that the code in this repository is not kept up-to-date; rather, it is an archive of our work.
for an up-to-date code base, check out: https://github.com/jovo/oopsi

This work may be cited as:

Vogelstein, Joshua T.; Watson, Brendon O.; Packer, Adam M.; Yuste, Rafael; Jedynak, Bruno; Paninski, Liam. Spike Inference from Calcium Imaging Using Sequential Monte Carlo Methods. Biophysical Journal, Volume 97, Issue 2, pp.636 - 655, 2009.

Repository Organization

• code: contains current stable version of code
• manuscript: contains everything related to [Vogelstein et al. 2009] 5, including ** code: version of code used for manuscript ** data: both raw data and data generated by simulations for figs ** figs: pdf's of figs in manuscript ** reviews&proofs: reviews and proofs from BJ to get this paper published ** tex: contains tex and compiled other files of submitted manuscript (and some minor revisions)
• general_background: a document containing some general background information that might be useful

Quick Tips

• final version of proof may be found [here] 1
• to learn how to use the algorithm, i suggest playing with scripts in the [script folder] 3
• any questions or issues, please use the [issues tab] 2
• any problems with that, or questions inappropriate for that forum, please contact me at: [email protected]
• please respect our policy of openness, and operate accordingly.

License

License stuff (apparently, this is useful for some people): Copyright (c) 2008 Joshua Vogelstein

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

smc-oopsi's Issues

Difference between smc_oopsi.m and ParticleFilterD.m, Get_Spt.m missing

Hi,
I have two minor comments:
Minor note 1: only nbar is computed in smc_oopsi but not the other additional Cbar, etc that are computed in ParticleFilterD.

Minor note 2: Get_Spt.m is missing from smc_oopsi repository, people can of course pull it from the other repo but maybe it would be good to post it here too.

holdTau does not work properly

when data quality is poor, the algorithm attempts to fit the noise with spikes. to do this, it pushed tau_c --> 0. many options are available to counteract this effect:

1. establish a strong prior on tau_c
2. fix tau_c

in the current version of the code, we chose to fix tau_c. technically, this could be considered a very strong prior, with all the mass at a single point. however, the way we do it is incorrect. we jointly solve for {A,tau_c,C_b}, and then simply let tau_c equal the fixed value. we should compute the sufficient stats for estimating {A,C_b}, and then maximize with respect to only them.

slow rise time of fluorescence

while organic sensors have effectively instantaneous rise times (around 1 ms), genetic sensors have rise times on the order of 10's of msec. thus, even when imaging quickly, when using genetic sensors, the rise might span multiple frames. perhaps the easiest way to incorporate this feature would be to introduce multiple time constants for calcium.

estimating noise terms is not done correctly

our noise model for the observations if:

F_t \sim Normal[ a S(C_t) + b, g S(C_t) + z ]

the proper way to estimate these parameters is to iterate:

1. estimate {a,b} while holding {g,z} fixed
2. estimate {g,z} while holding {a,b} fixed

until convergence (this is an example of coordinate descent). currently, we do something else that is not technically correct. more specifically, we estimate variance, and then we scale {g,z} by the appropriate factor.

multiple spikes per frame

obviously, this is a desiderata of many. even when frame rates approach 30 Hz and 60 Hz, multiple spikes are not impossible. there are two obvious ways to achieve this:

1. superresolution: allowing the number of time steps between observations to increase to the number of spikes observed is a distinct possibility. since computations are O(T), where T is the number of time steps, this makes things substantially slower. to solve that issue, we could introduce two stages. stage 1: no superresolution. stage 2: yes superresolution. stage 2 begins only after stage 1 is complete. further, we could only do superresolution for image frames in which the probability of spiking is high. in fact, the code is already set up to be able to choose for each frame whether to do superresolution, i think.

2. change the distribution from which spikes are sampled to Poisson. this requires a bit more work, given that everything done so far assumes bernoulli. its not clear that this option would be better in any way, although it might be more natural

M step of EM sometimes shows likelihood decreasing

this is a weird problem, given that the M step for EM is guaranteed to be log-concave, and should therefore monotonically increase. not sure what is going on here. my guess is that we don't compute the likelihood properly, because the results do seem to get better, but i can't be sure.

smc_oopsi is not "pretty"

smc_oopsi should merely call other functions, including

smc_oopsi_forward
smc_oopsi_backward
smc_oopsi_m_step

all the code associated with each of those functions should be embedded within those functions.

calibrating parameters using the true (or some estimated) spike train

we'd like to have a few lines of code that estimates the parameters given the true spike train (or some estimate of the true spike train). let 'n' be the guess for what the spike train is. the correct thing to do would be to skip the forward/backward, convolve 'n' with an exponential with time constant tau_c, yielding 'C', and then set the Mstep code to have only a single particle with hidden states {n,C}, and learn parameters from that. currently, we use a weird hack in GOOPSI_main_v1_0 after doing the forward/backward, that makes each particle be have 'n' (and i'm not even really sure where 'C' comes from.