Low dimensional spontaneous activity proceeds throughout stimulus presentation

We first considered two-photon calcium imaging data from the optic tectum of the developing larval zebrafish (Fig 1A). Fish expressing the genetically encoded calcium indicator GCaMP6s were embedded in agarose while small dark spots were presented at systematically varying angles across the visual field [17]. Onset of the visual stimulus evoked calcium transients in the tectum (Fig 1B) consistent with the topography of the retinotectal map. The presentation of a spot was followed by an interval of 19s without any stimulation, enough time for calcium levels to return to baseline before the next stimulus. The optic tectum was often highly active during these inter-stimulus intervals (Fig 1C) and calcium transients sometimes occurred spontaneously just before stimulus onset, elevating the recorded fluorescence levels associated with that stimulus.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Spontaneous activity in calcium imaging data.

(A) Two-photon calcium imaging of the larval zebrafish optic tectum. NP, neuropil; PVL, periventricular layer. (B) Fluorescence traces from 10 example neurons. Dashed vertical lines indicate stimulus onset; colour represents azimuth angle of presented stimulus. (C) Example fluorescence trace segment illustrating that spontaneous calcium transients can occur just before stimulus onset, inflating stimulus-response estimates.

https://doi.org/10.1371/journal.pcbi.1008330.g001

One approach to separating SA from EA would be to repeatedly present the same set of visual stimuli over many trials, compute the peri-stimulus time histogram (PSTH) across trials, and assign all calcium transients that deviate from the PSTH as “spontaneous”. However, this approach cannot handle stimulus sequences that occur within single trials or in a randomised order unless the PSTH is calculated over a short window surrounding the time of stimulus onset, in which case the estimated SA is biased by edge effects. To the authors’ knowledge there are no standard existing methods capable of separating SA from EA on a single-trial basis.

As our first attempt we therefore considered using a novel combination of existing tools. We characterised the stimulus-driven component of the population activity by simply expressing each fluorescence trace as a linear combination of regressors that specify the basic shape and timescale of calcium activity. We defined a basis of stimulus regressors [18, 19] by convolving a calcium impulse response kernel with the presentation times of each stimulus and performed a multivariate linear regression of the population data onto this basis set using non-negative least squares (S1A–S1C Fig). Any highly structured activity in the residual data would likely be driven primarily by latent sources of SA not directly related to the onset of the presented stimuli. We then applied non-negative matrix factorisation [20] (NMF) to search for low dimensional structure in the residuals. NMF attempted to reconstruct the residual population data as the product of two matrices with non-negative entries: a matrix whose rows are timeseries that capture patterns of SA shared across groups of neurons, and a matrix whose columns describe how neurons are coupled to such timeseries (S1D Fig). The NMF description of SA identified low dimensional structure that proceeded throughout the recording, largely independent of the stimulus (S1E Fig).

While this residual NMF approach is efficient due to highly optimised computational routines, it is limited by two characteristics. First, because the models of stimulus processing and shared SA are not inferred jointly, receptive field estimates are biased towards higher values by spontaneous calcium transients that coincide with stimulus presentations. This in turn introduces bias when applying NMF to the residual data, since some of the contribution of SA was already subtracted out at the receptive field estimation stage. Second, NMF has no model for the highly stereotyped structure of calcium transients, and therefore does not respect this structure in the components that it finds (S1F Fig).

CILVA simultaneously captures evoked responses and shared spontaneous activity

To overcome these difficulties, we instead developed a generative statistical model that describes evoked and spontaneous activity simultaneously rather than sequentially (Fig 2A and 2B). Our method works directly with filtered fluorescence traces of individual neurons rather than raw calcium imaging videos, and therefore can be used after segmenting cells with popular calcium imaging preprocessing packages such as CaImAn [21] and Suite2p [22]. While these packages provide spike deconvolution modules, we opted to work with fluorescence traces as many calcium imaging datasets (including those analysed in this paper) have low temporal resolution, rendering precise estimation of spike times difficult. With this in mind, the model specifies the observed fluorescence level f_n(t) for neuron n at time t in terms of the underlying calcium concentration c_n(t), where the scalars α_n and β_n determine the scale and baseline of the fluorescence signal respectively, and ϵ_n(t) represents Gaussian noise. Consistent with experimental data [6] and previous models for calcium imaging [7, 23], the calcium dynamics are assumed to be highly stereotyped and are defined by the convolution of a GCaMP impulse response kernel k with a vector of calcium influxes λ_n (analogous to an activity intensity function). The kernel k is a difference-of-exponentials function (see Methods), which includes both rise and decay time constants. We found that including an explicit rise time was essential as GCaMP6s activity is slow relative to the sampling rate of many optical imaging systems.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Overview of the CILVA approach for decoupling stimulus-evoked responses and latent sources of SA.

(A) Proposed generative architecture underlying multivariate calcium imaging data. Neurons are driven by sensory stimuli (red) and latent sources of SA (blue). These two sources are combined additively to define the underlying rate of calcium influx (λ_n), before being convolved with a GCaMP kernel. Calcium levels are subsequently reported through noisy fluorescence intensities. (B) The intensity of calcium influx λ_n encoding stimuli and shared SA is convolved with a GCaMP kernel k to generate observed calcium levels. (C) The learned encoding model provides a method for decoupling evoked responses from common patterns of SA.

https://doi.org/10.1371/journal.pcbi.1008330.g002

The key ideas of the model are that (i) evoked responses will tend to be locked to the onset of the stimulus, (ii) evoked responses typically have a simple impulse response structure in calcium imaging data, and (iii) neural activity not attributable to evoked activity should be explained as far as possible by SA with a specific structure. The rate of calcium influx λ_n(t) in each imaging frame t was thus assumed to be driven by the addition of two underlying non-negative sources: processing of the stimulus s(t) through a linear receptive field w_n, and a small number of unobserved or “latent” sources of SA x(t), Here x(t) is the low dimensional latent state at time t, b_n is a vector describing how neuron n is affected by these factors, and ⋅^⊤ denotes the transpose operation. In the event that neurons exhibit prolonged neural responses, the stimulus design matrix s can be straightforwardly modified by including copies of each stimulus shifted in time [24]. This SA model is inspired by the application of factor analysis methods to neural population data [9, 25], which posit that low dimensional structure arises from latent computations or brain states that concurrently affect subsets of neurons. However, the latent factors underlying spontaneous calcium transients in our model differ mathematically from classical factor analysis in that, due to non-negativity of the calcium levels, factor activity states and coupling between latent factors and neurons must be non-negative [26].

We fit the model by computing the maximum a posteriori estimate of the latent factor activity states. Because these activities were less constrained by the model compared to the time-locked evoked responses (and therefore likely to be more complex) we encouraged sparsity by placing a non-negative prior on the latent factors with high density near zero, and used a simple model selection procedure to estimate the sparsity penalty (see Methods).

Fluorescence signals can be decomposed into their evoked and spontaneous components

Our model can be used to analyse the separate contributions of evoked and spontaneous activity to the observed fluorescence levels (Fig 2C). The fitted model can be succinctly summarised by the equation where denotes an estimated variable, * denotes linear convolution, and 1_T is a vector of ones with length equal to the number of imaging frames T. Here we have dropped explicit dependence on the calcium levels c_n(t), which are deterministic given the other model parameters. The components of the signal driven purely by evoked or spontaneous activity can then be extracted from the convolution to give and We first verified the model on simulated data with known ground truth, modelling the properties of the zebrafish and mouse data that we subsequently consider (S1 and S2 Tables, S2 and S3 Figs). We then applied the model to our calcium imaging data from the zebrafish optic tectum to decouple the evoked and spontaneous calcium transients (Fig 3). These data were segmented using custom software, but we also verified that our results do not depend on the preprocessing package for source extraction by performing the same analysis on data processed with CaImAn (S4 Fig). Overlaying these decoupled calcium traces onto the experimental data, we found that they provided realistic descriptions of calcium activity (Fig 3A) and a close fit between the raw fluorescence trace and the model reconstruction (Fig 3B, S5 Fig). The time-locked responses to stimuli were well modelled by , while low dimensional SA was identified by the projected latent factor activity (Fig 3A, S6 Fig, S1 and S2 Videos). Neural activity in the residual data (i.e., after subtracting the model reconstruction from the raw fluorescence traces) arose primarily from spontaneous calcium transients that were independent of the latent sources of shared SA (and were therefore attributed to private, as opposed to shared, variability [27], S7 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Fitted model components for the zebrafish shown in Fig 1.

(A) Results of fitting CILVA and decoupling EA (red) and shared SA (blue) in an experimental recording. Inset numbers denote the Pearson correlation coefficient between raw fluorescence trace and model fit. The 10 neurons with the highest correlations between data and model fit are shown. (B) Distribution of correlation coefficients between data and model fits. Shuffled data obtained by cyclically permuting each trace by a random offset while preserving its temporal structure. (C) Inferred latent factor timeseries. Inset numbers denote the factor contribution indices, defined as the mean reduction in correlation coefficient across the population following deletion of the corresponding factor. (D) Cross-validated distributions of correlation coefficients for 1 to 15 latent factors. Shaded error bars indicate one standard deviation. For this fish adding additional latent sources of SA beyond three factors provides little improvement in the average correlation for both training and held-out test data. (E) Estimated factor coupling matrix shows that latent factors target distinct, non-overlapping sets of neurons. (F) Cumulative factor contribution indices for 0 to 3 latent factors. (G) Cross-correlograms show little interaction between latent factors and no long term structure in individual factor activity. While zebrafish can occasionally exhibit factors with secondary peaks in their autocorrelation plots (see e.g. S9H Fig), the location of these peaks varies between factors and fish, and do not align with the interstimulus interval. (H) Example estimated stimulus tuning curves (red). Tuning curves obtained by trial-averaging fluorescence levels over a window following stimulus onset (gray) provided for comparison (Methods). Shaded error bars denote one standard deviation. Full temporal traces for these neurons are given in S8 Fig. (I) Retinotopic maps obtained by trial-averaging fluorescence levels over a window following stimulus onset (left) and by CILVA-estimated tuning curves (right). The CILVA map is more refined since the estimated stimulus filters already account for ongoing SA.

https://doi.org/10.1371/journal.pcbi.1008330.g003

We fit the model with three latent factors, whose inferred activity timeseries were sparse (Fig 3C). Including additional latent factors beyond these resulted in better models of the SA of individual neurons or small subsets of neurons, but caused little improvement in the overall quality of model fit for this fish (Fig 3D). To understand the relative importance of each factor to the overall model fit, we defined a contribution index for a factor as the average reduction in the quality of model fit following its deletion (Methods). We found that each factor made a substantial contribution to the overall model fit by modulating shared SA across large groups of tectal neurons (Fig 3E and 3F).

The factor coupling matrix (defined by the vectors b_n) reports how neurons are affected by the latent factors. There are several possibilities for how this matrix could be structured. First, if there is a minimal presence of structured SA the coupling matrix may exhibit no coherent organisation at all. Second, neurons could require the coordinated activity of several latent factors to explain their SA. This would be the case if, e.g., neurons participated in multiple recurrently connected circuits driven by noise [28, 29], and would result in factors modulating overlapping groups of neurons. Third, latent factors may each drive their own distinct sets of neurons, with little cross-talk between them. This could occur if, e.g., latent factors were encoding unrelated streams of motor or non-visual sensory information [2, 30]. In our example zebrafish the estimated coupling matrix had a highly modular structure, with factors influencing largely non-overlapping sets of neurons (Fig 3E). Furthermore, the factor cross-correlograms showed no sign of dependence between factors, indicating that distinct sets of neurons were uniquely targeted by independent latent sources of SA (Fig 3G).

Since our model fits receptive fields jointly with latent sources of SA, the estimated tuning curve for each neuron already accounts for ongoing SA that may have inflated its responses to stimuli. Indeed, if spontaneous calcium transients coincide with the presentation of a stimulus, one could expect tuning curves obtained by simply averaging the fluorescence levels over a small window following stimulus onset to be spuriously larger, and exhibit higher variance than if these events did not occur. We plotted the tuning curves estimated by CILVA against tuning curves obtained by averaging (see Methods) and found that they confirmed this intuition (Fig 3H). Moreover, sorting the neurons according to their preferred stimulus revealed a more refined retinotopic map when explicitly accounting for SA (Fig 3I). While our data showed a variety of tuning types, of note were neurons that were unselective to visual stimuli (i.e., had relatively flat tuning curves), but that were highly active throughout the recording (e.g., the neuron marked by an asterisk in Fig 3H and S8 Fig).

Neurons are differentially driven by external stimuli and latent internal factors

To quantify the extent to which each neuron is driven by sensory stimuli versus shared SA, we derived an equation that expressed the variance of the reconstructed fluorescence levels in terms of three components: the variance attributable solely to EA, the variance attributable solely to shared SA, and the covariance (i.e., interaction) between EA and shared SA (Fig 4A, Methods). This revealed that across the population there was a continuous progression of responses, with some neurons being primarily driven by EA, some primarily by SA, and some by a mixture of both SA and EA (Fig 4A). To confirm that these effects were not artefacts of the model or the calcium indicator, we verified that the model does not overestimate the variance in the data (Fig 4B) and that there were interactions between EA and SA that were greater than expected by chance (Fig 4C). We defined a “drive ratio” to measure whether neurons were driven more by SA or EA (Methods). The distribution of drive ratios was largely bimodal (Fig 4D), indicating a preference to be dominated by either EA or SA rather than responding equally to both.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Analysis of the contribution of EA and SA to neural variability.

All data for the same fish as in Fig 3. (A) Top: composition of each neuron’s sample variance in terms of variance attributable solely to EA (red bars), solely to shared SA (blue bars), and their covariance (orange bars). Orange bars represent absolute values of covariances for ease of visualisation. Variance components are given as proportions of the total sample variance of the raw fluorescence signal var[f_n] (corrected for imaging noise, see Methods). Neurons sorted by the strength of their coupling to each factor (as in Fig 3E). Bottom: coupling between neurons and latent sources of SA suggests neurons with strong coupling are weakly driven by sensory stimuli. Maximum bar height of one. (B) Sample variance (corrected for imaging noise) vs variance of the statistical model indicates that the model does not overestimate variance. Each data point represents one neuron. (C) Covariances between evoked and spontaneous traces estimated by the model (vertical axis). Chance levels for a null model (horizontal axis) are 95th percentiles of shuffled data obtained by cyclically permuting evoked traces by random offsets 1000 times while preserving temporal structure. Sample covariances exceeding chance levels (orange circles above dashed identity line) cannot be attributed to the slow timescale of the calcium indicator. (D) Distribution of drive ratios across the population of neurons. (E) Correlation coefficient between raw fluorescence trace and evoked component of model fit (without SA) and full model fit (with SA). Neurons with strongly negative drive ratios show marked improvement in the quality of model fit. (F) Violin plots showing statistically significant improvement in the average correlation coefficient between experimental data and model fits after incorporating latent sources of SA (p < 0.001, Wilcoxon signed-rank test). (G) Spatial organisation of latent factors underlying SA. The three non-overlapping factors are spatially localised and tile the imaging plane. (H) Spatial organisation of the evoked variance components. Cell opacity is proportional to the fraction of variance attributable to EA for the given neuron. Neurons strongly driven by EA cluster in the middle region of the tectum. (I) Same as H, but for SA. Neurons strongly driven by SA cluster in the anterior and posterior tectum.

https://doi.org/10.1371/journal.pcbi.1008330.g004

We next quantified the improvement that resulted from incorporating SA into the model. Without SA, the model for each neuron consists of a simple linear filter convolved with a calcium kernel. This is a good description of neurons that possess high drive ratios (i.e., whose variances are dominated by EA) and thus these neurons show little improvement in how well the statistical model fits their fluorescence levels with the incorporation of SA (Fig 4E, neurons along the diagonal). In contrast, many neurons are poorly fit by a model that incorporates only stimulus responses, and show substantial improvement when shared sources of SA are taken into account (Fig 4E, neurons above the diagonal), leading to a significant increase in the average correlation between fluorescence traces and model fits (Fig 4F).

In the absence of sensory stimulation, SA in the optic tectum has previously been shown to exhibit a characteristic localised spatial structure [31]. We thus sought to determine whether this effect persisted when the tectum was being actively driven by sensory stimulation. The factors underlying SA identified by CILVA concentrated in the posterior, middle, and anterior regions of the tectum, together tiling the two dimensional imaging plane (Fig 4G). Interestingly, the evoked variance component was largely confined to the middle tectum, where coupling to latent factors was weakest (Fig 4H). Conversely, the spontaneous variance component was most strongly represented at the posterior and anterior ends of the tectum, where coupling to latent factors was strongest, with little SA in the middle tectum (Fig 4I). This spatial localisation was not imposed on the data by our model, and thus is a useful post-hoc verification that the SA our model identifies is likely to be biologically salient.

To determine how representative our results were, we fit the model to a dataset of seven additional zebrafish larvae. Example fits for two of these zebrafish are given in S9 and S10 Figs. For consistency of comparison between zebrafish we again fit the statistical model with three latent factors. Across the 8 fish the mean correlation coefficient was centered at ∼ 0.6 (S11A Fig), with latent factors having average individual contribution indices of 0.1 (S11B Fig). Factors were also mostly non-overlapping, with only a small fraction of neurons participating in multiple factors (S11C Fig). Incorporating all three factors increased the mean correlation coefficient between raw fluorescence data and model-fit by 34% on average compared to model-fits without the SA component (S11D Fig). Finally, we found that while EA and SA tended to be balanced at the population level, individual neurons mostly biased their activity towards being either stimulus-driven or spontaneous (S11E and S11F Fig). These results show that the basic statistical properties of the data are consistent across a set of different animals.

CILVA identifies low dimensional patterns of SA in visual cortex

We next explored the application of the model to publicly available data from mouse primary visual cortex [32]. In this case, stimuli of higher dimension were presented more rapidly than in our previous application. Briefly, head-fixed mice expressing the calcium indicator GCaMP6s (via viral injection) stood on an air-suspended ball while drifting gratings were presented across the visual field with 1 to 3 second intervals and at 8 orientations, 3 spatial frequencies, and 4 temporal frequencies (Fig 5A). We verified with simulated data that the model could accurately recover evoked and spontaneous components in this regime (S3 Fig, S2 Table), and then applied CILVA to decouple the evoked and spontaneous fluorescence components (Fig 5B–5D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Single-trial decoupling of EA and SA in visual cortex.

(A) Calcium imaging of mouse V1 during presentation of drifting gratings. (B) Raw data consists of 21 minutes of neural activity from 986 neurons. The fluorescence trace of each neuron is normalised to take values between 0 (dark) and 1 (light). Neurons sorted as in panel F. (C) Decoupled evoked component of neural activity. (D) Decoupled spontaneous component of neural activity, given a latent dimensionality of 10. (E) Distribution of correlation coefficients between data and model fits. Shuffled data obtained by cyclically permuting each trace by a random offset while preserving its temporal structure. (F) Learned factor coupling matrix showing that the inferred factors target largely non-overlapping sets of neurons.

https://doi.org/10.1371/journal.pcbi.1008330.g005

CILVA was able to extract low dimensional patterns of SA (Fig 5D, vertical bands of activity) that were much harder to discern in the raw data (Fig 5B). This included a spontaneous event that appeared to be triggered by stimulus onset (Fig 5D, first vertical band of activity) but that did not reoccur with subsequent stimulus presentations. The model reconstruction provided a good fit, with correlation coefficients much larger than in the case of shuffled data (Fig 5E). Similar to the zebrafish data, extracted latent factors were mutually independent (S12 Fig), targeted largely non-overlapping sets of neurons (Fig 5F), and were sparsely active (S12 Fig). CILVA is thus effective for discovering novel, interpretable patterns of neural activity in high dimensional cortical imaging data.

Zebrafish recordings

All procedures were performed with approval from The University of Queensland Animal Ethics Committee (approval certificate number QBI/152/16/ARC). Nacre zebrafish (Danio rerio) embryos expressing elavl3:H2B-GCaMP6s, of either sex, were collected and raised according to established procedures [40] and kept under a 14/10 hr on/off light cycle.

Zebrafish larvae were embedded in 2.5% low-melting point agarose, positioned at the centre of a 35 mm diameter plastic petri dish and overlaid with E3 embryo medium. Calcium imaging was performed at a depth of 70 μm from the dorsal surface of the tectal midline. Time-lapse two-photon images were acquired using a Zeiss LSM 710 inverted two-photon microscope. A custom-made inverter tube composed of a pair of beam-steering mirrors and two identical 60 mm focal length lenses arranged in a 4f configuration was used to allow imaging with a 40X/1.0 NA water-dipping objective (Zeiss) in an upright configuration. Samples were excited via a Spectra-Physics Mai TaiDeepSee Ti:Sapphire laser (Spectra-Physics) at an excitation wavelength of 940 nm and the emitted light was bandpass filtered (500–550 nm). Laser power at the sample ranged between 12 to 20 mW. Images of 416x300 pixels were obtained at 2.1646 Hz. To improve the stability of the recording, chambers were allowed to settle for three hours prior to start of two-photon imaging.

Visual stimuli were projected on white paper placed around the wall of a 35 mm diameter petri dish using a projector (PK320 Optoma, USA), covering a horizontal field of view of 174°. A red filter (Zeiss LP590 filter) was placed in the front of the projector to avoid interference of the projected image in the signal collected by the detector. Larvae were aligned with one eye facing the white paper side of the dish and with the body axis orthogonal to the projector. Visual stimuli were generated using custom software based on MATLAB (MathWorks) and Psychophysics Toolbox. Each trial consisted of 6° diameter black spots at nine different positions, separated by 15° intervals from 45° to 165°, where 0° was defined as the direction of the larva’s body axis. Their order was set to maximise spatial separation within a trial (45°, 120°, 60°, 135°, 75°, 150°, 90°, 165°, 105°). Spots were presented for 1 s, followed by 19 s of blank screen. We projected consecutive trials of nine spots with 25 s of inter-trial interval.

The cell segmentation procedure is described in ref. [41]. Briefly, custom MATLAB software was used to automatically detect the region-of-interest (ROI) of each active cell, i.e., the group of pixels defining each cell. The software searched for active pixels, i.e., pixels that showed changes in brightness across frames, resulting in an activity heatmap of all the active regions across frames. The activity map was then segmented into regions using a watershed algorithm, with a similar threshold applied to all movies. Within each segmented region, we computed correlation coefficients of all pixels in the region with the mean of the most active pixel and its eight neighbouring pixels. Correlation coefficients showed a bimodal distribution; one peak of highly correlated pixels representing pixels of the cell within the region, and a second peak of relatively low correlation coefficients representing nearby pixels within the region which were not part of the cell. Using a Gaussian mixture model, we found the threshold correlation which differentiated between pixels likely to form the active cell and neighbouring pixels that were not part of the cell. We also required that each detected active area covered at least 26 pixels (5.5 mm²). The software allowed visual inspection and modification of the parameter values where needed. All pixels assigned to a given cell were averaged to give a raw fluorescence trace over time.

Mouse recordings

We used publicly available data [32]. The experimental procedures are described in detail in refs. [42, 43]. Neurons were recorded simultaneously at 2.5Hz using calcium imaging and segmented using Suite2p. Visual stimuli were shown at approximately 1 Hz, with randomized inter-stimulus intervals. The stimuli were drifting gratings with 8 directions, 4 spatial frequencies and 3 temporal frequencies. Blank stimuli (gray screen) were also interleaved. While this data was originally recorded from multiple imaging planes, to avoid issues with timing of latent factor activity between planes we restricted our analysis to neurons from a single imaging plane, resulting in a population of 986 neurons.

Residual NMF method

As described in the main text, we first considered using a combination of non-negative least squares and non-negative matrix factorisation to decouple EA and SA (S1 Fig). We defined stimulus regressors for i = 1, …, K by convolving the binary stimulus timeseries with a calcium impulse response kernel k (a difference-of-exponentials function, defined below), giving ϕ_i = k ∗ s_i. We then estimated regression coefficients β_ni by solving the non-negative minimisation problem where β_n = (β_n1, …, β_nK)^⊤. The evoked component of the fluorescence signal can be defined in terms of the estimated regression coefficients as where . Here NNLS refers to the non-negative least squares algorithm used to perform the minimisation. Residual data e_n was then defined as where σ is a linear rectifier σ(x) = max(0, x) applied elementwise ensuring non-negativity of the residuals. We then applied NMF to the residual data by solving the minimisation problem where and . The L rows of are timeseries describing the evolution of low dimensional structure, and the columns of describe how neurons are coupled to such timeseries. The NMF-estimated SA for neuron n is given by the projection of the latent timeseries onto a single dimension by the corresponding row of The full fluorescence data can thus be approximately reconstructed as We use the NNLS routine in SciPy and the NMF routine in scikit-learn [44].

To evaluate the models of SA produced by NMF we then expressed the spontaneous components in terms of a basis of calcium impulses by solving the non-negative minimisation problem where each ψ_t is a calcium response following a unit impulse at time t, Here is a vector of length T that takes the value of 1 at t and 0 elsewhere. The expression of in the basis of calcium transients is then given by , where . Note that this differs from the basis of stimulus regressors used to model the stimulus-driven component of neural activity as we employ a regressor for each time point t in the entire trace.

CILVA model

Model fitting

We fit the model by maximising the posterior density of the latent variables, (7) where the parameters of the model are θ = ({α_n}, {β_n}, {w_n}, {b_n}). Ideally, one could perform this optimisation using the expectation-maximisation algorithm, which alternates between computing the posterior distribution over the latent factors q(x) = p(x|f, θ, γ), and maximising the posterior expectation . However, the E-step is not analytically tractable since our exponential prior on x_l(t) is non-conjugate for the likelihood model. Instead, we use a related “pseudo expectation-maximisation” approach [45] that alternately optimises Eq 7 according to the steps (8) (9) until numerical convergence or until i reaches a user-specified number of iterations. The alternating maximisations are each performed using the bounded BFGS algorithm with limited memory (L-BFGS-B), with exact gradients derived below.

The logarithm of the joint model probability density is where the constant term does not depend on the parameters of θ to be estimated. Let ℓ(x, θ) = ln p(f, x|θ, γ), and let denote the model reconstruction error for neuron n in imaging frame t, The derivatives of ℓ(x, θ) with respect to the parameters and latent variables are then where for a matrix the convolution is performed row-wise. In practice we vectorise the computation of the gradients to improve efficiency.

The imaging noise variance terms are estimated using the method in ref. [23]. Specifically, is estimated as the mean of the power spectral density of f_n over the range (F_s/4, F_s/2), where F_s is the imaging rate of the fluorescence microscope.

Simulated data

To generate simulated data we sampled the latent factors from a zero-inflated exponential distribution with probability ξ of a non-zero latent event, This ensured the latent factor activity was sparse. We also introduced a private SA term z_n(t) for neuron n at time t by sampling from a zero-inflated exponential with probability π of a non-zero private event, The intensity function was then given by with the fluorescence levels following the standard CILVA model with a common imaging noise variance σ², We sampled α_n from the discrete uniform distribution on {2, …, 10} and for simplicity set β_n = 0. The tuning curves w_n were defined as Gaussian functions x ↦ exp(−(x − μ_n)²/2ν). For the simulation of data in response to well-spaced, low dimensional stimuli (cf. Fig 3) we sampled the centres μ_n uniformly from the interval [0, K], where K is the number of stimuli, and sampled the widths ν_n uniformly from [0, K/2]. For the simulation of data in response to rapidly presented, high dimensional stimuli (cf. Fig 5) we chose our receptive fields to be more selective and sampled ν_n uniformly from .

Factor coupling vectors b_n were defined by evenly assigning the N neurons to L factors, and sampling b_nl ∼ U[q, 1] if neuron n is assigned to factor l, and b_nl ∼ U[0, 1 − q] otherwise. We found q = 0.85 provided simulations that appeared similar to the experimental data. We characterised the model reconstruction quality in terms of π and σ² in S2 and S3 Figs, with the associated model parameters provided in S1 and S2 Tables.