Spike to calcium models can be used to test our understanding of the process bridging spikes and calcium imaging.

In the context of simultaneously recorded spiking and imaging the synthetic calcium signal generated based on the recorded series of spikes can be compared directly to the recorded calcium signals.

In the context of non-simultaneously recorded data, spike to calcium models can test how much of the difference between the same analysis performed on a spiking data and on an imaging data can be explained by the indirect report of activity. A synthetic imaging dataset is generated from the recorded spike trains. Then the analysis of interest is applied on the synthetic data and how much of the difference between spike and imaging is quantified.

All of our spike to calcium models consists of two steps. First, spikes at times $\{t_k\}$ are converted to a latent variable, $c(t)$, by convolution with a double-exponential kernel: $$c(t)=\sum_{t>t_k}\exp(-\frac{t-t_k}{\tau_d})[1-\exp(-\frac{t-t_k}{\tau_r})]+n_i(t)$$ where $\tau_r$and $\tau_d$ are the rise and decay times, respectively. $n_i(t)\sim N(0, \sigma_i^2)$ is Gaussian distributed "internal" noise with negative values were set to zero. Second, c(t) was converted to a synthetic fluorescence signal through a general function $f(\cdot)$: $$F_{synth}(t) = f(c(t)) + n_e(t)$$ where $n_e(t)\sim N(0, \sigma_e^2)$ is Gaussian distributed "external" noise.
Note that since synthetic data contains no baseline activity, the synthetic imaging data itself $F_{synth}(t)$ is the signal that is relevant to compare to $\Delta F/F$ in data

## S2C Linear

$$f(c(t)) = F_{\max}\cdot c(t) + F_0$$ where $F_{\max}$ is the maximum possible fluorescence change, $F_{0}$ is the mimimun possible fluorescence change.

## S2C Hill

$$f(c(t)) = F_{\max}\frac{c(t)^n}{c(t)^n+K_d}$$ where $n$ is a non-linearity, $K_d$ is a half-activation parameter, $F_{\max}$ is the maximum possible fluorescence change.

## S2C Sigmoid

$$f(c(t)) = F_{\max}\frac{1}{\exp[-k(c(t)-c_{1/2})]+1}$$ where $k$ is a non-linearity, $c_{1/2}$ is a half-activation parameter, $F_{\max}$ is the maximum possible fluorescence change.

The inverse problem, inferring the times at which a neuron fired a spike given recorded calcium imaging data, is far more difficult.

We tested tens of calcium to spike model. In general, calcium to spike models posit a generative model for the imaging data given a spike time series, and then infer the spike time series by maximizing the likelihood of the data. The generative model typically follows a similar mathematical setup to that of the calcium to spike models:
First, the (unknown) spikes at time $t$ (usually in a small time bin, $dt$) are converted to a latent variable, c(t), using an n-step autoregressive, process, AR(n): $$c(t) = \sum_{i=1}^n a_i \cdot c(t-i) + spk(t)$$ where $n = 1$, $a_1 = 1-\frac{dt}{\tau_d}$, $spk(t)$ is the number of spikes in time bin $dt$ (usually zero or one). Second, c(t) is converted to a synthetic fluorescence signal through a general function $f(\cdot)$: $$F_{synth}(t) = f(c(t)) + n_e(t)$$. As before since synthetic data contains no baseline activity, the synthetic imaging data itself $F_{synth}(t)$ is the signal that is relevant to compare to $\Delta F/F$ in data

## C2S Nonnegative Wienner Filter

AR(1) process: $n = 1$.
Non-negative spike counts: $spk$ is non-negative (not strictly to be integer).
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: linear.
Reference Vogelstein et al. (2010) Fast nonnegative deconvolution for spike train inference from population calcium imaging. J. Neurophysiol.

## C2S FOOPSI

AR(1) process: $n = 1$.
Poisson spike counts: $spk$ follows "Poission" distribution, which is approximated with an Exponential (not strictly to be integer).
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: linear.
Reference Vogelstein et al. (2010) Fast nonnegative deconvolution for spike train inference from population calcium imaging. J. Neurophysiol.

## C2S Finite Rate of Innovation

AR(1) process: $n = 1$.
Spike counts: strictly to be integer and the number of spike is decided using theory of finite-rate innovation signals.
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: linear.
Reference Oñativia et al. (2013) A Finite Rate of Innovation algorithm for fast and accurate spike detection from two-photon calcium imaging. J. Neural Eng.

## C2S Constrained OOPSI AR(n)

AR(n) process: $n$ is taken from 1, 2, or 3.
Poisson spike counts: $spk$ follows "Poission" distribution.
Noise structures: $n_e$ follows Gaussian distribution; empirical estimation of noise prior.
$f(\cdot)$ function: linear.
Reference Pnevmatikakis et al. (2016) Simultaneous Denoising, Deconvolution, and Demixing of Calcium Imaging Data. Neuron.

## C2S Constrained OOPSI MCMC

AR(1) process: $n = 1$.
Poisson spike counts: $spk$ follows "Poission" distribution, which is sampled using discrete binary sampler using MCMC.
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: linear.
Reference Pnevmatikakis et al. (2014). Bayesian spike inference from calcium imaging data. Asilomar Conf. on Signals, Systems, and Computers.

## C2S Peel Linear

AR process: the same as spike-to-calcium model with $\tau_r$ and $\tau_d$
Spike counts: strictly to be integer and using template-matching to search the optimal number.
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: linear.
Reference Grewe et al. (2010) High-speed in vivo calcium imaging reveals neuronal network activity with near-millisecond precision. Nat. Methods

## C2S Peel Nonlinear

AR process: the same as spike-to-calcium model with $\tau_r$ and $\tau_d$
Spike counts: strictly to be integer and using template-matching to search the optimal number.
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: hill function with n=1.
Reference Grewe et al. (2010) High-speed in vivo calcium imaging reveals neuronal network activity with near-millisecond precision. Nat. Methods

## C2S MLSpike

AR process: the same as spike-to-calcium model with $\tau_d$
Spike counts: strictly to be integer and using template-matching to search the optimal number.
Noise structures: $n_e$ follows Gaussian distribution.
$f(\cdot)$ function: hill function with n=1.
Reference Deneux et al. (2016) Accurate spike estimation from noisy calcium signals for ultrafast three-dimensional imaging of large neuronal populations in vivo. Nat. Commun.

## Hidden dynamics in calcium imaging

Individual neurons often exhibit diverse temporal dynamics and even change selectivity. We classified neurons into three broad categories: ‘monophasic’ neurons showed consistent selectivity across the trial; ‘multiphasic’ neurons changed selectivity over time; ‘non-selective’ neurons responded similarly across trial types but in a task-modulated manner.

For two-alternative force choice task, neural selectivity for type A or type B trials was using two-sample t-tests, with neural activity binned over 67 ms, which corresponds to one imaging frame (~15Hz). A neuron was selective if it showed selectivity (p < .05) for >335 ms (5 continuous frames). A selective neuron was multiphasic if the polarity of selectivity switched, with periods of selectivity at least 335 ms long. Other selective neurons were monophasic. The fraction of each cell type is then computed for each recording and model condition.

## Peakiness of neural dynamics

Slow and nonlinear dynamics of imaging can lead to a distorted view of neuronal dynamics, where 'peaky' neural dynamics at transitions between behavioral epochs in ephys seemed to be delayed and jittered, causing a sequence-like appearance in imaging.

The ‘peakiness’ of the distribution of neuronal activity (‘s’) across recording modalities were measured as the difference between observed neural activity and temporally uniformly distributed neural activity ($P=\frac{1}{2T}$; T is the length of the recording/task): $$s=\frac{1}{P}\sqrt{P\sum_{i=A,B}\int_0^T dt(P_i(t)-P)^2}$$

## Decodability

Linear discriminant analysis (LDA) is applied on neural dynamics grouped into bins corresponding to single imaging frames (67 ms) to compute the instantaneous decodability of trial type. The optimal LDA decoder was computed separately for each time bin, using correct trials. We estimated performance for the instantaneous LDA decoder by sampling a subset of units (10-200), and averaging 100 samplings. We separated the trials of each neuron into non-overlapping training (70%) and testing (30%) sets. We then shuffled trial identity of each neuron to build up random collections of population activity at each time point for correct trial A (50% in training and testing) and trial B using all recorded neurons. The instantaneous decoder of trial type was computed from training set and its performance was evaluated on the testing set.

To achieve robust estimation given the large number of neurons, we applied a sparse version of LDA.

## Variance content

In this study, PCA was performed on the activity of neurons averaged across two alternative trial types ($s\in\{A, B\}$): $$\mathbf{r}(s,t)=\mathbf{Cx}(s,t)+\langle\mathbf{r}(s,t)\rangle_{s,t}$$ $\mathbf{r}$ is a $n\times 2T$ matrix, where n is the number of recorded units in each dataset and T is the number of time points for each trial type. $<\mathbf{r}(s,t)>_{s,t}$ a vector of the mean activity of each neuron across time and trial types. $\mathbf{x}(s,t)$ is an $n\times 2T$ PC score matrix, where the $i$th row corresponds to the $i$th PC score. We estimated the relative contribution to each PC of the different forms of variance: temporal dynamics, trial-type selectivity and other. Explained variance (EV) of temporal dynamics $EV_i(t)$ and trial-type selectivity $EV_i(s)$ for the th principal component (PC) were computed as: $$EV_i(t)=\langle\langle x_i(s,t)\rangle_s^2\rangle_t/\langle x_i(s,t)^2\rangle_{s,t}$$ $$EV_i(t)=\langle\langle x_i(s,t)\rangle_t^2\rangle_s/\langle x_i(s,t)^2\rangle_{s,t}$$ respectively.