Sleep State Analysis Using Calcium
Imaging Data by Non-negative Matrix
Factorization
Mizuo Nagayama1(B) , Toshimitsu Aritake1 , Hideitsu Hino2,3 ,Takeshi Kanda4 , Takehiro Miyazaki4 , Masashi Yanagisawa4 ,Shotaro Akaho3,5 , and Noboru Murata1,31

Waseda University, Shinjuku, Tokyo 169-0072, Japanmizuo.n@asagi.waseda.jp
2
The Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, Japan
3
RIKEN Center for Advanced Intelligence Project, Chuo, Tokyo 103-0027, Japan
4
University of Tsukuba, Tsukuba, Ibaraki 305-8573, Japan
5
National Institute of Advanced Industrial Science and Technology,
Tsukuba, Ibaraki 305-8568, Japan

Abstract. Sleep is an essential process for the survival of animals. However, its phenomenon is poorly understood. To understand the phenomenon of sleep, the analysis should be made from the activities ofa large number of cortical neurons. Calcium imaging is a recently developed technique that can record a large number of neurons simultaneously,
however, it has a disadvantage of low time resolution. In this paper, we
aim to discover phenomena which characterize sleep/wake states from
calcium imaging data. We made an assumption that groups of neurons
become active simultaneously and the neuronal activities of groups diﬀer
between sleep and wake states. We used non-negative matrix factorization (NMF) to identify those groups and their neuronal activities in timefrom calcium imaging data. NMF was used because neural activity canbe expressed by the sum of individual neuronal activity and ﬂuorescenceintensity data are always positive values. We found that there are certain
groups of neurons that behave diﬀerently between sleep and wake states.
Keywords: Sleep state analysis · Calcium imaging ·Non-negative matrix factorization

1

Introduction

Sleep is controlled by the brain [9] and is essential not only for the brain tofunction normally but also for the survival of animals. However, the phenomenon
itself is not well understood [10]. Currently, only electroencephalography (EEG),
an objective method, can detect sleep/wake states (wakefulness, non-rapideye-movement (NREM) sleep, and rapid-eye-movement (REM) sleep), whereas
behavioral methods, which record animals’ posture or movement, cannot [6].
c Springer Nature Switzerland AG 2019

I. V. Tetko et al. (Eds.): ICANN 2019, LNCS 11727, pp. 102–113, 2019.
https://doi.org/10.1007/978-3-030-30487-4_8

Sleep State Analysis Using Calcium Imaging Data by NMF

103

EEG recordings—signals from electrodes placed on the head—reﬂect extracellular electrical events across the cerebral cortex (the brain surface). Therefore,
the neurophysiological behavior of cortical neuronal populations are assumed to
diﬀer across sleep/wake states. Nevertheless, it is still unclear how cortical neural ensembles behave during sleep/wake states owing to the lack of measurement
and analysis technology. For example, EEG recordings reﬂect sleep/wake states;
however, individual neural activity cannot be observed using EEG. It is possible
to record neural activities at the cellular level with suﬃcient temporal resolution using electrophysiological techniques such as patch-clamp, intracellular, andextracellular unit recordings. However, it is diﬃcult to record a large number ofneurons or identify the recorded cell types using these techniques.
To overcome the technical limitations of electrophysiological methods, weuse calcium imaging techniques to observe the cortical neural activity during
sleep/wake states. Calcium imaging is a technique recently developed to record
neural activities at the cellular level. The ﬂuorescence intensity of calcium indicators such as GCaMP depends on the concentration of calcium ions. It reﬂects theneural activity because action potential generation (spike) increases the intracellular concentration of calcium. Another characteristic of GCaMP is that it
is genetically encodable and can be delivered into the target cells using virus
vectors and a Cre-LoxP system. Calcium imaging in the brain can be performedin vivo with two-photon laser scanning microscopy. Calcium imaging has the
disadvantage of a low time resolution (e.g., individual spikes cannot easily be
captured owing to the slow kinetics and low sampling rates.); however, it has
the considerable advantages of (1) high spatial resolution, (2) large recording
ﬁeld, and (3) ease of combination with genetic methods, where the activity ofvarious identiﬁed neurons can be obtained simultaneously.
In this study, we aim to understand sleep/wake-dependent neural ensembles in the cerebral cortex, which hopefully increase the understanding of sleep.
Calcium imaging data were acquired from identiﬁed excitatory and inhibitory
neurons in layer 2/3 of the primary motor cortex (M1)—a part of the cerebral
cortex—of a sleeping/waking mouse at 8 Hz. It is well known that individual
neurons in the M1 exhibit sleep/wake state-dependent activity patterns [8] andthe M1 contributes to memory consolidation of acquired motor skills during
sleep [5]. Their population behavior during sleep/wake states, however, is still
unexplored, despite its close involvement in learning and memory.
Two approaches can be considered to analyze neural ensembles; estimationof correlation and causality of neurons. We take the ﬁrst approach, estimating
the correlation of observed cortical neurons. To estimate causality, we need to
distinguish whether a presynaptic or postsynaptic cell is ﬁred from the data.
There are studies of inference on spikes [16] and neuronal connectivity [13] from
calcium imaging data. However, it is claimed in [13] that the calcium imaging
data sampled with a frequency lower than 30 Hz do not provide meaningful
results. Therefore, we can estimate only correlation from the calcium imaging
data sampled with low frequency.

104

M. Nagayama et al.

There are few reports on the statistical analysis of calcium imaging data inthe brain during sleep. One example of the sleep research using calcium imaging
is an analysis of neural activity in and near the lateral dorsal tegmental nucleus(LDT) of the brainstem that is a component of the REM-regulatory circuits [4].
Neuronal activity in the LDT for 20 s before and after the state transition was
analyzed by a principal component analysis (PCA). Then, k-means clustering
was performed for the 2D plot of the ﬁrst and second principal components.
A few LDT neurons were observed to be more active in the wake state thanthe REM state and vice versa, and all LDT neurons were less active in the
NREM state than in the REM or wake states. This research indicates thatneurons can be divided into a few groups and the behavior of each group diﬀers
across sleep/wake states. This study focuses on the state transition and extracted
the principal components of the ﬂuorescence intensity during state transition.
However, neuronal groups throughout sleep/wake states should be analyzed forour aim.
To better understand sleep/wake-dependent neural ensembles in the cerebral
cortex, we assume that groups of cortical neurons become active simultaneously
and the activated groups change over time. We used non-negative matrix factorization (NMF [11]) to estimate those groups and how they are activated intime. NMF was not performed to extract the ﬂuorescence intensity like Cox [4]did, but for the neuronal groups. We used NMF because neural activity can be
expressed by the sum of individual neuronal activity; further, ﬂuorescence intensity data are always positive values. We found that there are certain groups ofneurons that behave diﬀerently between sleep and wake states.
We proposed a protocol to analyze the calcium imaging data with a low
sampling rate. We showed that the well-known NMF algorithm could eﬀectively
analyze such data by comparing its performance with those of PCA and independent component analysis (ICA).

2

Methodology

We acquire calcium imaging data of neurons in layer 2/3 from the primary motor
cortex (M1) of a sleeping/waking mouse. The imaging lasted for 15 min at 8 Hz.
A total of 154 neurons were observed. The ﬂuorescence signal in each frame foreach neuron was extracted by evaluating the mean intensity of the pixels withineach region of interest (ROI) after subtracting the background signal. EEG and
EMG signals were also recorded during the experiments, and NREM, REM, andwake states were scored.
We made an assumption regarding the generative model of the data, which isthat there are K groups of observed neurons that are activated simultaneously,
and the ﬂuorescence intensity of observed neurons is generated by the combined
activities of those groups. Groups can overlap and neurons can belong to multiple
groups. We assume that the behavior of the groups vary over time and theactivities of those groups diﬀer between NREM, REM, and wake states. The
schematic is shown in Fig. 1.

Sleep State Analysis Using Calcium Imaging Data by NMF

Fig. 1. Observed neurons are assumed toform K-grouped and the recorded ﬂuorescence intensity is generated by the combined activities of those groups. In thisﬁgure, K = 3 groups are assumed andthe activated neurons are colored.

105

Fig. 2. Mathematical model designed
based on the assumption Fig. 1. The ﬂuorescence intensity of the observed neurons at time t, x(t), is represented as aweighted sum of K groups dk . The ﬁgureshows a schematic of the model whenK = 3.

Based on the assumption above, we designed a mathematical model for the
calcium imaging data. Let R+ be a set of non-negative real numbers. Let x(t) ∈RN+ , (t = 1, 2, · · · , T ) be the ﬂuorescence intensity of the observed neurons attime t and dk ∈ RN
+ , (k = 1, 2, · · · , K) be K groups of neuronal activities where
N is the number of observed neurons and T is the observation time. Fluorescenceintensity data are all positive values; therefore, x(t) is a positive vector. Then,x(t) can be modeled by the weighted sum of the groups {dk ; k = 1, · · · , K} as
follows:K
ck (t)dk + η(t),(1)
x(t) =
k=1

where ck (t) is a positive coeﬃcient of dk and η(t) ∈ RN is a noise vector at timet. The schematic of the model is shown in Fig. 2. We note that the groups dk areconstant over time and time dependency of the observed data x(t) is represented
by the coeﬃcient ck (t) and noise η(t). The noise of calcium imaging is mostly
photon shot noise, which obeys a Poisson distribution. For high photon counts,
shot noise can be approximated by a Gaussian distribution [15]; therefore, η(t)can be represented as(2)
η(t) ∼ N (0, σ 2 ),where σ 2 is the variance.
Equation (1) can be written in a matrix form as follows:
X = DC + H,

(3)

×T
×K
where X ∈ RN
, D ∈ RN
, C ∈ RK×T
, and H ∈ RN ×T . The columns of X+++
are the observation x(t), columns of matrix D are the groups dk , (k, t)-element

106

M. Nagayama et al.

of matrix C are the coeﬃcients ck (t), and the columns of H are noise vectors
η(t). The correspondence between Eqs. (1) and (3) is shown in Fig. 3.

Fig. 3. Matrix form of Eq. (1).

We estimate the matrices D and C from the data matrix X by using
NMF [11], which is an unsupervised learning method similar to PCA and ICA,
as one of the approaches in the linear generalized component analysis [3]. NMF
decomposes a non-negative matrix X into a product of non-negative matrices
D and C. By decomposing the data matrix X by NMF, frequent signals will
be obtained as the columns of D. The matrix D is called a dictionary in theliterature of signal processing and machine learning. The columns of the dictionary D are interpreted as basis vectors or also called atoms. The columns of thematrix C are coeﬃcients of atoms and we call C a coeﬃcient matrix henceforth.
Because of the non-negativity of the coeﬃcient matrix C, X is represented by anon-negative weighted sum of atoms of D. The restriction of non-negativity also
induces sparsity of C. The number of groups K ∈ N should be given a priori fordecomposition.
From Eqs. (1) to (3), Xij is i.i.d. Gaussian random variable with mean [DC]ij
and variance σ 2 :(4)p(Xij ; [DC]ij ) = N ([DC]ij , σ 2 ).
The dictionary D and coeﬃcient matrix C is estimated by the following maximum log-likelihood problem under the non-negative constraints on D and C:log p(Xij |[DC]ij ).
(5)maximizeD≥0,C≥0i,j

By calculating Eq. (5), we can see that Eq. (5) is equivalent to minimizing the
Frobenius norm of the diﬀerence between X and DC under the non-negative
constraints on D and C:
minimize ||X − DC||F ,D≥0,C≥0where the Frobenius norm of matrix A ∈ Rm×n is deﬁned as follows:


n
m 
||A||F = 
a2ij , aij is the (i, j)-element of A.
i=1 j=1(6)(7)

Sleep State Analysis Using Calcium Imaging Data by NMF

3

107

Results

We used two-photon calcium imaging data of a transgenic mouse during
sleep and wakefulness. To distinguish between excitatory (glutamatergic) andinhibitory (GABAergic) neurons, Vgat-tdTomato mice were generated. Red ﬂuorescence protein tdTomato was expressed particularly in inhibitory neurons
in layer 2/3 of the M1 of Vgat-tdTomato mice, which were histologically conﬁrmed. A genetically encoded green ﬂuorescence calcium indicator GCaMP6s
was used to monitor neural activity. GCaMP6s was delivered to M1 neurons
using adeno-associated virus vectors under the control of a neuron-speciﬁchuman synapsin 1 promoter. Thus, we deﬁned GCaMP6s-positive/tdTomatonegative and GCaMP6s-positive/tdTomato-positive neurons as excitatory andinhibitory neurons, respectively. Fluorescence imaging was performed with acustom-designed upright two-photon microscope (based on Axio Examiner
Z1/LSM780) and a trackball-treadmill system [10]. The mice were acclimated tosleeping on the trackball-treadmill prior to imaging experiments. GCaMP6s andtdTomato were two-photon excited at 910 and 1040 nm, respectively, using a tunable Ti:Sa laser (Maitai DeepSee, Spectra-Physics). Fluorescence of GCaMP6s
and tdTomato was detected with a non-descanned GaAsP detector in the rangeof 500–550 nm and >555 nm, respectively (BiG, Zeiss). A 1040-nm excitation
was used only to identify inhibitory neurons prior to time-lapse calcium imaging.
Two-photon time-lapse images were intermittently acquired six times between
12:00 and 17:00 (the rest time for mice) at 8 frames/s with 128 × 256 pixels of 16bit depth. Each imaging lasted for 15 min and we refer to each 15 min data as adataset. EEG and EMG signals were also recorded during imaging experiments.
EEG signals were ampliﬁed 40,000x and ﬁltered with a pass-band of 0.5–500 Hz,
and EMG signals were ampliﬁed 4000x and ﬁltered with a pass-band of 1.5–1000 Hz using an analog ampliﬁer (MEG-5200, NIHON KOHDEN). EEG/EMG
signals were digitized at 2000 Hz using a 16-bit analog-to-digital converter (Digidata 1440A, Molecular Devices), and they were acquired with Clampex 10.3 software (molecular devices). Scan timing signals for two-photon imaging were also
digitized with the same systems to temporally match the EEG/EMG data to thetwo-photon images. EEG and EMG signals were downsampled at 250 Hz, divided
into 4 s epochs, and analyzed for scoring NREM/REM/wake states using a homemade Matlab program. The ﬂuorescence signal in each frame for each neuron was
extracted by evaluating the mean intensity across the pixels within each region ofinterest (ROI) after subtracting the background signal using Fiji/ImageJ. ROIs
were set such that the neuropil signals were not included.
We calculated ΔF/F from each dataset and the data matrix X was createdby concatenating six datasets. Then, the data matrix X was normalized such thatthe sum of each column was one. By minimizing Frobenius norm cost function,
the dictionary D and coeﬃcient matrix C were obtained. We tried 29 diﬀerent
values for the number of atoms K, increasing 10 to 150 in increments of 5, because
the appropriate value for K is unknown. Larger K can reconstruct X with small
error; however, the atoms would have small information. We want atoms that

108

M. Nagayama et al.

have rich information; thus, we assumed a smaller K than the number of neurons.
We also used PCA and ICA to validate the performance of NMF.
We used random forest [1] to solve the state classiﬁcation problem to compare
three methods and diﬀerent K. The input of random forest are columns of thecoeﬃcient matrix C for NMF. For PCA and ICA, matrices equivalent to C wereused for the input of random forest. High accuracy in random forest means thatthe extracted features of the corresponding method can distinguish states. We
are comparing three methods and K by the accuracy of random forest. Threestates, NREM, REM, and wake, were classiﬁed using 3-fold cross-validation.
Train and test data were prepared after dividing the input data into 15 s blocks.
The results are shown in Fig. 4. The accuracy of each validation is shown bymarkers, and the average of the three validations is shown by lines.

Fig. 4. Accuracies of state classiﬁcation by the random forest with three methods anddiﬀerent numbers of atoms. Features used in the random forest are columns of thecoeﬃcient matrix C and matrices correspond to that. Three states, NREM, REM, andwake, were classiﬁed and the classiﬁcation accuracy is calculated using 3-fold crossvalidation. The accuracy of each validation is shown with markers and the average ofthree validations is shown with lines.

Sleep State Analysis Using Calcium Imaging Data by NMF

109

The result shows that NMF exhibits better performance than PCA and ICA.
This suggests that NMF is more suitable for the generative process of the datawith its non-negativity. The result also shows that K = 65 had the highest
accuracy for NMF, suggesting that this is the best K that extracts atoms that
diﬀer across sleep/wake states.
We then analyzed the atoms of K = 65. To evaluate the atoms used diﬀerently between sleep and wake states, we used the importance of random forest
and Jensen–Shannon divergence [12] (JSD). We used JSD because analyzing theimportance of random forest might not be suﬃcient. The k-th row of C, ck (t),is a time series of the coeﬃcient of the atom dk . When values of ck (t) diﬀer
between states, the corresponding atom dk is considered to be used diﬀerently
between states. To measure how diﬀerently an atom was used in each state, we
calculated the JSD of the coeﬃcient of each atom between NREM and wakestates. REM state was not used because it rarely appears in our dataset.
The importance from random forest and the JSD of each atom are presented
in Figs. 5 and 6 respectively. Top 12 atoms that had high importance and the
JSD is shown in Table 1. The 58th atom is the highest and the 18th atom is thesecond highest in both importance and the JSD. However, the order of atoms
after third place diﬀered in importance and the JSD. This diﬀerence should be
carefully considered; however, we could not determine the same, and thus, wewould leave it for our future work.

Fig. 5. The importance of K = 65 atoms
from random forest.

Fig. 6. JSD calculated from coeﬃcients of K = 65 atoms.

Our aim is to analyze sleep/wake-dependent neural ensembles. Therefore, weused Otsu’s method to obtain the number of neurons responsible in each atom.
The result is shown in Fig. 7. Top 12 atoms from Table 1 is colored with light
blue. 58th and 23rd atoms had more than 60 responsible neurons. Both atoms
had high importance and the JSD.

110

M. Nagayama et al.
Table 1. Top 12 atoms that had high importance and the JSD.
MethodAtom

The importance 58 18 17 36 43 53 23 10 48 4 33
JSD

58 18 11 48 41 40

1

1 23 17 8 53 36

Fig. 7. Number of neurons responsible in each atom using Otsu’s method. Top 12
atoms from Table 1 is colored with light blue. (Color ﬁgure online)

We picked up 6 atoms from Table 1 that had more than one responsible
neuron. The spatial plots of those atoms are shown in Fig. 8. Red markers represent excitatory neurons and blue markers represent inhibitory neurons. The
size of a marker is the intensity of the corresponding neuron. The boxplot of thecoeﬃcients of those atoms in NREM and wake is shown in Fig. 9.
From Fig. 8, 58th , 23rd , 17th , and 4th atoms do not have evident spatial
clusters; however, 48th and 11th atoms do. From Fig. 9, the former four atoms
contribute to NREM, whereas the latter two atoms somehow contribute to wake.
Co-activation of nearby neurons emerges in the M1 for stereotyped simple behaviors [7] but not for learned voluntary movements [14]. In our analysis, sleeppreferred neuronal ensembles exhibited no evident spatial clustering, suggestingthat the M1 is devoted to higher-order information processing for acquiring complex behaviors during sleep.

Sleep State Analysis Using Calcium Imaging Data by NMF

111

Fig. 8. The spatial plot of atoms. Each marker represents a neuron. The size of amarker is the intensity of the corresponding neuron.

Fig. 9. The boxplot of the coeﬃcients of the atoms in NREM and wake.

112

4

M. Nagayama et al.

Conclusion

In this paper, we proposed a protocol to analyze the calcium imaging data witha low sampling rate. We decomposed the data matrix X with NMF and validated its performance by comparing the results obtained with those of PCA andICA. The results showed that NMF performed better than PCA and ICA by
extracting atoms whose activities diﬀer between states. This suggests that NMF
is more suitable for the generative process of the data with its non-negativity.
Results also showed that there are groups of neuronal activities that become
active simultaneously and some groups are related to sleep/wake states. Thisis a notable result because NMF does not use state information to decompose
the data matrix X. Furthermore, atoms contributing to NREM had no evident
spatial clusters, which suggests that the M1 is devoted to high-order information
processing for acquiring complex behavior during sleep.
For further investigation, two things should be considered. First, a method
to determine the number of atoms K should be developed. NMF requires thatthe number of atoms K be determined beforehand. In this paper, we used arandom forest to determine the number of atoms K; however, further consideration should be made from a biological perspective. Second, the characteristicof calcium imaging should be considered. Calcium imaging has a slow responseto neuronal activities. The decay halftime of the calcium indicator GCaMP6s isapproximately 0.5s [2], which is much higher than the sampling rate of the data
we use. Therefore, we can expect better results if this property is incorporated
into the decomposition model of NMF. For example, adding a constraint thatmakes two adjacent coeﬃcients ck (t) and ck (t + 1) similar would help eﬀectively
obtain a better dictionary D and its coeﬃcient C.
Acknowledgements. This work was supported by Grants-in-Aid for ScientiﬁcResearch (KAKENHI), Japan Society for the Promotion of Science (JSPS) (Grant
Number 16K18358 to T.K.; 26220207 to T.K. and M.Y.; 19K12111 to H.H; 17H06095
to M.Y.); World Premier International Research Center Initiative (WPI), the Ministry of Education, Culture, Sports, Science and Technology (MEXT) (to M.Y.); Core
Research for Evolutional Science and Technology (CREST), Japan Science and Technology Agency (JST) (Grant Number JPMJCR1761 to H.H.; JPMJCR1655 to M.Y.);
Yamada Research Grant (to T.K.), Takeda Science Foundation (to M.Y.), and Uehara
Memorial Foundation (to M.Y.).

References
1. Breiman, L.: Random forests. Mach. Learn. 45(1), 5–32 (2001). https://doi.org/
10.1023/A:1010933404324
2. Chen, T.W., et al.: Ultrasensitive ﬂuorescent proteins for imaging neuronal activity.
Nature 499(7458), 295–300 (2013). https://doi.org/10.1038/nature12354
3. Cichocki, A., Zdunek, R., Phan, A.H., Amari, S.: Nonnegative Matrix and Tensor
Factorizations. Wiley, Chichester (2009)
4. Cox, J., Pinto, L., Dan, Y.: Calcium imaging of sleep-wake related neuronal activity
in the dorsal pons. Nat. Commun. 7(1), 10763 (2016). https://doi.org/10.1038/
ncomms10763

Sleep State Analysis Using Calcium Imaging Data by NMF

113

5. Dayan, E., Cohen, L.G.: Neuroplasticity subserving motor skill learning. Neuron
72(3), 443–54 (2011). https://doi.org/10.1016/j.neuron.2011.10.008
6. Deboer, T.: Technologies of sleep research. Cell. Mol. Life Sci. 64(10), 1227 (2007).
https://doi.org/10.1007/s00018-007-6533-0
7. Dombeck, D.A., Graziano, M.S., Tank, D.W.: Functional clustering of neurons inmotor cortex determined by cellular resolution imaging in awake behaving mice.
J. Neurosci. 29(44), 13751–13760 (2009). https://doi.org/10.1523/JNEUROSCI.
2985-09.2009
8. Evarts, E.V.: Temporal patterns of discharge of pyramidal tract neurons during
sleep and waking in the monkey. J. Neurophys. 27, 152–71 (1964). https://doi.org/
10.1152/jn.1964.27.2.152
9. Hobson, J.A.: Sleep is of the brain, by the brain and for the brain. Nature
437(7063), 1254–1256 (2005). https://doi.org/10.1038/nature04283
10. Kanda, T., et al.: Sleep as a biological problem: an overview of frontiers in sleep
research. J. Physiol. Sci. 66(1), 1–13 (2016). https://doi.org/10.1007/s12576-0150414-3
11. Lee, D.D., Seung, H.S.: Learning the parts of objects by non-negative matrix factorization. Nature 401(6755), 788–791 (1999). https://doi.org/10.1038/44565
12. Lin, J.: Divergence measures based on the Shannon entropy. IEEE Trans. Inform.
Theory 37(1), 145–151 (1991). https://doi.org/10.1109/18.61115
13. Mishchencko, Y., Vogelstein, J.T., Paninski, L.: A Bayesian approach for inferring
neuronal connectivity from calcium ﬂuorescent imaging data. Ann. Appl. Stat.
5(2B), 1229–1261 (2011). https://doi.org/10.1214/09-AOAS303
14. Peters, A.J., Chen, S.X., Komiyama, T.: Emergence of reproducible spatiotemporal
activity during motor learning. Nature 510(7504), 263–267 (2014). https://doi.
org/10.1038/nature13235
15. Sjulson, L., Miesenböck, G.: Optical recording of action potentials and other discrete physiological events: a perspective from signal detection theory. Physiology
22(1), 47–55 (2007). https://doi.org/10.1152/physiol.00036.2006
16. Vogelstein, J.T., Watson, B.O., Packer, A.M., Yuste, R., Jedynak, B., Paninski,
L.: Spike inference from calcium imaging using sequential Monte Carlo methods.
Biophys. J. 97(2), 636–655 (2009). https://doi.org/10.1016/J.BPJ.2008.08.005