Imperial College LondonDepartment of Electrical and Electronic EngineeringDetecting Cells and Cellular Activity fromTwo-Photon Calcium Imaging DataStephanie ReynoldsSupervised by Pier Luigi Dragotti and Simon SchultzSubmitted in part fulfilment of the requirements for the degree ofDoctor of Philosophy of Imperial College London2Declaration of OriginalityI declare that this thesis and the research it contains are the product of my own workunder the guidance of my thesis supervisors, Professors Pier Luigi Dragotti and SimonSchultz. All material that is not my own work has been properly acknowledged.Stephanie Reynolds34Copyright DeclarationThe copyright of this thesis rests with the author and is made available under a CreativeCommons Attribution Non-Commercial No Derivatives licence. Researchers are freeto copy, distribute or transmit the thesis on the condition that they attribute it, thatthey do not use it for commercial purposes and that they do not alter, transform orbuild upon it. For any reuse or redistribution, researchers must make clear to othersthe licence terms of this work.56AbstractTo understand how networks of neurons process information, it is essential to monitortheir activity in living tissue. Information is transmitted between neurons by elec-trochemical impulses called action potentials or spikes. Calcium-sensitive fluorescentprobes, which emit a characteristic pulse of fluorescence in response to a spike, areused to visualise spiking activity. Combined with two-photon microscopy, they enablethe spiking activity of thousands of neurons to be monitored simultaneously at single-cell and single-spike resolution. In this thesis, we develop signal processing tools fordetecting cells and cellular activity from two-photon calcium imaging data.Firstly, we present a method to detect the locations of cells within a video. In ourframework, an active contour evolves guided by a model-based cost function to identifya cell boundary. We demonstrate that this method, which includes no assumptionsabout typical cell shape or temporal activity, is able to detect cells with varied prop-erties from real imaging data.Once the location of a cell has been identified, its spiking activity must be inferredfrom the fluorescence signal. We present a metric that quantifies the similarity betweeninferred spikes and the ground truth. The proposed metric assesses the similarity ofpulse trains obtained from convolution of the spike trains with a smoothing pulse,whose width is derived from the statistics of the data. We demonstrate that the pro-posed metric is more sensitive than existing metrics to the temporal and rate precisionof inferred spike trains.Finally, we extend an existing framework for spike inference to accommodate a widerclass of fluorescence signals. Our method, which is based on finite rate of innovationtheory, exploits the known parametric structure of the signal to infer the unknownspike times. On in vitro imaging data, we demonstrate that the updated algorithmoutperforms a state of the art approach.78AcknowledgementsFirst and foremost, I would like to thank my primary supervisor, Pier Luigi Dragotti.When I first decided to pursue a PhD, countless older and wiser people told me thatthe most important factor is finding the right supervisor and, four years later, I findthat I agree. I am grateful for his clear guidance, the many opportunities he gave meto teach, learn and attend conferences, and his understanding. I would like to thankmy second supervisor Simon Schultz, whose infectious enthusiasm made him great funto work alongside. I am grateful for the all the support and advice he has given meduring this time.I would like to thank our collaborators Per Jesper Sjöström, Therese Abrahamsson andRenaud Schuck, who collected data that was essential to this thesis. I have learnt alot from interactions with members of the Schultz lab, thanks in particular to AmandaFoust, Caroline Copeland and Peter Quicke. I am also thankful to my examiners,Timothy Constandinou and Dimitri Van De Ville for their valuable input regardingimprovements to my thesis.I am grateful to my friends in the C&SP group for brightening up the more difficultPhD days. I will miss the fiery lunchtime debates, office gossips, coffees and cocktails.Thanks to you all for being such great company. Special thanks go to Tricia, for beingmy organising partner and for keeping the more rowdy members of the group in check;to John, whose big smile and words of wisdom made him an awesome desk neighbour;to Alex and Michael, who shared the journey with me; and Charlie, Doris, Ilia, Sithanand Wilhelm, for the good times.Finally, I would like to express my gratitude to my parents for their endless support.To my friends and family for reminding me what is important. And to Leo, for listeningto every presentation, sharing both my difficulties and my successes, and for makingme smile every day.910Contents1 Introduction 231.1 Original contribution and outline of thesis . . . . . . . . . . . . . . . . . 241.2 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252 Background 272.1 Neuronal information processing . . . . . . . . . . . . . . . . . . . . . . 272.1.1 The action potential . . . . . . . . . . . . . . . . . . . . . . . . . 282.2 Monitoring neuronal activity at cellular resolution . . . . . . . . . . . . 292.3 Calcium imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.1 Signal model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.4 Two-photon microscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.4.1 Fluorescence excitation . . . . . . . . . . . . . . . . . . . . . . . 342.4.2 Scattering and absorption in biological tissue . . . . . . . . . . . 362.4.3 Fluorescence collection . . . . . . . . . . . . . . . . . . . . . . . . 362.4.4 Shot noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.4.5 Neuropil contamination . . . . . . . . . . . . . . . . . . . . . . . 392.5 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40113 ABLE: an activity-based level set segmentation algorithm 433.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.2 Active contours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.2.1 External energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3 Level set method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.3.1 External velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.3.2 Regularisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4 Proposed method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.4.1 Dissimilarity metrics . . . . . . . . . . . . . . . . . . . . . . . . . 553.4.2 Neighbouring cells . . . . . . . . . . . . . . . . . . . . . . . . . . 553.4.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.5 Experimental methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.5.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.5.2 Two-photon calcium imaging of hippocampal slices . . . . . . . . 603.5.3 In vivo imaging datasets . . . . . . . . . . . . . . . . . . . . . . . 613.5.4 Performance metrics . . . . . . . . . . . . . . . . . . . . . . . . . 613.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.6.1 Robustness to heterogeneity in cell properties . . . . . . . . . . . 623.6.2 Demixing overlapping cells . . . . . . . . . . . . . . . . . . . . . 633.6.3 Results on hippocampal brain slices . . . . . . . . . . . . . . . . 653.6.4 Algorithm comparison on manually labelled dataset . . . . . . . 653.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67124 CosMIC: A consistent metric for spike inference from calcium imag-ing 694.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 694.2 Fuzzy set theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.2.1 Sørensen-Dice coefficient . . . . . . . . . . . . . . . . . . . . . . . 744.3 Constructing the metric . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3.1 Ancestor metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.4 Temporal error tolerance . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.4.1 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . 804.4.2 Cramér-Rao bound . . . . . . . . . . . . . . . . . . . . . . . . . . 814.4.3 Bound for spike inference . . . . . . . . . . . . . . . . . . . . . . 824.4.4 Pulse width . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.5 Experimental methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.5.1 Success rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.5.2 Spike train correlation . . . . . . . . . . . . . . . . . . . . . . . . 874.5.3 Two-photon calcium imaging of cortical slices . . . . . . . . . . . 874.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.6.1 CosMIC penalises overestimation . . . . . . . . . . . . . . . . . . 884.6.2 CosMIC rewards high temporal precision . . . . . . . . . . . . . 904.6.3 Application to real imaging data . . . . . . . . . . . . . . . . . . 924.6.4 CosMIC discriminates precision and fall-out rate of spike trains . 944.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95135 Spike inference using Finite Rate of Innovation theory 975.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.2 Sampling theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.2.1 Signals with finite rate of innovation . . . . . . . . . . . . . . . . 1015.2.2 Sampling kernels . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.2.3 Reconstruction of a stream of Diracs . . . . . . . . . . . . . . . . 1045.3 Spike inference from calcium imaging data . . . . . . . . . . . . . . . . . 1065.3.1 Mapping samples to the FRI setting . . . . . . . . . . . . . . . . 1085.3.2 Amplitude estimation . . . . . . . . . . . . . . . . . . . . . . . . 1115.4 Spike inference from noisy data . . . . . . . . . . . . . . . . . . . . . . . 1125.4.1 Pre-whitening . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1135.4.2 Model order estimation . . . . . . . . . . . . . . . . . . . . . . . 1155.4.3 Spike inference algorithm . . . . . . . . . . . . . . . . . . . . . . 1165.5 Experimental methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1175.5.1 Simulated data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1175.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.6.1 Comparison with instantaneous-rise algorithm . . . . . . . . . . 1185.6.2 Model order estimation . . . . . . . . . . . . . . . . . . . . . . . 1195.6.3 Comparison with state of the art . . . . . . . . . . . . . . . . . . 1215.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1226 Conclusion 1236.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1236.2 Future research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12414A ABLE: supplemental analytical results 127A.1 External velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127A.2 Euclidean dissimilarity metric . . . . . . . . . . . . . . . . . . . . . . . . 129B CosMIC: proof of analytical results 131Bibliography 13515List of Figures2.1 Example of a whole-cell patch-clamp recording. . . . . . . . . . . . . . . 292.2 Model of the neuronal fluorescence signal. . . . . . . . . . . . . . . . . . 322.3 Procedure of fluorescence excitation by a two-photon microscope. . . . . 352.4 Example of in vivo two-photon calcium imaging data. . . . . . . . . . . 372.5 Illustration of neuropil contamination of cellular signal. . . . . . . . . . 403.1 Example of an active contour segmenting a piecewise smooth image. . . 473.2 The level set function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.3 The Chan-Vese active contours model. . . . . . . . . . . . . . . . . . . . 523.4 Flow diagram of the proposed segmentation algorithm. . . . . . . . . . . 543.5 Evolution of multiple active contours. . . . . . . . . . . . . . . . . . . . 573.6 Segmentation results on in vivo imaging data. . . . . . . . . . . . . . . . 623.7 Segmentation results on simulated data containing overlapping cells. . . 633.8 Segmentation results on in vitro imaging data containing synchronouslyspiking cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.9 Algorithm comparison on a manually labelled dataset. . . . . . . . . . . 664.1 Inconsistency in usage of spike inference metrics. . . . . . . . . . . . . . 704.2 Membership functions of the union and intersection of fuzzy sets. . . . . 73164.3 Flow diagram of the proposed metric. . . . . . . . . . . . . . . . . . . . 774.4 Proposed metric example. . . . . . . . . . . . . . . . . . . . . . . . . . . 784.5 Metric pulse width. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.6 Illustration of common metrics . . . . . . . . . . . . . . . . . . . . . . . 864.7 Results on spike train estimates with varying overestimation ratio. . . . 894.8 Results on spike train estimates with varying temporal precision. . . . . 914.9 Results on real imaging data . . . . . . . . . . . . . . . . . . . . . . . . 934.10 Results on spike train estimates with varying precision and fall-out rate. 965.1 Example of a slow-rise fluorescence signal. . . . . . . . . . . . . . . . . . 995.2 Example of signals with finite rate of innovation. . . . . . . . . . . . . . 1025.3 Exponential reproduction with an E-spline of order 2. . . . . . . . . . . 1035.4 Example of sampling and reconstruction of a stream of Diracs. . . . . . 1055.5 Example of spike inference with the FRI algorithm. . . . . . . . . . . . . 1095.6 Model order estimation procedure. . . . . . . . . . . . . . . . . . . . . . 1145.7 Example of simulated data. . . . . . . . . . . . . . . . . . . . . . . . . . 1175.8 Comparison of the temporal error of the FRI algorithm with instantaneous-rise and slow-rise models. . . . . . . . . . . . . . . . . . . . . . . . . . . 1195.9 Accuracy of the model order estimation procedure on simulated data. . 1205.10 Algorithm comparison on real data . . . . . . . . . . . . . . . . . . . . . 12217List of Tables2.1 Time constants of four commonly used calcium indicators. . . . . . . . . 343.1 Runtime of the segmentation algorithm (minutes) on synthetic data. . . 583.2 Algorithm success rate on a manually labelled dataset. . . . . . . . . . . 675.1 Amplitude of simulated calcium transients. . . . . . . . . . . . . . . . . 11818NomenclatureAcronyms2PCI Two-photon calcium imagingAP Action potentialCRB Cramér-Rao boundFRI Finite rate of innovationGECI Genetically encoded calcium indicatori.i.d. Independent identically distributedICC Intracellular calcium concentrationLSF Level set functionMSE Mean square errorPDF Probability density functionRMSE Root mean square errorROI Region of interestSDC Sørensen-Dice coefficientSNR Signal-to-noise ratioSTC Spike train correlationNumber SetsR Real numbers19R>0 Positive real numbersZ IntegersC Complex numbers: x+ iy for x, y ∈ RVariablesg Vector: g ∈ Cm with m > 1G Matrix: G ∈ Cm×n with m,n > 1g(t) Continuous-time signalg[n] Discrete-time signalFunctionsg : A→ B A function, g, defined on domain A with image a subset of B〈g1(t), g2(t)〉 Inner product:∫R g1(t)g∗2(t) dtg1(t) ∗ g2(t) Continuous-time convolution:∫R g1(τ)g2(t− τ) dτg1[n] ∗ g2[n] Discrete-time convolution:∑m∈Z g1[m]g2[n−m]1A(t) Indicator function, equal to one if t ∈ A and 0 otherwise1t>a Indicator function of set A = {t ∈ R : t > a}δ(t) Dirac delta functionδ[n] Discrete Dirac delta functionGH Hermitian (conjugate transpose) of matrix GH(x) Heaviside functionProbabilityE [U ] Expectation of random variable UP(A) Probability of event AN (µ, σ2) Normal distribution with mean µ and standard deviation σNorms20|a| Absolute value of a real or complex number a‖g‖ Euclidean norm of vector g:√g[1]2 + g[2]2 + ...+ g[m]2‖g(t)‖1 L1-norm of continuous-time function g(t):∫R|g(t)| dt|A|c Cardinality of a fuzzy setAsymptotic notationO(h(x)) Referred to as ‘Big-O’ notation, this indicates the set of functions, suchthat{g(x) : ∃ c, x0 ∈ R>0 such that 0 ≤ g(x) < ch(x) ∀x ≥ x0}2122Chapter 1.IntroductionThe human brain is a remarkably complex organ responsible for a vast array of tasks.It is believed to consist of approximately 100 billion neurons, whose activity is organ-ised in hierarchical networks. Information is sent between neurons by electrochemicalimpulses called action potentials (APs) or spikes. Coordinated spiking activity acrossmultiple spatial scales allows the brain to, for example, process incoming sensory in-formation, encode memories and initiate motor activities.To understand how networks of neurons process information, it is essential to observebrain activity in living tissue. Crucially, this can provide insight into the differencebetween healthy and unhealthy brains. For example, functional magnetic resonanceimaging (fMRI) and electroencephalography (EEG) have been successfully used toidentify impairment in connectivity in patients with brain disorders [17]. Such findingscan have clinical relevance both as a diagnostic tool and as a focus for new treatments.While fMRI and EEG are effective for assessing larger-scale connectivity, the study ofmicrocircuit activity has been limited by the challenges of recording at cellular resolu-tion in living tissue. Traditionally, electrophysiological approaches have been used tomonitor a cell’s spiking activity. These approaches are, however, limited to samplingfrom a small number of unidentifiable cells. As an alternative, fluorescent sensors canbe used to visualise cellular activity via changes in their fluorescence intensity. Onesuch approach is calcium imaging, whereby calcium-sensitive probes monitor the levelsof intracellular calcium, which, in turn, reflects a cell’s spiking activity.To perform calcium imaging, a light source is required that excites fluorescence fromthe probes. In lightsheet microscopy (LSM), for example, fluorescence is excited from23Chapter 1. Introductionall neurons in a thin sheet of tissue simultaneously [1, 12]. This technology can beused to monitor neuronal activity at high speeds in transparent organisms, such as thelarval zebrafish. When imaging through light scattering media, however, the spatialresolution of LSM suffers. In contrast, two-photon microscopy overcomes light scatter-ing through the usage of a highly localised focal spot that is scanned across the tissueto generate an image. Due to its relatively high spatial resolution, this technology isunrivalled in its ability to image hundreds of microns deep in living tissue. State of theart two photon calcium imaging is able to monitor the activity of thousands of cellssimultaneously in living animals [124, 100]. The prospects for this technology are thusvast; it has already been used, for example, to study Alzheimer’s disease in animalmodels [112] and to determine the spatial organisation of receptive fields [86, 10].Before neuronal activity from two-photon calcium imaging (2PCI) videos can be anal-ysed, two challenging signal processing problems must be tackled. Firstly, the locationsof cells within the video must be identified. This task is complicated by heterogene-ity in cellular properties, such as their shape and typical temporal activity, and thefrequency of cellular overlap. Secondly, one must recover a cell’s spiking activity fromits fluorescence signal. This problem, which we refer to as spike inference, is an in-verse problem that can be solved by considering the known parametric structure ofthe fluorescence signal. In this thesis, we focus on these problems and outline ourcontributions below.1.1 Original contribution and outline of thesisThe remainder of the thesis proceeds as follows.In Chapter 2, we provide relevant background information on 2PCI data. We discusshow the data is generated and outline the sources of noise.In Chapter 3, we present a novel segmentation algorithm for detecting the locationof cells from 2PCI data. In our framework, multiple coupled active contours evolve,guided by a model based cost function, to identify cell boundaries. On in vivo imag-ing data, we demonstrate the versatility of the algorithm, which includes no priors ona cell’s stereotypical morphology or temporal activity. We also demonstrate the ro-bustness of the algorithm on challenging datasets exhibiting highly correlated cellularactivity and overlapping cells, respectively.241.2. PublicationsIn Chapter 4, we consider the metrics that are used to compare the locations of inferredspikes with ground truth spiking activity. We highlight the limitations of existingmetrics and propose an alternative. Rather than operating on the true and estimatedspike trains directly, the proposed metric assesses the similarity of the pulse trainsobtained from convolution of the spike trains with a smoothing pulse, whose widthis derived from the statistics of the imaging data. On real and simulated data, wedemonstrate the sensitivity of the proposed metric to the temporal and rate precisionof inferred spike trains.In previous work, it was established that the neuronal fluorescence signals that areinduced by 2PCI belong to the class of signals with finite rate of innovation (FRI).Using FRI theory and a simplified model of the fluorescence signal, an algorithm wasdeveloped that inferred the locations of spikes with high accuracy [87]. In Chapter 5, weextend this framework to encompass a broader class of fluorescence signals, which arenot well-matched by the original framework. We also introduce procedures to improvethe algorithm’s robustness to noise. Using the spike inference metric presented inChapter 4, we demonstrate that the proposed algorithm outperforms a state of the artapproach on in vitro imaging data.Finally, in Chapter 6 we summarise the contributions of the thesis and discuss somepossible directions for future research.1.2 PublicationsThe materials presented in this thesis have led to the following publications:Peer-reviewed journalsS. Reynolds, T. Abrahamsson, P. J. Sjöström, S. R. Schultz, and P. L. Dragotti.CosMIC: A consistent metric for spike inference from calcium imaging, Neural Com-putation, in press.S. Reynolds, T. Abrahamsson, R. Schuck, P. J. Sjöström, S. R. Schultz, and P. L.Dragotti. ABLE: An activity-based level set segmentation algorithm for two-photoncalcium imaging data, eNeuro, 4(5):ENEURO.0012-17.2017, 2017.25Chapter 1. IntroductionR. Schuck, M. A. Go, S. Garasto, S. Reynolds, P. L. Dragotti, and S. R. Schultz.Multiphoton minimal inertia scanning for fast acquisition of neural activity signals,Journal of Neural Engineering , 15(2):025003, 2018.Peer-reviewed conferencesS. Reynolds, C. S. Copeland, S. R. Schultz, and P. L. Dragotti. An extension of theFRI framework for calcium transient detection. In 13th International Symposium onBiomedical Imaging (ISBI), pages 676 – 679, 2016. IEEE.S. Reynolds, J. Oñativia, C. S. Copeland, S. R. Schultz, and P. L. Dragotti. Spikedetection using FRI methods and protein calcium sensors: performance analysis andcomparisons. In 11th international conference on Sampling Theory and Applications(SampTA), pages 533 – 537, 2015. IEEE.26Chapter 2.BackgroundTo understand function and dysfunction in the mammalian brain it is crucial to ob-serve its components in action. In particular, long-term studies of brain activity inliving animals (‘in vivo’) could be key to identifying how the brain encodes stimulior orchestrates behaviours. Whilst there are many viable techniques for observingbrain activity at larger spatial scales, two-photon calcium imaging is unparalleled inits ability to monitor cellular activity, deep in living tissue, at single-cell resolution.To emphasise the importance of this technology, in Section 2.1, we discuss some ba-sic aspects of neuroscience. In Section 2.2, we discuss the relative merits of existingtechniques to monitor neuronal activity at single-cell resolution. Finally, in Sections2.3 and 2.4, we provide further detail on calcium imaging and two-photon microscopy,respectively.2.1 Neuronal information processingThe mammalian cortex is responsible for a wide array of tasks, such as motor control,integration of sensory inputs and decision making [151]. It is believed that commu-nication occurs across multiple spatial scales, with local networks of interconnectingneurons constituting microcircuits whose activity is aggregated in larger macrocircuitsspread across multiple brain areas. Microcircuits, which consist of 1000 - 10,000 neu-rons, are hypothesised to be computational units that undertake specific informationprocessing tasks [81]. For example, one ‘canonical’ microcircuit has been observedto repeat throughout the cortex and is presumed to be responsible for fundamental27Chapter 2. Backgroundcomputations [37]. Despite this apparent hierarchical organisation, local networks areheavily interconnected with other brain areas.To understand information processing it is, therefore, necessary to monitor networkactivity at multiple spatial scales. At the macroscopic scale, technologies such asfMRI can be used to study coordinated activity across brain areas. At the microscopicscale, cellular resolution functional imaging has the potential to elucidate the dynamicsof microcircuits. The action potential (AP), an electro-chemical impulse propagatedbetween connected neurons, is the core method of information transmission betweenneurons [70]. It is not yet fully understood how neuronal networks encode informationthrough coordinated spiking activity. Accordingly, many technologies are aimed atmonitoring neuronal population activity at single-cell and single-spike resolution.2.1.1 The action potentialInputs from other neurons or sensory organs can excite a neuron and cause it to firean electro-chemical impulse called an action potential. The resting potential of a neu-ron, which reflects the difference in ion concentrations inside and outside the cell, isapproximately -70mV. Inputs from other neurons arrive in the form of inward currentinjections, causing the neuron’s membrane potential to become less negative. Thisis referred to as depolarisation. If the stimulus is sufficiently large to depolarise theneuron above a (neuron-dependent) threshold, an AP is fired. Over a period of ap-proximately 1ms, the membrane potential rapidly rises to +30mV, subsequently fallingback to the resting potential. This rapid rise in membrane potential is referred to asan AP, or spike. An AP is thought to be an all-or-nothing event: sub-threshold inputsdo not illicit an AP and all supra-threshold inputs illicit APs of the same amplitude.So, although sub-threshold inputs can be important in identifying connections betweenneurons, it is thought that APs are key to neural computation. A neuron’s activitycan therefore be characterised by the set of times at which that neuron fired a spike:{tk}Kk=1 with all tk ∈ R. This activity can be represented as a spike trainx(t) =K∑k=1δ(t− tk), (2.1)where δ(·) is the Dirac delta function.Individual neurons encode information through the rate and timing of the spikes they282.2. Monitoring neuronal activity at cellular resolutionFigure 2.1: Example of a whole-cell patch-clamp recording. This recording,from an online dataset [132, 26], was sampled at 10kHz. Action potentials are visibleas peaks in the membrane potential.fire [70]. For example, the firing rate of motor neurons is related to the force ofcontraction of associated muscle fibres. In contrast, visual cortex neurons encodeinformation about visual stimuli using precise timing of individual spikes [66]. In thecat visual system, temporal coding has been observed on the millisecond time-scale. Todecipher the neural code, it is therefore necessary to record from the brain at sufficienttemporal resolution to detect individual spikes.2.2 Monitoring neuronal activity at cellular resolutionTo enable analysis of neuronal microcircuit activity, recording techniques aim to max-imise the number of neurons that can be monitored at single-cell and single-spikeresolution. The patch-clamp technique is the current gold standard for monitoringthe voltage of a single cell [50, 115] and is commonly used as a benchmark for newtechnologies. To obtain a patch-clamp recording, an electrode is inserted into braintissue and attached to a cell membrane. This configuration allows the cell’s voltageto be recorded at a sampling rate on the order of 10kHz, which is sufficient to resolvesingle APs (Fig. 2.1). As one electrode is required per cell, the patch-clamp techniqueis limited by the number of neurons that can be monitored simultaneously. Althoughstate of the art implementations can record from up to 8 neurons at the same time[16, 150], this is not sufficient for analyses of population activity.In contrast to the patch-clamp technique, which monitors the intracellular voltagesignal, extracellular approaches take advantage of the fact that an AP can be detectedfrom voltage fluctuations outside of the cell body. Instead of pairing one electrode29Chapter 2. Backgroundwith one cell, these techniques insert electrodes with multiple recording sites intobrain tissue [18]. Each recording site detects transmembrane potentials, which areindicative of APs, from several cells. Activity must then be assigned to the cell fromwhich it originated by a process known as spike sorting. At each site, spikes may stemfrom neurons as far as 150-200µm away (a neuron’s cell body has diameter 10µm),although typically only activity from cells within 50µm will have sufficient amplitudeto be detected [51]. Despite their relatively high temporal resolution and signal-to-noise ratio (SNR), extracellular electrophysiological approaches are not well-suited tolong-term studies, which are necessary to understand learning, memory and plasticityin the living brain. This limitation is due to the instability in the monitored neuronalpopulations. In one study, statistical analysis of spike waveforms and distribution wasused to show that the number of units recording consistently from the same neuronalpopulations fell to 57% over one week [33]. In contrast, fluorescence imaging, whichcan provide both functional and anatomical information, can be used to monitor thesame neuronal populations over weeks and months [98].Whereas an electrode must be in contact with the tissue that it monitors, fluorescenceimaging can be performed non-invasively to monitor cellular activity in vivo. In par-ticular, functional imaging utilises sensors in the brain whose fluorescence intensityvaries in response to cellular activity. One option is to use voltage-sensitive fluorescentindicators [101], which visualise changes in a cell’s membrane potential. Due to thefast time scale of APs (on the order of 1ms), to visualise cellular activity with a volt-age sensor it is necessary to have sub-millisecond temporal resolution, which, in turn,limits the signal-to-noise ratio of the acquired signals. This limitation means that,in practice, voltage imaging is mostly used to perform widefield imaging of cellularaggregates [71]. Alternatively, a secondary variable can be imaged, from which spikingactivity can be inferred. In this thesis, we focus on calcium imaging data, which usesfluorescent sensors to visualise a cell’s intracellular calcium concentration (ICC). AnAP causes the ICC to rise to levels that are 10 to 100 times the levels at baseline [9],this high dynamic range makes the ICC a favourable signal to monitor. The durationof the ICC elevation after an AP is an order of magnitude slower than the durationof an AP. As a consequence, cells can be monitored with calcium imaging at lowersampling rates whilst still capturing high-frequency spiking activity.Fluorescence imaging approaches require a light source, with which to excite fluores-cence from the sample (the biological tissue), and a detector, with which to collect thefluorescence. Approaches that excite fluorescence from a large volume simultaneously,302.3. Calcium imagingsuch as confocal [97] and lightsheet [1] microscopy, treat excited photons differentlydepending on how they arrive at the detector. Confocal microscopes reject all photonsoutside of a detection pinhole, thereby ideally only collecting photons from a localisedregion. In contrast, lightsheet microscopes use detector arrays, such as CMOS cam-eras. Each detector represents a pixel in the output image and, in theory, collectsexcited photons from a localised region of the sample. While appropriate for trans-parent specimens or superficial layers of tissue, these approaches are not well-suitedto imaging through light-scattering media. Changes in the refractive index, which canoccur every nanometre in the mammalian brain [151], divert light from its path. Thisphenomenon is referred to as scattering (see Section 2.4.2). It is hard to determine thesource in the sample from which scattered photons originate, which complicates theuse of detection pinholes and camera arrays. Two-photon microscopy, on the otherhand, uses a highly localised laser to excite fluorescence [31]. This allows all photonsemitted from the sample to be collected by a single detector, regardless of the angleat which they arrive at the detector, as they reliably stem from a small volume.2.3 Calcium imagingCalcium ions are a crucial component in many aspects of cellular function [9]. Forexample, they trigger the contraction of heart muscle cells and regulate both cell deathand proliferation. In neurons, the ICC is tightly linked to spiking activity. When anAP is fired in a neuron, the depolarisation causes voltage-gated calcium ion channelsto open. As a consequence, there is a brief influx and a subsequent slower efflux ofcalcium ions. This transient elevation in the ICC is referred to as a calcium transient.Compared to the rapid timescale of an AP, a calcium transient is relatively slow; thecalcium ion influx lasts approximately 1ms [11] and the efflux resembles an exponentialdecay with time constant 0.5s - 1s [46].The development of calcium-sensitive fluorescent indicators (‘calcium indicators’) hasenabled the visualisation of the ICC and, indirectly, neuronal spiking activity. Theindicators consist of a fluorophore and a domain to which calcium ions can bind. Cal-cium binding alters the properties of the indicator and, in turn, elicits a stereotypicalchange in the fluorescence, the nature of which varies depending on the indicator used[14]. Calcium indicators can be categorised into two groups: synthetic indicators andgenetically encoded indicators (GECIs). The former group, which has been in devel-31Chapter 2. BackgroundFigure 2.2: Model of the neuronal fluorescence signal. The fluorescenceemitted from a neuron can be modelled as that neuron’s spike train convolved with acharacteristic pulse. The rise and decay rate of the pulse varies with the fluorescentindicator used to generate the imaging data (F). We display simulated fluorescencetraces generated from the spiking activity in A and the characteristic pulse shapes ofthe indicators GCaMP6s (B), GCaMP6f (C), OGB-1 (D) and Cal-520 (E).opment since the early 1980s [138], can be used to label large cell populations in intacttissue through bolus injection [130]. They are, however, limited to acute imaging ex-periments (on the order of hours) [71]. GECIs, on the other hand, can be encoded intothe DNA of transgenic animals [80, 2, 26]. In that case, they are naturally expressed inneurons and can be targeted to specific cell types and subcellular compartments. Thismakes GECIs ideal for chronic imaging experiments over multiple weeks and months[55]. Moreover, recent advances have produced GECIs whose sensitivity and signalquality surpasses those of synthetic indicators [26], which was not the case as recentlyas 2009 [136].In the following, we introduce a mathematical model relating the fluorescence intensityof somatic calcium indicators to spiking activity. The mathematical model introducedbelow is prevalent in the literature, in particular it is commonly used by algorithmsthat infer spike timing from calcium imaging data [148, 147, 30, 73, 87, 109, 105]. For amore comprehensive coverage of the biophysical processes underlying the mathematicalmodel, the reader is referred to [14].322.3. Calcium imaging2.3.1 Signal modelWhen an AP is fired from a neuron, there is a rapid influx and subsequent slower effluxof calcium ions, which is referred to as a calcium transient. This calcium response ismodelled with a rapid rise and an exponential decay. A corresponding fluorescencepulse is emitted from the calcium-sensitive fluorescent indicators in the neuron. Thepulse emitted as a result of an AP at time t = 0 is proportional top(t) = c(1− e−t/τon)e−t/τoff 1t>0, (2.2)where τon and τoff are parameters which define the rate of the rise and decay, respec-tively, and c is a normalisation constant that ensures that p(t) has peak value equal to1.The rise time is predominantly influenced by the rate at which ions bind to the fluores-cent indicator (the association rate). The rate at which the fluorescence pulse decaysis affected by cellular properties, such as the extrusion rate of calcium channels, andproperties of the indicator, including the association and disassociation rates. In Table2.1, we indicate values of τon and τoff obtained from the literature for the followingindicators: Cal-520 AM [133], OGB-1 AM [73], GCaMP6f and GCaMP6s [26]. Onnew data generated by those indicators, these values can be used as a guideline. Thecorresponding characteristic pulses are plotted in Fig. 2.2F.Each time a neuron fires an AP, a stereotypical pulse of fluorescence is emitted. Writ-ing ak as the amplitude of the pulse emitted in response to an AP at time tk, thefluorescence signal, f : R→ R, isf(t) = b(t) +(K∑k=1akδ(t− tk))∗ p(t) (2.3)= b(t) + cK∑k=1ak(1− e−(t−tk)/τon)e−(t−tk)/τoff 1t>tk , (2.4)where b(t) reflects the baseline fluorescence. The baseline component is typically mod-elled as either a low-frequency drift [30, 106] or a piecewise constant function [44, 104].This can reflect ICC changes in response to other cellular processes [48] or a steadydecline in baseline fluorescence due to bleaching of the fluorescent indicator.The amplitude of calcium transients varies with the local spike rate of neurons. The33Chapter 2. BackgroundFluorescent indicator τon (ms) τoff (ms)GCaMP6f 18 205GCaMP6s 72 794OGB-1 10 667Cal-520 32 314Table 2.1: Time constants of four commonly used calcium indicators. Thepulse shape is defined by two parameters, τon and τoff, which represent the characteristicdelays. These values can be used as a guideline for each calcium indicator; in practice,they will vary with the indicator expression level as well as the cell type.co-operative calcium binding of some GECIs [14] results in increasing amplitudes asthe spike rate increases [2, 26]. In contrast, synthetic indicators are known to saturateat high spike rates [64], producing spikes of decreasing amplitude. The amplitude ofspikes can also steadily decrease as a result of indicator bleaching. In Fig. 2.2, wedisplay simulated fluorescence signals generated by the same spiking activity for fourcommon calcium indicators.2.4 Two-photon microscopyTwo-photon microscopy is a powerful technique for imaging deep in light-scatteringmedia. A beam of light is used to excite fluorescence from a concentrated focal spot inthe sample, this process is discussed in Section 2.4.1. The effectiveness of the excitationof fluorescence in biological tissues is affected by scattering and absorption, which isdiscussed in Section 2.4.2. Finally, in Section 2.4.3, the acquisition process of a two-photon microscope image is described and sources of noise are discussed in Sections2.4.4 and 2.4.5.2.4.1 Fluorescence excitationFluorescence microscopy relies on the use of fluorophores (fluorescent molecules) thatcan be excited by the absorption of photons. When a fluorophore absorbs a photon ofsufficient energy, it transitions from a ground state to an excited state. The fluorophoresubsequently relaxes to the ground state by either a radiative or non-radiative process.In the latter case, light energy is absorbed as thermal energy. In the former, a photon342.4. Two-photon microscopyA BFigure 2.3: Procedure of fluorescence excitation by a two-photon micro-scope. A beam of light is focused through the microscope objective such thatfluorescence is only excited in the vicinity of the focal spot (A). The focal spot isscanned sequentially through the field of view to generate one frame of an image (B).The beam moves discretely along the scan path with increment ∆x, spending a fixedtime at each point.is emitted by the fluorophore. This process, which occurs on the nanosecond timescale,is known as fluorescence.In standard fluorescence microscopy, a single photon is absorbed by the fluorophore,some energy is transformed into thermal energy and a lower-energy photon is subse-quently emitted. In two-photon microscopy, excitation only occurs when two photonshave been absorbed essentially simultaneously. With the combined energies of the twoabsorbed photons, the fluorophore transitions into an excited state. To relax fromthis state, some energy is absorbed as thermal energy and a higher-energy photon isemitted. Typically, the excitation beam has wavelength in the near-infrared range andthe light emitted from the fluorophore is in the visible spectral range [151].A two-photon excitation event can only occur when two photons arrive at the samemolecule within approximately 0.5 femtoseconds of each other [52]. To enable thisevent, the photons in the excitation light beam must be highly concentrated bothspatially and temporally. To concentrate photons temporally, the excitation lightbeams consist of ultrashort pulses (less than a picosecond in length) with high peakintensities [52]. Spatial density is achieved by focusing the excitation light through amicroscope objective with high numerical aperture. As a result, fluorescence is onlyexcited in a highly localised focal spot (Fig. 2.3A). Indeed, the probability of excitationscales with the square of the light intensity, meaning that the excitation decays sharplyboth axially and laterally.35Chapter 2. Background2.4.2 Scattering and absorption in biological tissueThe scattering and absorption of light in biological tissue have a large impact on thesignal generated in fluorescence microscopy. Biological tissues contain numerous en-dogenous molecules that absorb light of different wavelengths [151]. Before photonscan reach the focal spot, they may be absorbed by these molecules. Although thisabsorption does attenuate the light at the focal spot, the greater impact stems fromthe associated tissue damage. This damage, which limits the viability of confocalmicroscopy for long-term imaging of biological tissue, is mitigated in two-photon mi-croscopy [128]. In the latter case, the localised focal spot spares the majority of thetissue from high-intensity light. Moreover, the excitation beam uses near-infrared light,which is absorbed by few endogenous molecules [131].A greater obstacle is presented by the scattering of light by biological tissue. Lightis deflected from its original path as the refractive index of a medium changes. Thisphenomenon is referred to as scattering. Before light reaches its focal spot, it passesthrough multiple components, such as cells and blood vessels, each of which has adifferent refractive index. Even within these components, the refractive index maychange on the nanometer scale [151]. Consequently, the intensity of non-scatteredlight decreases exponentially with penetration depth [52]. This, in turn, reduces theintensity of fluorescence excited at the focal spot, with greater effect as the imagingdepth increases. In vivo fluorescence microscopy is thus fundamentally limited byscattering. As longer wavelength light is affected less in this regard, the near-infraredbeam used in two-photon microscopy mitigates some of this effect. The maximumimaging depth, which depends on the sparseness and brightness of fluorophores, ishard to quantify. Two-photon microscopy is, however, sufficiently robust to scatteringthat it is routinely used to image hundreds of microns deep in tissue [59].2.4.3 Fluorescence collectionAs the spatial density of scattered photons is generally low, it can be assumed thatthey do not excite significant fluorescence. The majority of emitted photons, therefore,can be assumed to originate from the vicinity of the focal spot. Consequently, a two-photon microscope collects all emitted photons, regardless of the angle from whichthey appear to originate [52].To generate an image, the focal spot is typically scanned in a raster pattern across362.4. Two-photon microscopyFigure 2.4: Example of in-vivo two-photon calcium imaging data. Theexcitation laser scans the field of view once to produce a single frame of the video(A). In the mean image, which displays the mean of each pixel’s activity, bright spotscorrespond to cells with high baseline fluorescence (B). In C, we display the correlationimage, which, at each pixel, is computed as the mean of the correlation coefficientsbetween that pixel’s temporal activity and the activity of neighbouring pixels. Brightspots in the correlation image correspond to active cells. In D, E and F, we focus onthe three regions highlighted in the mean and correlation images. The regions visiblein D, E and F correspond to a cell, blood vessel and background region, respectively.In each subfigure, we display the zoomed in mean and correlation images. We alsoplot the average temporal activity from all pixels in the interior of the object and theactivity from a single pixel within the object. The data presented in this figure wascollected by Svoboda et al. [99].the field of view (Fig. 2.3B). The scanning is performed discretely; after dwelling ∆tseconds in one location, the laser either moves ∆x microns in the x-direction or startsa new row. Generally the field of view is square; the between-row spacing is equal tothe within-row (column) spacing and there are the same number of rows as columns.A single image is generated by completing the scan path once. A video of N frames is37Chapter 2. Backgroundgenerated by completing the scan path N times, each completed scan path generatesthe data for one frame of the video (Fig. 2.4A).The spatial increment, ∆x, must be sufficiently small that single cells can be resolved.Meanwhile, the area covered — (M∆x)2 with M the number of rows/columns — mustbe sufficiently large to contain enough cells for a meaningful analysis. To maximisethe number of monitored cells, therefore, one would wish to maximise M . Similarly,the dwell-time at each spatial location, ∆t, must be long enough to obtain a usableSNR (see Section 2.4.4). However, the sampling period (the temporal spacing betweensamples at the same spatial location), which is approximately equal to M2∆t, mustbe sufficiently small that individual APs can be observed. To minimise the samplingperiod, therefore, the dwell-time or number of monitored cells must be compromised.Typically, ∆x and ∆t are in the range 1-5µm and 10-100ms, respectively.The trade-off between temporal resolution and other factors is lessened by the useof calcium-sensitive fluorescent probes, as calcium transients are significantly slowerthan the APs that induce them. To further loosen the constraints on the data, signalprocessing methods are required that are robust in low-noise settings. This will reducethe necessary dwell-time and further increase the number of cells that can be monitoredat single-cell and single-spike resolution.2.4.4 Shot noiseIn the following, we introduce a mathematical model for the noise introduced duringthe acquisition process. This model relates to the statistical fluctuations in photoncount and electron emission and is standard in the literature [97, 43, 27]. The photonsemitted from the sample are detected by a photomultiplier tube, which is extremelysensitive to low levels of light. On detecting photons, a photocathode emits electrons,whose number is subsequently multiplied by a sequence of dynodes [97]. The pixelintensity value is obtained by analogue-to-digital conversion of current emitted fromthe final dynode.As photons arrive independently at the photocathode, the number arriving during theacquisition period is modelled as a Poisson random variable. The number of electronsemitted in response to these photons, E0, can be shown to be Poisson distributed [27]with mean dependent on the quantum efficiency of the photocathode and the rate ofphoton arrival. At each subsequent dynode, the number of electrons is multiplied by382.4. Two-photon microscopyan amount dependent on an amplification factor, δ. The emission processes instigatedby different incident electrons are physically independent [27]. As such, the numberof electrons emitted from the ith dynode, Ei, is a sum of Ei−1 independent randomvariables, which model the output in response to the electrons emitted by the previousdynode. Given ND dynodes in total, we haveEi =Ei−1∑m=1ηi,m, (2.5)for i ∈ {1, ..., ND} and ηi,m ∼ Poisson(δ).It has been shown empirically that, when the rate of photon influx is not too low,the distribution of END can be approximated by a Gaussian, whose parameters aredependent upon the rate of photon influx, the amplification factor of the dynodes andthe quantum efficiency of the photocathode. As the rate of photon influx varies at eachpixel depending on the level of activity at the source, the statistics of the noise are notconstant spatially. Rather, the noise is approximated as a heteroskedastic Gaussian.The signal at a pixel x is writtenĨ(x) = I(x) + σ(I(x))ξ, (2.6)where I(x) is the underlying, deterministic signal and ξ ∼ N (0, 1). Assuming that theanalogue-to-digital conversion process is linear, we have [27]σ(I(x)) = aI(x) + b, (2.7)i.e. the noise level at a pixel increases linearly with the level of activity at that pixel.In practice, this model is simplified by assuming that σ is constant at each pixel[147, 105, 106, 44, 127].2.4.5 Neuropil contaminationAs discussed in Section 2.4.3, in two-photon microscopy, fluorescence is excited froma localised focal spot within the sample. Photons collected at one point on the scanpath stem from a small volume centred around this focal spot. The size and shapeof this volume are characterised by the point spread function of the microscope [52].Particularly when imaging at depth, as scattering of light becomes an increasing prob-lem, the resolution of this focal spot suffers. At best, in-plane (lateral) resolution is39Chapter 2. BackgroundABCDEFigure 2.5: Illustration of neuropil contamination of cellular signal. Theaverage temporal signal from pixels inside the cell delineated by the dark orange con-tour in A are plotted in C. The neighbouring pixels that do not belong to another cell(highlighted in B) exhibit neuropil activity (D), which originates from densely packedcellular processes. Due to the relatively low axial resolution of the focal volume, thisactivity is present in the cellular signal. Subtracting the weighted neuropil signal fromthe raw cellular signal removes this artefact (E). The data presented in this figure wascollected by Svoboda et al. [99].sub-micron and axial (perpendicular to the imaging plane) resolution is on the orderof several microns [51].The effect of the relatively low axial resolution is that cellular signals are contaminatedwith signals from neighbouring objects. In the cortex, where cells are relatively sparselydistributed, the main confounding factor is neuropil contamination [46]. Neuropilsignal, which corresponds to the average activity of densely packed cellular processes[46], is present in ‘background’ pixels that do not contain cell bodies. It has beenfound to correlate strongly with the electrocortigram [64] and be tuned to behaviouralactivity [100]. Pixels in cell somata are often contaminated with neuropil activity (Fig.2.5). To decontaminate cellular signals, it is common to subtract the weighted localneuropil signal from the cellular signal. The value of the weight parameter chosenvaries depending on the imaging set-up [100, 26, 63].2.5 OutlookTwo-photon microscopy is a powerful tool for imaging with cellular resolution severalhundred microns deep in living tissue. New generations of calcium indicators have402.5. Outlooksufficient dynamic range to visualise single APs. Moreover, genetically encoded in-dicators can be expressed chronically, through non-invasive methods, in specific cellpopulations. Paired together, these technologies allow the same neuronal populationsto be monitored at single-cell and single-spike resolution over periods of weeks andmonths [98]. They have been used to study the spatial organisation of receptive fields[86, 10] and to monitor changes in neural activity during learning [55, 100].At present, cellular resolution imaging tends to occur over areas less than 1mm2 on asingle plane; even state of the art microscopes are limited to imaging thousands of cellssimultaneously [124]. To allow more complex behaviours to be studied, the expansionof a microscope’s field of view is the focus of constant development [59]. The trade-offbetween SNR and cell number is a limiting factor in this regard — a larger field ofview necessitates a compromise on dwell-time and, therefore, noise level at each pixel.To further push the boundary of this technology, robust signal processing methods arerequired, which can maximise the information gleaned from noisy data and, in turn,relax the constraints on the imaging systems.41Chapter 2. Background42Chapter 3.ABLE: an activity-based level set seg-mentation algorithmA state of the art two-photon calcium imaging video may contain activity from thou-sands of neurons. Before neuronal activity can be analysed, the location of each cellwithin the video must be identified. This task is complicated by the heterogeneityin cellular properties, both spatial and temporal. For example, the calcium indicatorused to generate a video affects a cell’s stereotypical temporal activity and its apparentshape — some genetically encoded indicators are excluded from the nucleus and there-fore produce fluorescent ‘doughnuts’. Segmentation algorithms are further impeded bythe frequency of cellular overlap, which occurs due to the relatively low axial resolu-tion of the two-photon microscope. These challenges necessitate flexible segmentationalgorithms with minimal assumptions on the properties of regions of interest (ROIs).In this chapter, we present an algorithm for detecting the location of cells from two-photon calcium imaging data. In our framework, multiple coupled active contoursevolve, guided by a model-based cost function, to identify cell boundaries. An activecontour seeks to partition a local region into two subregions, a cell interior and ex-terior, in which all pixels have maximally ‘similar’ time courses. This simple, localmodel allows contours to be evolved predominantly independently. When contours aresufficiently close, their evolution is coupled, in a manner that permits overlap. Weillustrate the ability of the proposed method to identify the boundaries of overlappingcells. We also demonstrate the versatility of our framework, which includes no priorson a cell’s morphology or stereotypical temporal activity, on mouse in vivo imagingThe work presented in this chapter led to the following publication [107].43Chapter 3. ABLE: an activity-based level set segmentation algorithmdata.3.1 IntroductionA broad spectrum of methods have been proposed to overcome the challenges inherentin segmenting calcium imaging data. One confounding factor in the development ofsuch algorithms is the heterogeneity in neuroscientists’ selection criteria for ROIs.While some are only interested in somata [86], others also require the detection ofdendrites [120]. While some require only active cells [78], others analyse the proportionof cells that are active [100]. A set of such selection criteria, which can be hard to definemathematically, can instead be implicitly incorporated through the use of supervisedlearning from manually labelled data [143, 4, 65]. Whereas an earlier approach useda non-binary decision tree on pixel features, such as mean, variance and covariance[143], later approaches have utilised convolutional neural networks (CNNs) [4, 65]. OneCNN has been found to outperform most state of the art unsupervised methods on abenchmark dataset [65]. However, an obstacle for all supervised approaches is that thestatistics of imaging videos vary significantly with factors such as the brain area, themicroscope zoom and objective. In this case, the benchmark test dataset contained no‘unseen’ data — all test samples had a corresponding training sample acquired underthe same imaging conditions. It is, therefore, not yet clear how well supervised learningapproaches generalise to new data from different imaging conditions.A second subset of approaches bypass the (2+1)-D imaging video, utilising insteada 2-D summary image [122, 92, 60, 65]. This approach is practical; by significantlyreducing the dimension of the input data, it allows videos to be processed that cannotbe stored in the RAM of a typical laboratory workstation. Typically, one of twosummary images is used: the mean image [92], which displays the mean of each pixel’sactivity, and the correlation image [122, 60, 65], which displays the average correlationbetween the activity of a pixel and that of its neighbours. While efficient, the usageof a 2D summary image limits the types of cells that can be detected. Only activecells are visible in the correlation image, whereas sparsely firing cells have low visibilityin the mean image. Indeed, when a CNN was trained on manually labelled (2+1)-Dimaging videos rather than the corresponding mean images, it was found to detectmore transiently active neurons [4]. A further impediment to such approaches is theprevalence of cells in 2D summary images that appear to overlap. To identify the true443.1. Introductionboundaries of overlapping cells it is necessary to consider the (2+1)-D imaging data.The (2+1)-D imaging video, which consists of two spatial dimensions and one temporaldimension, is often prohibitively large to work on directly. One family of approachestherefore reshapes the (2+1)-D imaging video into a 2-D matrix. The resulting matrixadmits a decomposition, derived from a generative model of the imaging video, into twomatrices, each encoding spatial and temporal information. Early approaches appliedindependent component analysis to find the spatial components that have maximallystatistically independent temporal activity [118, 82]. The constraint of temporal inde-pendence, however, causes algorithm performance to deteriorate in the case of overlap-ping cells [105, 154]. By instead including a constraint on non-negativity of the matrixentries, Maruyama et al. estimated the component matrices using non-negative ma-trix factorization (NMF) [76]. Pnevmatikakis et al. further developed this approachinto a widely-used constrained NMF (CNMF) algorithm [105, 45]. The additionalconstraints on neuronal temporal activity, which is assumed to be an auto-regressiveprocess derived from calcium dynamics, and sparsity constraints on the size of spa-tial components, have been shown to improve the performance of CNMF over plainNMF [105]. Related methods suggest performance gains can be achieved by includinga more accurate model of the background activity [94, 56] or learning a dictionary ofcellular features from the data [34]. By expressing the (2+1)-D imaging video as a 2-Dmatrix, this type of approach can achieve relatively high processing speeds. This does,however, come at the cost of discarded spatial information, which occurs when the twospatial dimensions are reduced to one. To mitigate this, either morphological filters[82, 94, 105, 45] or manual curation [118, 76] are required in post-processing steps.In this chapter, we present a segmentation method in which cell boundaries are detectedby multiple coupled active contours. To evolve an active contour we use the level setmethod, which is a popular tool in bioimaging due to its topological flexibility [28]. Toeach active contour, we associate a higher-dimensional function, referred to as a levelset function, whose zero level set encodes the contour location. We implicitly evolvean active contour via the level set function, whose evolution is driven by a local modelof the imaging data temporal activity. The model includes no assumptions on a cell’smorphology or stereotypical temporal activity. We demonstrate that this versatilityallows us to segment a variety of cell types and shapes and both active and inactivecells. For convenience, we refer to our method as ABLE (an Activity-Based LEvel setmethod). In Sections 3.2 and 3.3, we discuss relevant background on active contoursand the level set method, respectively. Then, in Section 3.4, we present the proposed45Chapter 3. ABLE: an activity-based level set segmentation algorithmsegmentation algorithm. Finally, in Sections 3.5 and 3.6, we describe the data on whichwe evaluate ABLE and assess ABLE’s segmentation performance, respectively.3.2 Active contoursThe goal of image segmentation is to partition the domain of an image, Ω, into multiplesubregions. When the boundary of a single object is sought, the domain is partitionedinto two disjoint subregions: Ωin, the interior, and Ωout, the background. The objectboundary, δΩin, is then obtained from the closed curve delineating the border of Ωin(Fig. 3.1A). Although in our application the video has three dimensions, we only wishto partition the spatial domain of the video, which is two-dimensional. This is due tothe assumption that the objects are static during an imaging experiment. Otherwise,it would be necessary to partition the whole spatiotemporal domain [79, 110]. Theimage is denoted with a function, I : Ω → Rm, that associates a vector of length mwith every pixel, x, in the image domain, Ω ⊂ R2. When m = 1, we have a single-frameimage, when m > 1, we have a video consisting of m frames.An active contour is a commonly-used tool for identifying an object’s boundary. Start-ing from a seed location, the contour evolves towards the object boundary, see Fig. 3.1.Its evolution is guided by a cost function that is designed to be minimised when thecontour is located at the true boundary. We denote the active contour, which evolvesin artificial time, τ ∈ R, with Γ(τ). At each timestep, Γ(τ) represents the most recentestimate of δΩin. Active contours are versatile; they can accommodate a broad range ofcost functions and evolution frameworks (for a review of the bioimaging literature see[28]). Consequently, since their conception in 1988 [61], they have found widespreadusage in image processing, both in segmentation and, for example, in motion detection[58] and multi-view scene modelling [8].There are two main approaches to modelling an active contour: parametric and im-plicit. In the former approach, a curve is represented by a continuous parameter, adiscrete set of coefficients and a method of interpolation from the discrete coefficientsto the continuous curve [61, 57]. For example, curves can efficiently be constructedas linear combinations of B-spline basis functions1 [42, 15]. While effective for objectsthat have stereotypical shapes, this approach is not ideal for irregularly-shaped objects1A B-spline is a type of piecewise polynomial function (a ‘spline’) with minimal support for a givendegree of smoothness.463.2. Active contoursFigure 3.1: Example of an active contour segmenting a piecewise smoothimage. The image consists of an object, Ωin, and a background, Ωout (A). An activecontour, Γ, evolves guided by a cost function that is minimised when it is located atthe true object boundary, δΩin (B - D). After Nmax iterations, the contour convergesat the true object boundary.whose shape cannot easily be represented by a combination of basis functions.The second class of approaches evolves the contour implicitly via a higher-dimensionalfunction [21, 90]. In this framework, the active contour is defined as the zero level setof a manifold, φ : Ω× R→ R, that also evolves in time, such thatΓ(τ) = {x ∈ Ω : φ(x, τ) = 0}. (3.1)Unless shape constraints are explicitly included in the cost function, this class of ap-proaches allows a greater degree of topological flexibility than parametric active con-tours. A further advantage of this framework is that it can easily be extended to imageswith more than two spatial dimensions. This is fitting for calcium imaging data, whichis starting to be collected in three spatial dimensions [1]. This is the framework weadopt for our segmentation algorithm, which we introduce in more detail in Section3.3.The evolution of an active contour is guided by a cost function, which is designedto be minimised when the contour is located at the true object boundary. The costfunction is typically formed of two components: one relating to the image, which isreferred to as the external energy, and another relating to the properties of the contouritself, which we refer to as the regulariser. Both components are weighted by scalarparameters, µ, λ ∈ R, to form the cost function:E(Γ) = λ Eext(Γ) + µ Ereg(Γ), (3.2)47Chapter 3. ABLE: an activity-based level set segmentation algorithmwhere Eext is the external energy and Ereg is the regulariser. The external energy, whichwe discuss in the following section, incorporates assumptions about the properties ofthe object to be detected. The regulariser, which we discuss in Section 3.3.2, encouragesthe object boundary to conform to a set of pre-defined conditions. For example, jaggededges or excessively long boundaries may be penalised.3.2.1 External energyTo construct the external energy, quantitative image features are computed. Thefeatures are typically one of two types: boundary-based and region-based features.As objects in calcium imaging videos tend not to possess sharp edges (due to theblurring induced by the point spread function), we focus on the latter approach. Twofeatures f in ∈ Rm and fout ∈ Rm are derived from values of the image in the interiorand exterior of Γ, respectively. These features, which are re-computed each time thecontour is updated, are then used to classify pixels into the interior or exterior regions.Using a function, D : Rm × Rm → R≥0, that computes the discrepancy between theimage at pixel x and a feature, we haveEext(Γ) = λin∫inside ΓD(I(x), f in)dx + λout∫outside ΓD(I(x), fout)dx, (3.3)where λin and λout are real-valued weight parameters. The external energy thus pe-nalises discrepancies between the value of an image at pixel x and the feature of thesubregion to which it belongs. The relative impact of discrepancies in the interior andexterior is controlled by the weights, λin and λout.Chan and Vese proposed a region-based active contours approach for piecewise smooth,single-frame images containing a single object and a background [22, 24]. The externalenergy, which is inspired by the Mumford-Shah functional for the segmentation ofpiecewise smooth images [83], is equal toEext(Γ) = λin∫inside Γ‖I(x)− f in‖2dx + λout∫outside Γ‖I(x)− fout‖2dx, (3.4)where f in and fout are the average intensity of pixels in the interior and exterior,respectively. This energy penalises discrepancies between a pixel’s intensity and theaverage intensity of the subregion to which it belongs. The initial formulation wasdesigned for single-frame images, where I(x), f in, fout ∈ R, and later extended to vector-valued images [23] and images of more than two pieces [25].483.3. Level set methodFigure 3.2: The level set function (LSF). The active contour (A) is encodedas the zero level set of the LSF (B). The LSF is positive inside the contour, negativeoutside and zero on the contour.3.3 Level set methodTo update the active contour we use the level set method of Osher and Sethian [91].This method was first introduced to image processing by [20] and [74]; it has sincefound widespread use in the field. We represent the active contour by a function φ,which is positive for all x inside the contour, negative for x outside the contour andequal to zero for x on the boundary (see Fig. 3.2). We refer to φ as the level setfunction (LSF), as its zero level set encodes the contour of interest. We note that wehave reversed the convention of negative φ inside the contour and positive φ outsidein order to improve the clarity of figures representing the LSF.We aim to find the active contour that minimises the energy in Eq. (3.2). We thusseek the LSF that minimises the analogous cost function:E(φ) = λ Eext(φ) + µ Ereg(φ), (3.5)where we have replaced Γ with φ. As it is in general not possible to minimise this costfunction directly, a standard way to obtain the LSF that minimises Eq. (3.5) is to findthe steady-state solution to the gradient flow equation [5]:∂φ∂τ= −∂E∂φ. (3.6)As artificial time, denoted by τ , advances, the LSF is thus evolved in the image planein a manner that yields the sharpest decrease in the cost function. From Eq. (3.5),49Chapter 3. ABLE: an activity-based level set segmentation algorithmthis partial differential equation (PDE) becomes∂φ∂τ= −(λ∂Eext∂φ+ µ∂Ereg∂φ). (3.7)The evolution of the LSF is thus driven by two components: the former providesthe impetus from the image data, and is referred to as the external velocity, and thelatter the impetus from the regulariser. In Sections 3.3.1 and 3.3.2, we discuss thecomputation of each component.To solve the PDE numerically, it is necessary to discretise Eq. (3.7). We approximateφ by storing and updating its value on the uniformly-spaced image grid. The evolutionparameter τ is also discretised, such thatφ(τ + ∆τ) = φ(τ)−∆τ(λ∂Eext∂φ+ µ∂Ereg∂φ). (3.8)At every timestep τ , the LSF is updated until a maximum number of iterations, Nmax,has been reached or the active contour has converged. Using one or both of theseconditions is common in the active contour literature, see, for example, [29, 69]. Weretain µ∆τ < 0.25 in order to satisfy the Courant-Friedrichs-Lewy condition [69] — anecessary condition for the convergence of a numerically-solved PDE. This conditionrequires that the numerical waves propagate at least as fast as the physical waves [89].3.3.1 External velocityTo compute the external velocity, we rewrite the external energy from Eq. (3.3) interms of the Heaviside function, such thatEext(φ) =∫ΩAin(x)H(φ(x)) +Aout(x)(1−H(φ(x)))dx, (3.9)where Ain(x) = λinD(I(x), f in) and Aout(x) = λoutD(I(x), fout). Here, we are im-plicitly using H(φ) and (1−H(φ)) as indicator functions of the subregions inside andoutside the active contour, respectively. To calculate the derivative with respect toφ, we use the standard variational formulation [23, 90, 5]. The Heaviside function is503.3. Level set methodapproximated by a smooth function, Hε(x), whereHε(x) =1 x > ε12(1 + xε +1π sin(πxε))|x| ≤ ε0 x < −ε.(3.10)The derivative of the Heaviside function, the delta function, is approximated by δε(x),whereδε(x) =12ε(1 + cos(πxε))|x| ≤ ε0 |x| > ε.(3.11)We note that δε(x) =ddxHε(x). Replacing the Heaviside function in Eq. (3.9) with itsapproximation and computing the derivative, the external velocity at x ∈ Ω is∂Eext∂φ(x) = δε(φ(x))(Ain(x)−Aout(x))(3.12)= δε(φ(x))(λinD(I(x), f in)− λoutD(I(x), fout)). (3.13)The term δε(φ(x)), which is only non-zero at pixels on or near the active contour,acts as a localization operator, ensuring that the velocity only impacts the LSF in thevicinity of the active contour. The parameter ε determines the radius of the non-zeroband.An example of the external velocity for the single-image Chan-Vese active contourmodel is shown in Fig. 3.3. From Eq. (3.4), setting λin = λout = 1, we have∂Eext∂φ(x) = δε(φ(x))(‖I(x)− f in‖2 − ‖I(x)− fout‖2), (3.14)where, in this case, f in and fout are the average intensity values within each subregion.The contour is thus moved towards the pixels whose values are more similar to f inthan fout, and repelled by other pixels.3.3.2 RegularisationThe external velocity updates the LSF only in the vicinity of the zero level set. As theLSF is updated, the position of the zero level set is altered and the LSF is subsequentlyupdated at new locations. Consequently, it is beneficial to maintain an LSF that is51Chapter 3. ABLE: an activity-based level set segmentation algorithmFigure 3.3: The Chan-Vese active contours model. In A, we display the initialestimate of the object boundary, Γ(0). The external velocity at each pixel is equal toa data-based term, Ain(x) − Aout(x), multiplied by a localisation operator. In B, wedisplay the former, which is negative at pixels that are more similar to the interior thanthe exterior of the active contour, and positive otherwise. The localisation operatorsuppresses the contribution of the data-based term at pixels distant from the activecontour (C). In D, we plot the updated active contour, whose evolution is driven bythe external velocity in C.smoothly varying in that region. We use the regulariser of Li et al. [69], which has beenimplemented in numerous medical imaging applications [68, 72, 146]. The goal of theregulariser is to ensure that φ approximates a signed distance function in the vicinityof the contour and is flat elsewhere. In this context, a distance function satisfies:d(x) = minxc∈Γ‖x− xc‖, (3.15)i.e. the value of d(x) is the distance from x to the closest point on the contour. Itfollows that, for all x on the contour, the distance is equal to 0. A signed distancefunction encodes both distance to the contour and the subregion that a point belongsto; it is equal to d(x) for pixels in the curve interior and −d(x) for those in the exterior,see Fig. 3.2B. A signed distance function, therefore, satisfies the properties of an LSFand is sufficiently smoothly varying. The regularisation energy isEreg(φ) =∫Ωp (‖∇φ‖) dx, (3.16)wherep(s) =1(2π)2(1− cos(2πs)) if s ≤ 1,12 (s− 1)2 otherwise.(3.17)523.4. Proposed methodThe partial derivative is then [69]∂Ereg∂φ= div(∇φp′(‖∇φ‖)‖∇φ‖), (3.18)where div(·) is the divergence operator.3.4 Proposed methodIn the approaches presented thus far, an active contour partitions the whole imagedomain into two distinct subregions: the object and the background. In the case ofcalcium imaging data, the objects to be detected (cells or cellular components) aremultiple orders of magnitude smaller than the spatial domain. We therefore adopta more local approach to the segmentation problem. In the following, we start byconsidering the problem of detecting the boundary of a single, isolated cell.In our framework, an active contour seeks to partition a local region of the video intoa cell interior and exterior. At each iteration, the cell exterior is defined as the setof pixels exterior to the contour within a fixed distance, see Fig. 3.4B. The defaultdistance is taken to be two times the expected radius of a cell. We refer to this exteriorregion as the narrowband to emphasise its proximity to the contour of interest. As anactive contour is updated, so is the corresponding narrowband (Fig. 3.4F). The regionof the video for which the optimal partition is sought is therefore not static; rather, itevolves as an active contour evolves.Given an instance of the active contour, Γ(τ), the corresponding LSF is defined on thecontour interior and narrowband (see Fig. 3.4C). We compute f in ∈ Rm and fout ∈ Rmas the average cell interior and narrowband time courses, respectively, where m is thenumber of frames in the video. As before, the external energy isEext(φ) =∫φ>0D(I(x), f in)dx +∫φ≤0D(I(x), fout)dx. (3.19)where we take λout = λin = 1.In the following section, we introduce the dissimilarity metrics used. Then, in Section3.4.2, we extend the framework to multiple active contours.53Chapter 3. ABLE: an activity-based level set segmentation algorithmFigure 3.4: Flow diagram of the proposed segmentation algorithm. Thethree stages are: initialization (A-C), iterative updates of the estimate (D-G) andconvergence (H-J). When cells are sufficiently far apart we can segment them inde-pendently — in this example we focus on the isolated cell in the dashed box on themaximum intensity image in A. We make an initial estimate of the cell interior, fromwhich we form the corresponding narrowband (B) and level set function φ (C). Basedon the discrepancy between a pixel’s time course and the time courses of the interiorand narrowband regions (D), we calculate the velocity of φ at each pixel (E). φ evolvesaccording to this velocity (G), which updates the location of the interior and narrow-band (F). Final results for: one cell (H), the average signals from the correspondinginterior and narrowband (I) and segmentation of all four cells (J).543.4. Proposed method3.4.1 Dissimilarity metricsDue to the heterogeneity of calcium imaging data, we do not use a universal dissimi-larity metric. When both the pattern and magnitude of a pixel’s temporal activity areinformative, as is typically the case for synthetic dyes, we useDE(I(x), f)=1m‖I(x)− f‖2, (3.20)for f ∈ Rm. When we have an image not a video (i.e. I(x) and f are one-dimensional)this dissimilarity metric reduces to the fitting term introduced by Chan and Vese [24].For datasets in which the fluorescence expression level varies throughout cells and, asa consequence, pixels in the same cell exhibit the same pattern of activity at differentmagnitudes, we use a metric based on the correlation, such thatDC(I(x), f)= 1− corr(I(x), f), (3.21)where corr represents the Pearson correlation coefficient.3.4.2 Neighbouring cellsWe now extend the cost function presented in Eq. (3.19) to one suitable for parti-tioning a region into multiple cell interiors and a global exterior, which encompassesthe narrowbands of all the cells. We denote with φm the LSF that corresponds tothe mth cell and f in,m the corresponding interior time course. Due to the relativelylow axial resolution of a two-photon microscope, fluorescence intensity at one pixelcan originate from multiple cells in neighbouring z-planes. Therefore, in contrast toprevious approaches to coupling active contours [156, 40], we do not penalise overlapof contour interiors. Instead, we allow them to overlap when this best fits the data. Inparticular, we assume that a pixel in multiple cells would have a time course well fitby the sum of the interior time courses for each cell. The external energy in the caseof M > 0 cells is thusEext(φ1, φ2, ..., φM ) =∫outsideD(I(x), fout)dx +∫insideD(I(x),∑m∈C(x)f in,m)dx (3.22)where ‘outside’ denotes the global exterior, the union of the narrowbands excludingany cell interior, and ‘inside’ denotes the union of the contour interiors. The function55Chapter 3. ABLE: an activity-based level set segmentation algorithmC(x) identifies all cells whose interior contains x ∈ Ω. We denote the local spatialdomain, which incorporates the global exterior and the union of contour interiors,with Ω̄. When the region to be partitioned contains only one cell, the external energyin Eq. (3.22) reduces to that in Eq. (3.19).We can rewrite the RHS of Eq. (3.22) in terms of the LSFs using the Heaviside functionimplicitly as an indicator function, such thatEext(φ1, φ2, ..., φM ) =∫Ω̄A(x)∏k∈C(x)H(φk(x))∏n6∈C(x)(1−H(φn(x))) dx, (3.23)where we simplify the integrand by writingA(x) =D(I(x),∑m∈C(x)f in,m)if C(x) 6= ∅,D(I(x), fout)otherwise.(3.24)To obtain the external velocity of the mth active contour we compute the derivative ofEq. (3.23) with respect to φm. We state the result here and provide the derivation inAppendix A.1. As in Eq. (3.12), the velocity is the product of a data-based term anda localisation operator, such that∂Eext∂φm(x) = δε(φm(x))(D(I(x), f in,m)−D(I(x), fout)), (3.25)if x is not in a neighbouring cell, i.e. H(φk(x)) = 0 for all k 6= m. Otherwise∂Eext∂φm(x) = δε(φm(x))(D(I(x), f in,m +∑j∈C(x)j 6=mf in,j)−D(I(x),∑j∈C(x)j 6=mf in,j)). (3.26)The external velocity can be interpreted as follows: if a pixel, not in another cell, hastime course more similar to that of the contour interior than the narrowband, thenthe contour moves to incorporate that pixel. If a pixel in another cell has time coursebetter matched by the sum of the interior time courses of cells containing that pixelplus the interior time course of the evolving active contour, then the contour moves toincorporate it. Otherwise, the contour is repelled from that pixel.Although the external velocity contains contributions from all cells in the video, thelocalisation operator, δε, ensures that the problem remains local. We only evaluate theimage-based terms in Eq. (3.25) and (3.26) at pixels where the localisation operator is563.4. Proposed methodFigure 3.5: Evolution of multiple active contours. Each active contour seeksto partition its local region into a cell interior and exterior. When contours becomesufficiently close, their evolution is coupled.non-zero. As noted in Section 3.3.1, this is a band of pixels around the active contour,with width on the order of a single pixel. Consequently, many active contours can beevolved independently. Mathematically, fout is the average time course correspondingto the global exterior, which is the union of all contours’ narrowbands. In practice,when we compute the velocity of each active contour we utilise a local fout, which iscomputed from pixels in that contour’s narrowband.The active contours are evolved in parallel, so that, at each time step, τ , we computeφm(τ + ∆τ) = φm(τ)−∆τ(λ∂Eext∂φm+ µ∂Ereg∂φm). (3.27)for all m ∈ {1, ...,M}, where the regulariser is the same as in Section 3.3.2. Then, thefeatures, f in,m and fout, are updated and each contour is evolved again. This processis repeated until a maximum number of iterations has been reached or the contourshave converged. In the next section, we discuss details of the implementation, such asthe conditions for convergence and the initialisation procedure.3.4.3 ImplementationInitialisationWe devised an automatic initialization algorithm which selects connected areas of peaklocal correlation and/or peak mean intensity as initial ROI estimates. The correlation57Chapter 3. ABLE: an activity-based level set segmentation algorithmNum. cells25 125 225Num. frames (m)100 1.1 6.5 11.21000 1.3 6.5 12.7Table 3.1: Runtime of the segmentation algorithm (minutes) on syntheticdata. On data with dimensions 512 × 512 × m, the runtime of ABLE (minutes)increases linearly with the number of cells and is not significantly affected by increasingnumber of frames, m. Runtime was measured on a PC with 3.4GHz Intel Core i7 CPU.image displays, at each pixel, the average correlation between that pixel’s time courseand those of neighbouring pixels. Local peaks in this image and the mean intensityimage are identified (by a built-in MATLAB function, ‘imextendedmax’) as candidateROIs. The selectivity of the initialization is set by a tuning parameter, which definesthe relative height with respect to neighbouring pixels (in units of standard deviationof the input image) of the peaks that are suppressed. The higher the value of thisparameter, the more conservative the initialization. Empirically, we have found valuesin the range 0.2 - 0.8 to be suitable. This typically overestimates the number ofROIs; redundant estimates are automatically pruned during the update phase of thealgorithm.ConvergenceWe stop updating a contour estimate when a maximum number of iterations, Nmax,has been reached or the active contour has converged — using one or both of theseconditions is common in the active contour literature, see, for example, [29, 69]. Acontour is deemed to have converged if, in Ncon consecutive iterations, the number ofpixels that are added to or removed from the interior is less than ρ. As default, wetake Nmax = 100, Ncon = 40 and ρ = 2.The complexity of the level set method is intrinsically related to the dimensionality ofthe active contour; the number of frames of the video is only relevant to the evaluationof the external velocity, which accounts for a small fraction of the computational cost.In Table 3.1, we demonstrate that increasing video length by a factor of 10 has onlya minor impact on processing time. Updating an active contour is generally a localproblem — consequently, we observe that algorithm runtime increases linearly withthe total number of cells, see Table 3.1. Due to the independence of spatially separate583.5. Experimental methodsROIs in our framework, further performance speed-ups are achievable by parallelisingthe computation.Merging and pruning regionsABLE automatically merges two cells if they are sufficiently close and their interiorssufficiently correlated — a strategy previously employed in the CNMF algorithm ofPnevmatikakis et al. [105]. When two contours are merged, their respective level setfunctions are replaced with a single level set function, initialized as a signed distancefunction, with a zero level set that represents the union of the contour interiors. Therequired proximity for two cells to be merged is one cell radius (the expected cell radiusis one of two required user input parameters) and the default correlation threshold is0.8. A contour is automatically removed (‘pruned’) at the end of the update phase ifits area is smaller or greater than adjustable minimum or maximum size thresholds,which, as default, are set at 3 and 3πr2 pixels, respectively, where r is the expectedradius of a cell.3.5 Experimental methodsIn Section 3.6, ABLE’s performance is assessed on a range of real and simulated data.In the following, we describe the experimental methods used to obtain this data andthe performance metrics used.3.5.1 SimulationsCellular spike trains are generated from mutually independent Poisson processes. Acell’s temporal activity is the sum of a stationary baseline component, the value ofwhich is selected from a uniform distribution, and a spike train convolved with astereotypical calcium transient pulse shape. Cells are ‘doughnut’ (annulus) shaped tomimic videos generated by genetically encoded calcium indicators, which are excludedfrom the nucleus. To achieve this, the temporal activity of a pixel in a cell is generatedby multiplying the cellular temporal activity vector by a factor in [0, 1] that decreasesas pixels are further from the cell boundary. When two cells overlap in one pixel,59Chapter 3. ABLE: an activity-based level set segmentation algorithmwe sum the contributions of both cells at that pixel. Spatially and temporally vary-ing background activity, generated independently from the cellular spiking activity, ispresent in pixels that do not belong to a cell.3.5.2 Two-photon calcium imaging of hippocampal slicesThis dataset was collected by Dr Renaud Schuck, our collaborator in Imperial CollegeLondon’s Department of Bioengineering. Here, we provide an overview of the pro-cedures used. For further details, we refer the reader to [117]. All procedures wereperformed in accordance with national and institutional guidelines and were approvedby the UK Home Office under Project License 70/7355 to Simon R Schultz. Juvenilewild-type mice of either sex (C57Bl6, P13-P21) were anaesthetised using isofluraneprior to decapitation procedure. Brain slices (400 µm thick) were horizontally cut in1-4◦C ventilated (95% O2, 5% CO2) slicing Artificial Cerebro-Spinal Fluid (sACSF:0.5 mM CaCl2, 3.0 mM KCl, 26 mM NaHCO3, 1 mM NaH2PO4, 3.5 mM MgSO4,123 mM Sucrose, 10 mM D-Glucose). Hippocampal slices containing Dentate Gyrus,CA3 and CA1 were taken and resting in ventilated recovery ACSF (rACSF: 2 mMCaCl2, 123 mM NaCl, 3.0 mM KCl, 26 mM NaHCO3, 1mM NaH2PO4, 2mM MgSO4,10mM D-Glucose) for 30min at 37◦C. After this the slices were placed in an incuba-tion chamber containing 2.5 mL of ventilated rACSF and ‘painted’ with 10 µL of thefollowing solution: 50 µg of Cal-520 AM (AAT Bioquest), 2 µL of Pluronic-F127 20%in DMSO (Life Technologies) and 48 µL of DMSO (Sigma Aldrich) where they wereleft for 30 min at 37◦C in the dark. Slices were then washed in rACSF at room tem-perature for 30 min before imaging. Dentate Gyrus granular cells were identified usingoblique illumination prior to being imaged using a standard commercial galvanometricscanner based two-photon microscope (Scientifica Ltd) coupled to a mode-locked MaiTai HP Ti Sapphire (Spectra-Physics) laser system operating at 810 nm. Functionalcalcium images of granular cells were acquired with a 40X objective (Olympus) byraster scanning a 180× 180 µm2 square Field of View at 10 Hz. Electrical stimulationwas accomplished with a tungsten bipolar concentric microelectrode (WPI) where thetip of the electrode was placed into the molecular layer of the Dentate Gyrus (20 pulseswith a pulse-width of 400µs and a 60µA amplitude were delivered into the tissue witha pulse repetition rate of 10 Hz, repeated every 40 sec).603.5. Experimental methods3.5.3 In vivo imaging datasetsWe demonstrate algorithm performance on two imaging videos made freely availablein online datasets [99, 85]. In the former dataset, the vibrissal cortex was imaged atvarious depths, from layer 1 to deep layer 3, whilst the mouse performed a pole local-ization task [100, 49]. Imaging was performed using the genetically-encoded indicatorGCaMP6s. Fluorescence was collected through raster scans of a 600×600µm2 regionat 7Hz. The latter dataset was recorded from the vibrissal primary somatosensorycortex of an awake head-fixed mouse in a virtual reality setup [123]. Images were ob-tained using indicator GCaMP6s. Fluorescence was collected through raster scans ofa 600×600µm2 region at 8Hz.3.5.4 Performance metricsTo quantify segmentation accuracy we use two sets of metrics. On simulated data,we compute the success rate of algorithm output compared to the ground truth celllocations. For each cell in the ground truth dataset, we identify the correspondingcell in the algorithm output. We then compute the precision, which is the percentageof pixels in an estimated cell that are also in the corresponding ground truth cell.Similarly, we compute the recall, which is the percentage of pixels in the ground truthcell that are in the estimated cell. The success rate (%) is the harmonic mean of theprecision and recall, such thatsuccess rate = 2precision ∗ recallprecision + recall. (3.28)For each simulated imaging video, the output of the metric is a single score, which isthe average of the success rates over all cells in the video.On real data with ground truth stemming from manually labelled cells, we use adifferent approach. In this case, the ground truth cell locations are not likely tobe accurate on the order of single pixels. Therefore, we identify an estimated cellas matching a ground truth cell if their centres are within 5 pixels of one another.We compute the recall as the proportion of ground truth cells that are matched bya detected cell. We also compute the fall-out-rate, which is the complement of theprecision. This is the proportion of detected cells that do not match a ground truthcell.61Chapter 3. ABLE: an activity-based level set segmentation algorithmFigure 3.6: Segmentation results on in vivo imaging data. ABLE detectscells with varying size and shape from mouse in vivo imaging data. The detected cellboundaries are superimposed on the correlation image of the imaging video (A). Ex-tracted neuropil-corrected time series and corresponding cell boundaries are displayedfor a subset of the detected regions: both cell bodies (B, C and D) and cross-sectionsof dendrites (E, F and G) are detected .3.6 Results3.6.1 Robustness to heterogeneity in cell propertiesABLE detected region of interest (ROIs) with diverse properties from the publiclyavailable mouse in vivo imaging dataset of [99], see Fig. 3.6. A video demonstrationof segmentation on this data along with the code to generate the results is available atgithub.com/StephanieRey.To maintain a versatile framework we included no priors on cellular morphology in thecost function that drives the evolution of an active contour. This allowed ABLE todetect irregularly-shaped ROIs, such as cell bodies with dendrites attached (Fig. 3.6B,C and D). No prior information on stereotypical neuronal temporal activity is includedin our framework. As such, ABLE detected cross-sections of dendrites, which do notdisplay stereotypical calcium transient activity (Fig. 3.6E, F and G). Activity fromdendritic cross-sections is an essential indicator of neuronal activity in lower corticallayers, at depths where imaging is not possible [98]. We note that, due to ABLE’streatment of overlapping cells, it was possible to identify the overlapping dendritesbelonging to the cells in Fig. 3.6F and G.623.6. ResultsFigure 3.7: Segmentation results on simulated data containing overlappingcells. ABLE detected the true boundaries of overlapping cells from noisy simulateddata. The detected contours for one realization of noise with standard deviation (σ)equal to 60 are plotted on the correlation image in A. Given an initialisation on afixed grid, the true cell boundaries are detected with success rate of at least 99% forσ < 90 (B). The central marker and box edges in B indicate the median and the 25thand 75th percentiles, respectively. For noise level reference, we plot the average timecourse from inside the green contour in A at various levels (C). In D, we plot theinitialisation, which was on a fixed grid, on the mean image of the video.3.6.2 Demixing overlapping cellsWhen imaging through scattering tissue, a two-photon microscope can have relativelylow axial resolution (on the order of ten microns) in comparison to its excellent lateralresolution. As a consequence, the photons collected at one pixel can in some casesoriginate from multiple cells in a range of z-planes. For this reason, cells can appear tooverlap in an imaging video. It is crucial that segmentation algorithms can delineatethe true boundary of apparently overlapping cells, so that the functional activity ofeach cell can be correctly extracted and analysed.On synthetic data containing 25 cells, 17 of which had some overlap with anothercell, we measured the success rate of ABLE’s segmentation compared to the groundtruth cell locations (Fig. 3.7). For more detail on the performance metric, see Section3.5.4. Cells were simulated with uneven brightness to mimic the ‘doughnut’ cellsgenerated by some genetically encoded indicators that are excluded from the nucleus.Consequently, pixels within the same cell had the same pattern of temporal activity63Chapter 3. ABLE: an activity-based level set segmentation algorithmFigure 3.8: Segmentation results on in vitro imaging data containing syn-chronously spiking cells. The boundaries of the detected ROIs are superimposedon the thresholded maximum intensity image (A). Panel B displays the time coursesof the detected ROIs (each plotted as a row of the matrix) along with the video meanraw activity and the electrical stimulations. In B, the time courses are median filteredfor visual clarity. ABLE detected both active (C, D and E) and inactive ROIs (F, Gand H). We display the contours of the two detected ROIs on the correlation image(C and F), the mean image (D and G) and the cellular time courses (E and H).but vastly different baseline fluorescence, as is apparent in the mean image in Fig.3.7D. To segment the data, therefore, we applied ABLE with the correlation-baseddissimilarity metric. As can be seen in the example in Fig. 3.7A, despite unevenfluorescence and the prevalence of cellular overlap, pixels were reliably incorporatedinto the correct cell. We tested segmentation performance as the noise level of thedata varied, measuring performance over 10 realisations of noise at each noise level. Toeliminate variability due to the quality of the initialisation, we initialised cell estimateson a fixed grid (see Fig. 3.7D). On average, ABLE achieved success rate greater than99% when the noise standard deviation was less than 90 (Fig. 3.7B). Segmentationperformance deteriorated as the noise standard deviation increased; ABLE achievedan average success rate of 79% when the noise standard deviation was equal to 200.643.6. Results3.6.3 Results on hippocampal brain slicesABLE detected synchronously spiking, densely-packed cells from mouse in vitro imag-ing data (Fig. 3.8). As the brain slice was electrically stimulated (at rate 10Hz for2s every 40s) during imaging, cellular activity in this dataset is highly correlated withthat of other cells and the background. To segment the dataset, therefore, we appliedABLE with the Euclidean dissimilarity metric. When the cell interior and narrow-band time courses are highly correlated, the external velocity of an active contour inthis version of the algorithm is predominantly driven by the discrepancy in baselineintensities. For the mathematical reasoning behind this statement, see Appendix A.2.As a consequence, the algorithm was able to segment cells whose activity was highlycorrelated with the background activity (Fig. 3.8B). Furthermore, inactive cells weredetected when their baseline fluorescence allowed them to be identified from the back-ground (Fig. 3.8 F-H). In Fig. 3.8C-E, we display an example of an active cell thatwas detected. We note that each of the four visible activity peaks in this figure corre-sponds to twenty calcium transients that were induced by the electrical stimulations.The relatively low amplitude of these peaks reflects the low SNR of this data. As aresult, only a small amount of spontaneous activity is visible in the activity plot ofFig. 3.8B.3.6.4 Algorithm comparison on manually labelled datasetWe compared the performance of ABLE with two state of the art calcium imaging seg-mentation algorithms, CNMF [105] and Suite2p [94], on a manually labelled dataset,see Fig. 3.9. We apply ABLE with the correlation-based dissimilarity metric, Eq.(3.21), which is well suited to neurons with low baseline fluorescence and unevenbrightness. As the dataset is large enough (512x512x8000 pixels) to present mem-ory issues on a standard laptop, we ran the patch-based implementation of CNMF,which processes spatially-overlapping patches of the dataset in parallel. We optimisedthe performance of each algorithm by selecting a range of values for each of a set oftuning parameters and generating segmentation results for all combinations of the pa-rameter set. The results were visualised on the correlation image and the parameterset that presented the best match to the correlation image was selected. This processis representative of what a user may do in practise when applying an algorithm to anew dataset.65Chapter 3. ABLE: an activity-based level set segmentation algorithmFigure 3.9: Algorithm comparison on a manually labelled dataset. Wecompare the segmentation results of ABLE, CNMF and Suite2p on a manually labelleddataset. On the correlation image we plot the boundaries of the manually labelled cellscolour-coded by the combination of algorithms that detected them (A), undetectedcells are indicated by a white contour. Suite2p detected the highest proportion ofmanually labelled cells (B), whereas ABLE had the lowest fall-out rate (C), which isthe percentage of detected regions not present in the manual labels. Some algorithm-detected ROIs that were not present in the manual labels are detected by multiplealgorithms (D) and have time courses which exhibit stereotypical calcium transientactivity (E). The correlation image in D is thresholded to enhance visibility of localpeaks in correlation. In E, we plot the extracted time courses of the ROIs in D.ABLE achieved the highest success rate (67.5%) when compared to the manual labels,see Table 3.2. For a definition of the success rate and other performance metricsused, see Section 3.5.4. ABLE achieved a lower fall-out rate than Suite2p and CNMF(Fig. 3.9C) — 67.5% of the ROIs it detected matched with the manually labelledcells. Some of the ‘false detections’ were consistent among algorithms (Fig. 3.9C)and corresponded to local peaks in the correlation image (Fig. 3.9D), whose extractedtime courses displayed stereotypical calcium transient activity (Fig. 3.9E). A subset ofthese ROIs may thus correspond to cells omitted by the manual operator. The highestproportion of the manually labelled cells were detected by Suite2p, which detected thegreatest number of cells not detected by any other algorithm (Fig. 3.9B). A smallproportion (13.2%) of cells were detected by none of the algorithms. As can be seenfrom Fig. 3.9A, these do not correspond to peaks in the correlation image, and mayreflect inactive cells detected by the manual operator.663.7. SummarySuccess rate (%) Precision (%) Recall (%)ABLE 67.5 67.5 67.5CNMF 63.4 60.7 66.5Suite2p 63.7 56.5 73.1Table 3.2: Algorithm success rate on a manually labelled dataset. Wecompare the performance of three segmentation algorithms: ABLE, CNMF [105] andSuite2p [94], using the manual labels as ground truth.3.7 SummaryIn this chapter, we presented a novel approach to the problem of detecting cells fromcalcium imaging data. Our approach uses multiple coupled active contours to identifycell boundaries. The core assumption is that the local region around a single cell canbe well-approximated by two subregions, the cell interior and exterior. The averagetime course of the respective subregions is used as a feature with which to classifypixels into either subregion. We assume that pixels in which multiple cells overlaphave time courses that are well-approximated by the sum of each cell’s time course.We form a cost function based on these assumptions that is minimised when the activecontours are located at the true cell boundaries.ABLE is a flexible method: we include no priors on a region’s morphology or stereo-typical temporal activity. Due to this versatility, ABLE segmented cells with varyingsize and shape from a mouse in vivo dataset (Fig. 3.6) and both active and inactivecells from brain slices (Fig. 3.8). Furthermore, we demonstrated its ability to detectthe true boundaries of overlapping cells on both real and simulated data.67Chapter 3. ABLE: an activity-based level set segmentation algorithm68Chapter 4.CosMIC: A consistent metric for spikeinference from calcium imagingTwo-photon calcium imaging does not directly read out a neuron’s spiking activity.Rather, it reads out a secondary variable — the intracellular calcium concentration— from which spiking activity must be inferred. After an experiment, neural codinghypotheses are investigated by comparing the rate and timing of inferred spike trainsto properties of stimuli and behavioural variables. To justify such analysis, the abilityof spike inference algorithms to generate accurate spike train estimates must first beverified.Inferred spikes are compared to the electrophysiological ground truth (see Section 2.2)in order to assess their accuracy. Whilst there are a rich collection of ground truthdatasets [7], the metrics currently used to compare ground truth and estimated spiketrains are not ideal. In this chapter, we highlight the limitations of the two mostcommonly used spike inference metrics, the spike train correlation and success rate,and propose an alternative, which we refer to as CosMIC. In the following, we introduceCosMIC and demonstrate its application to real and simulated spike trains.4.1 IntroductionAlthough the development of spike inference algorithms has received a lot of recentattention, few researchers have examined the metrics used to assess an algorithm’sThe work presented in this chapter led to the following publication [108].69Chapter 4. CosMIC: A consistent metric for spike inferenceFigure 4.1: Inconsistency in usage of spike inference metrics. This chartillustrates the number of papers, of the 16 presenting a spike inference algorithm, thatuse each metric to assess their algorithm. Of the surveyed papers, 44% use a metricunique to that paper and 57% use more than one metric. Under the umbrella of thesuccess rate, we include the true and false positive rates, from which the success rateis computed, and the error rate — the complement of the success rate.performance. At present, there is no consensus regarding the best choice of metric.In fact, 44% of papers presenting a new method assess its performance using a metricunique to that paper (Fig. 4.1). This inconsistency impedes progress in the field— algorithms are not directly comparable and, consequently, data collectors cannoteasily select the optimal algorithm for a new dataset. In light of this problem, a recentbenchmarking study tested a range of algorithms on a wide array of imaging data [7].Although informative, the study, which relied heavily on the spike train correlation(STC) to assess algorithm performance, may not provide the full picture.Despite being the most commonly used metric, the STC is not well-suited to the task.It is invariant under linear transformations of the inputs and is therefore unable toevaluate the similarity of the rates of two spike trains [95]. The second most commonlyused metric, the success rate, does not reward increasing temporal precision abovea given threshold. Consequently, it is not an appropriate metric for evaluating analgorithm’s performance when the end goal is, for example, to investigate the synchronyof activations within a network.Due to limiting factors, such as noise and model mismatch, it is improbable that an es-timate will match a true spike with infinite temporal precision. In order to incorporatean implicit temporal tolerance, metrics manipulate the sets of spikes before compar-ing their similarity. One such technique is to discretise the temporal interval into afinite number of ‘bins’ and to compare the corresponding vectors of spike counts. Thismethod, which is employed by the STC, has its limitations. In particular, estimates704.1. Introductionthe same absolute distance from a true spike can fall into different bins and thus receivedifferent scores. Furthermore, estimates in the same bin are treated equally and, as aconsequence, temporal precision above the bin width is not rewarded. An alternativestrategy is to convolve spike trains with a smoothing pulse and to assess the similarityof the resulting pulse trains. This type of continuous approach is preferred by metricsassessing the relationship between spike trains from different neurons [144, 116].In the following, we propose a metric that we refer to as CosMIC (a Consistent Metricfor spike Inference from Calcium imaging). CosMIC first convolves the spike trainswith a smoothing pulse whose width reflects the temporal precision that an estimateis able to achieve given the limitations of the data. To set this width, we calculate alower bound on the temporal precision of the estimate of one spike, the Cramér-Raobound (CRB). In the context of calcium imaging, this bound has been previously usedto evaluate detectability of spikes under different imaging modalities [117]. We usethe CRB as a performance benchmark; an estimator that achieves the CRB should beawarded a relatively high score.We formalise the comparison of the resulting pulse trains by leveraging results fromfuzzy set theory. In contrast to classical sets, to which an element either belongs ordoes not belong, fuzzy sets contain elements with a level of certainty represented bya membership function — the higher the value of the membership function, the morecertain the membership. We view the pulse trains as the membership functions ofthe fuzzy sets of true and estimated spikes. These sets represent the original sets ofspikes with a level of temporal tolerance set by the pulse width. CosMIC computesthe similarity of the original sets of spikes from the size of the intersection of the fuzzysets as a fraction of their average size. In this context, CosMIC resembles a continuousSørensen-Dice coefficient — an index commonly used to assess the similarity of discrete,presence/absence data [32, 125].In Section 4.2, we give a brief introduction to fuzzy set theory. Then, in Section 4.3, wedetail the computation of CosMIC. The relationship between the pulse width and theCRB is explained in Section 4.4. In Section 4.6, we demonstrate CosMIC’s ability todiscriminate spike train similarity on real and simulated data. CosMIC’s discriminativeperformance is compared to that of the success rate and the STC, whose definitionswe provide for completeness in Section 4.5.71Chapter 4. CosMIC: A consistent metric for spike inference4.2 Fuzzy set theoryFuzzy set theory was introduced by Zadeh in 1965 to extend classical set-theoreticconcepts to areas in which abrupt set boundaries are hard to define [155]. As an illus-trative example, Zadeh considered a universe consisting of all men and a set containingall tall men. As with many semantic definitions, it is hard to identify where the setof tall men ends and the set of ‘not tall’ men begins. Instead of forcing elements toeither belong or not belong to a set, as in classical set theory, it was proposed thatsets should have graded membership. Accordingly, a non-negative value is associatedwith each element of the universe, representing the degree of certainty to which thatelement belongs to a set. Here, we provide a brief introduction to fuzzy sets. Werefer the interested reader to [39] for an extensive theoretical overview and to [157] formore recent developments. To differentiate between classical sets and fuzzy sets, weuse plain font for the former (S) and script-like font for the latter (S).In classical set theory, membership of a set, S, is defined for all elements of a universe,X, using a characteristic function:1S(x) =1 if and only if x ∈ S,0 if and only if x 6∈ S. (4.1)In fuzzy set theory, the range of the analogous function, µS(x), which we refer to asa membership function, is the non-negative, real numbers. The value of µS increasesas the certainty of membership does. A fuzzy set is characterised by a set of orderedpairs:S = {(x, µS(x)) : x ∈ X}, (4.2)each representing an element of the universe and the certainty with which that elementbelongs to the set. It is not the absolute value of a membership function that is ofimportance. Rather, it is the implicit ordering induced by a membership function thatwe are interested in. For example, we may want to compare the values of µS(x) andµS(y) or, alternatively, µS1(x) and µS2(x). We note that the characteristic functionof a classical set is a special case of a membership function that only takes the values{0, 1}.The operations of intersection and union were introduced to fuzzy sets by Zadeh [155]and, at a later date, justified axiomatically by Bellman and Giertz [6]. The membership724.2. Fuzzy set theoryBAFigure 4.2: Membership functions of the union and intersection of fuzzysets. In A, the membership functions are continuous-valued. In B, they are binaryand are therefore analogous to the characteristic functions of classical sets.function of an intersection of fuzzy sets is computed asµS1∩S2(x) = min (µS1(x), µS2(x)) , ∀x ∈ X. (4.3)Similarly, the membership function of the union isµS1∪S2(x) = max (µS1(x), µS2(x)) , ∀x ∈ X. (4.4)These definitions coincide with the classical operations of intersection and union in thespecial case in which the membership functions are characteristic functions. We showan illustrative example of the above operations in Fig. 4.2.When the universe is a finite set, the cardinality of a fuzzy set, |S|c, is [39]|S|c =∑x∈XµS(x). (4.5)When the universe is infinite, the cardinality of a set may not exist. In the case thatX is a measurable set and P is a measure on X [113], then|S|c =∫XµS(x) dP(x). (4.6)In particular, when the universe is the real line, we have|S|c =∫RµS(x) dx. (4.7)73Chapter 4. CosMIC: A consistent metric for spike inferenceFinally, the support of S is: supp(S) = {x ∈ X : µS(x) > 0}.4.2.1 Sørensen-Dice coefficientThe Sørensen-Dice coefficient (SDC) was designed to quantify similarity between sam-ples of discrete, presence/absence data [125, 32]. This metric, which is also known asthe F1-score, is widely used in many fields, including ecology [102] and image segmen-tation [159]. Although its initial justification was empirical rather than theoretical, itcan be interpreted as a similarity metric for fuzzy sets with binary membership func-tions belonging to a finite universe [96, 137]. This connection was made after the SDChad become established in numerous fields. For example, an ecologist, Roberts, usedfuzzy sets as a theoretical basis for a popular technique for environmental analyses[111]. In the following, we detail the computation of the SDC from the perspective offuzzy set theory.Let S1 and S2 be two fuzzy sets defined on a finite universe, X = {x1, x2, ..., xm},represented by the membership functions µS1 and µS2 , respectively. We assume themembership functions take values in the range {0, 1}. Then, the SDC isSDC(S1,S2) =|S1 ∩ S2|c(|S1|c + |S2|c) /2(4.8)=∑mi=1 min(µS1(xi), µS2(xi))(∑mi=1 µS1(xi) + µS2(xi)) /2, (4.9)where Eq. (4.9) follows from Eq. (4.3) and Eq. (4.5).When S1 and S2 represent a set of ground truth items and a set of estimates, respec-tively, the coefficient has two alternative formulations in terms of the types of errors.It is perhaps the intuitive appeal of these expressions that has led to the SDC’s pop-ularity. To illustrate these alternative formulations, we consider the example of imagesegmentation, where we want to quantify the proportion of pixels shared by a groundtruth mask and an estimate. The binary formulation of the SDC is well-suited to thissituation; if the ground truth mask contains the ith pixel, then µS1(xi) = 1, otherwiseµS1(xi) = 0. Similarly, if the ith pixel is in the estimated mask then µS2(xi) = 1,otherwise µS2(xi) = 0. We write the sets of true positives (TP), false positives (FP)744.2. Fuzzy set theoryand false negatives (FN) as:TP = {xi : µS1(xi) = 1 ∩ µS2(xi) = 1}, (4.10)FP = {xi : µS1(xi) = 0 ∩ µS2(xi) = 1}, (4.11)FN = {xi : µS1(xi) = 1 ∩ µS2(xi) = 0}. (4.12)The number of pixels in the ground truth mask can be calculated from the number oftrue positives and false negatives, such thatm∑i=1µS1(xi) = |{xi : µS1(xi) = 1}|c (4.13)= |{xi : µS1(xi) = 1 ∩ µS2(xi) = 1} ∪ {xi : µS1(xi) = 1 ∩ µS2(xi) = 0}|c (4.14)= |TP|c + |FN|c. (4.15)Eq. (4.15) follows from the fact that the cardinality of the union of disjoint sets is thesum of the cardinalities. By analogous logic, it follows thatm∑i=1µS2(xi) = |TP|c + |FP|c, (4.16)andm∑i=1min(µS1(xi), µS2(xi)) = |TP|c. (4.17)Then, substituting these identities into Eq. (4.9), the SDC can be expressed in equiv-alent form:SDC(S1,S2) =2 |TP|c2 |TP|c + |FP|c + |FN|c. (4.18)Alternatively, defining the precision as the proportion of pixels (or, more generally,items) in the estimated mask that are true positives and the recall as the proportionof ground truth pixels that are true positives, the SDC is, equivalently, the harmonicmean of the precision and recall:precision ∗ recall(precision + recall) /2=|TP|c|TP|c+|FP|c ∗|TP|c|TP|c+|FN|c(|TP|c|TP|c+|FP|c +|TP|c|TP|c+|FN|c)/2(4.19)=2|TP|c2|TP|c + |FP|c + |FN|c(4.20)= SDC(S1,S2). (4.21)75Chapter 4. CosMIC: A consistent metric for spike inferenceA generalised SDC that quantifies common abundance is used to assess similarity ofnon-binary data [13]. The only difference in the calculation of this coefficient fromEq. (4.8) is that the entries of µS1 and µS2 can assume any non-negative value. Itis common that the ith entries of µS1 and µS2 represent counts of the correspondingobjects at the ith location. Interpreting these frequency counts as probabilities ofobject occurrence is consistent with an empirical approach to probability. In the fuzzyset literature, this coefficient has been further generalised to compare non-quantitativedata [137].4.3 Constructing the metricA spike inference metric compares the similarity between two sets of spikes: a groundtruth set, S = {tk}Kk=1, and a set of estimates, Ŝ = {t̂k}K̂k=1. In the following, wedefine two fuzzy sets, Sε and Ŝε, which represent the original sets of spikes, S and Ŝ,with a level of temporal tolerance defined by a parameter ε. We set ε to reflect thetemporal precision that an estimate is able to achieve given the statistics of a dataset(see Section 4.4.4). The corresponding membership functions are defined on the realline and are denoted µ(t) and µ̂(t). The membership functions are obtained throughconvolution of the spike trains,x(t) =K∑k=1δ(t− tk) and x̂(t) =K̂∑k=1δ(t− t̂k), (4.22)with a triangular pulse. The resulting functions have local maxima at the locationsof the respective sets of spikes (Fig. 4.3A). As x(t) and x̂(t) are analogous to themembership functions of the classical sets of spikes, we can think of the convolutionas a temporal smoothing of the membership. The pulse that we employ is a triangularB-spline (Fig. 4.3B),qε(t) =ε−|t|ε |t| ≤ ε,0 otherwise.(4.23)Using this triangular pulse means that, the further a time point is from a spike, theless weight the membership function receives at that point. Past a certain distance,ε, the membership function receives no weight. Many pulse shapes could be chosento introduce this grading of temporal precision, we select a triangular pulse as it is764.3. Constructing the metricFigure 4.3: Flow diagram of the proposed metric. The ground truth spiketrain and estimated spike train are convolved with a triangular pulse (B), whose widthis determined by the statistics of the data. The metric compares the difference betweenthe resulting pulse trains (A). Metric scores are in the range [0,1] — a perfect estimateachieves score 1 and an empty spike train is scored 0 (C).straightforward to examine analytically and implement computationally.We design the proposed metric to quantify the size of the intersection of the fuzzy setsof true and estimated spikes with respect to the average size of the sets, such thatM(S, Ŝ)=|Sε ∩ Ŝε|c(|Sε|c + |Ŝε|c)/2. (4.24)From Eq. (4.8), the analogy between this metric and the SDC is clear. Whereasthe SDC operates on binary membership function defined on finite universes, CosMICoperates on smooth membership functions defined on the real line. From the definitionsof intersection, Eq. (4.3), and cardinality, Eq. (4.7), we haveM(S, Ŝ)= 2∫Rmin(µ(t), µ̂(t)) dt∫Rµ(t) dt+∫Rµ̂(t) dt= 2‖min(µ, µ̂)‖1‖µ‖1 + ‖µ̂‖1, (4.25)77Chapter 4. CosMIC: A consistent metric for spike inferenceFigure 4.4: Proposed metric example. The proposed metric quantifies thecommonalities of the sets of true and estimated spikes as a proportion of the averagesize of those sets. Commonalities are found by taking the minimum of the pulse trains— as such, spikes that appear in only one pulse train are excluded (A) and estimateswith lower temporal precision receive a lower score (B and C).where we denote the continuous L1-norm with ‖f‖1 :=∫R |f(t)| dt. We note thatthe absolute value is not necessary in the integrand of Eq. (4.25) as the membershipfunctions are non-negative. Taking the minimum of the membership functions producesa conservative representation of the intersection of two sets; in our context, spikes thatappear in one spike train and not in the other are removed (Fig. 4.4A) and spikesthat are detected with poor temporal precision are assigned less weight (Fig. 4.4B and4.4C).Defining D+ = {t ∈ R : µ(t) > µ̂(t)} and D− = {t ∈ R : µ(t) ≤ µ̂(t)}, we have‖µ− µ̂‖1 =∫R|µ(t)− µ̂(t)| dt =∫D+µ(t)− µ̂(t) dt+∫D−µ̂(t)− µ(t) dt, (4.26)and‖µ‖1 + ‖µ̂‖1 − ‖µ− µ̂‖1 = 2(∫D+µ̂(t) dt+∫D−µ(t) dt)(4.27)= 2∫Rmin(µ(t), µ̂(t)) dt (4.28)= 2‖min(µ, µ̂)‖1. (4.29)784.3. Constructing the metricIt follows that an alternative expression for the metric isM(S, Ŝ) =2‖min (µ, µ̂)‖1‖µ‖1 + ‖µ̂‖1(4.30)=‖µ‖1 + ‖µ̂‖1 − ‖µ− µ̂‖1‖µ‖1 + ‖µ̂‖1(4.31)= 1− ‖µ− µ̂‖1‖µ‖1 + ‖µ̂‖1. (4.32)From Eq. (4.32), it is clear that the maximal score of 1 is achieved when the mem-bership functions, and therefore the sets of true and estimated spikes, are equivalent.The minimal score of 0 is achieved when the support of the membership functions donot overlap, i.e. no estimates are within the tolerance of the metric.4.3.1 Ancestor metricsLike the SDC, CosMIC can alternatively be derived from a pair of metrics, which werefer to as ancestor metrics. The first of these metrics measures the proportion ofground truth spikes that were detected within the precision of the pulse width, suchthatRCosMIC =|Sε ∩ Ŝε|c|Sε|c=‖min(µ, µ̂)‖1‖µ‖1. (4.33)This score is analogous to the recall of a spike train estimate, one of the ancestormetrics from which the SDC is formed. The second of CosMIC’s ancestor metricsmeasures the proportion of estimated spikes that detect a ground truth spike withinthe precision of the pulse width, such thatPCosMIC =|Sε ∩ Ŝε|c|Ŝε|c=‖min(µ, µ̂)‖1‖µ̂‖1. (4.34)This is analogous to the precision, the second metric used to compute the SDC. Finally,computing the harmonic mean of the two ancestor metrics and rearranging, we obtainCosMIC:2RCosMIC ∗ PCosMICRCosMIC + PCosMIC= 2|Sε ∩ Ŝε|c|Sε|c|Sε ∩ Ŝε|c|Ŝε|c|Sε ∩ Ŝε|c|Sε|c+|Sε ∩ Ŝε|c|Ŝε|c= 2|Sε ∩ Ŝε|c|Sε|c + |Ŝε|c= M(S, Ŝ). (4.35)79Chapter 4. CosMIC: A consistent metric for spike inference4.4 Temporal error toleranceThe width of the triangular pulse with which the spike trains are convolved reflectsthe accepted tolerance of an estimated spike’s position with respect to the groundtruth. To set this width, we calculate a lower bound on the temporal precision of theestimate of one spike — the Cramér-Rao bound (CRB) — from the statistics of thedata. The CRB is a lower bound on the mean square error of any unbiased estimator.It is therefore useful as a benchmark; an estimator that achieves the CRB should beawarded a relatively high metric score. In the following, we discuss the mean squareerror of a parameter estimator (Section 4.4.1), introduce the CRB (Section 4.4.2) andcalculate this bound in the context of calcium imaging (Section 4.4.3). We relate thisbound to CosMIC’s pulse width in Section 4.4.4.4.4.1 Parameter estimationLet g(t, θ) be a signal evolving in time, t, which is dependent on a univariate, determin-istic parameter, θ ∈ R. We consider the situation where the parameter θ is unknownand the goal is to infer its value from noisy samples of g(t, θ). We haveg̃[n] = g[n, θ] + ξ[n] for n ∈ {0, ..., N − 1}, (4.36)where g[n, θ] = g(nT, θ) are deterministic samples of g(·, θ) at time resolution T andξ[n] are samples of the noise.Denoting the vector of data samples with g̃ = (g̃[0], g̃[1], ..., g̃[N − 1]) , we now considerthe error of the estimator, θ̂(g̃), constructed to infer the value of θ. The averageaccuracy of the estimator can be represented by the mean square error (MSE):MSE(θ̂) = E[(θ̂ − θ)2]. (4.37)The MSE can be decomposed into two distinct types of error, in what is known as thebias-variance tradeoff [62], such thatMSE(θ̂) =(E[θ̂]− θ)2+ var(θ̂). (4.38)The bias, E[θ̂] − θ, is a systematic error in the estimator’s output compared to thetrue value. The variance, var(θ̂), represents the degree of variability in the estimator’s804.4. Temporal error toleranceoutput. We restrict our attention to unbiased estimators, which are those that, onaverage, attain the correct value for all values of θ:E[θ̂]= θ. (4.39)Within this class of estimators, the one that achieves the minimum MSE is the onethat has the smallest variance. One approach to identifying an optimal estimator is,therefore, to identify within the class of unbiased estimators the one with the smallestvariance. To do this, one can compute a lower bound on the variance of an unbiasedestimator. This bound can subsequently be used as a benchmark; the quality of anestimator is judged by how close it is to attaining the bound.4.4.2 Cramér-Rao boundEach sample, g̃[n], is the sum of a deterministic signal component whose value isdependent on the parameter θ and a stochastic noise component, see Eq. (4.36). Thevector of samples, g̃, can thus be modelled by a PDF parametrised by θ: p(g̃; θ). Fromthese components, a lower bound on the variability of an unbiased estimator, θ̂, can becalculated. There are many lower bounds on the variance of unbiased estimators, forexample [158, 77]. Of all the bounds, the CRB is the most straightforward to computeand is accompanied by a set of conditions under which it is satisfied [62]. In Theorem4.1, we present, without proof, the CRB for a scalar, deterministic parameter, θ ∈ R.For a more detailed exposition of the theory, including the result for vector-valuedparameters and a proof, we refer the reader to [62].Theorem 4.1. The variance of an unbiased estimator satisfiesvar(θ̂)≥ 1−E[∂2 ln p(g̃; θ)∂θ2] , (4.40)if the following regularity condition holds:E[∂ ln p(g̃; θ)∂θ]= 0, ∀θ. (4.41)Here, the expectation is taken with respect to g̃ and the derivatives are evaluated at thetrue value of θ.We now derive the specific form of the CRB when the samples are independent and81Chapter 4. CosMIC: A consistent metric for spike inferenceidentically distributed (i.i.d.) and are corrupted by white Gaussian noise. In the caseof independence, we havep(g̃; θ) =N−1∏i=0pi(g̃[i]; θ), (4.42)where pi is the PDF of g̃[i]. When the samples are identically distributed, we havepi = q for all i ∈ {0, 1, ..., N − 1} and some PDF q. In particular, when ξ[n] are i.i.d.samples of white Gaussian noise with variance σ2, we haveln p(g̃; θ) = ln(N−1∏i=01√2πσ2exp(−(g̃[i]− g[i; θ])22σ2))(4.43)= −N2ln(2πσ2)− 12σ2N−1∑i=0(g̃[i]− g[i; θ])2 . (4.44)It follows that∂ ln p(g̃; θ)∂θ=1σ2N−1∑i=0(g̃[i]− g[i; θ]) ∂g [i ; θ]∂θ. (4.45)As E [g̃[i]− g[i; θ]] = 0, it is clear from Eq. (4.41) and the linearity of the expectationthat the regularity condition is satisfied when the derivative, ∂g[i ;θ]∂θ , exists. Differenti-ating again, we have∂2 ln p(g̃; θ)∂θ2=1σ2N−1∑i=0(g̃[i]− g[i; θ]) ∂2g[i; θ]∂θ2−(∂g [i ; θ]∂θ)2. (4.46)Taking the expectation removes the first term of the above formula, which has expec-tation equal to 0. As a result, the variance of an unbiased estimator in this settingsatisfies:var(θ̂)≥ σ2N−1∑i=0(∂g[i ;θ]∂θ)2 , (4.47)as long as the regularity condition in Eq. (4.41) is satisfied.4.4.3 Bound for spike inferenceWe consider the problem of estimating the location of one spike, t1, from noisy calciumimaging data. We assume that we have access to noisy samples, f̃ [n] = f [n] + ξ[n], forn ∈ {0, 1, ..., N − 1}. Here, f [n] = f(nT ) are samples of the fluorescence signal withtime resolution T and ξ[n] are i.i.d samples of white Gaussian noise with variance σ2.824.4. Temporal error toleranceThe fluorescence signal is modelled asf(t) =Ae−(t−t1)/τoff(1− e−(t−t1)/τon)for t > t10 otherwise,(4.48)where τon and τoff define the speed of the rise and decay of the fluorescence pulse,respectively, and A is the amplitude.From Eq. (4.47), the CRB of the parameter estimator, t̂1, isvar(t̂1) ≥σ2N−1∑i=0(∂f [n]∂t1)2 . (4.49)To compute the CRB, we must therefore evaluate the partial derivative of f(t) withrespect to t1 at t ∈ {0, T, 2T, ..., (N − 1)T}. For t > t1, we have∂f∂t1= Ae−(t−t1)/τoff(1τoff−(1τoff+ 1τon)e−(t−t1)/τon). (4.50)For t < 0, the derivative is equal to 0. To evaluate the derivative at t = t1, we considerthe left and right derivatives [134], starting with the former:limt→t−1f(t1)− f(t)t1 − t= limt→t−10− 0t1 − t= 0. (4.51)The right derivative islimt→t+1f(t)− f(t1)t− t1= limh→0+f(t1 + h)− f(t1)h(4.52)= limh→0+Ae−h/τoff(1− e−h/τon)h(4.53)= A limh→0+(1− hτoff +O(h2))(hτon+O(h2))h(4.54)= A limh→0+hτon+O(h2)h(4.55)=Aτon. (4.56)As the left and right derivatives of f(t) are not equal at t = t1, the derivative does notexist at that point. When we compute the CRB, therefore, we ensure that no samplesare taken at nT = t1.83Chapter 4. CosMIC: A consistent metric for spike inferenceFigure 4.5: Metric pulse width. The pulse width is set to reflect the temporalprecision achievable given the statistics of the dataset. We calculate the Cramér-Raobound (CRB), σ2CRB, a lower bound on the mean square error of the estimated locationof one spike from calcium imaging data (A). This bound decreases as the scan rate(Hz) and peak signal-to-noise ratio (squared calcium transient peak amplitude/noisevariance) increase. We set the pulse width to ensure that an estimate of one spike atthe temporal precision of the CRB achieves, on average, a score of 0.8. This results ina pulse width of approximately 7.3 σCRB (B).Substituting the value of the derivative into Eq. (4.49), we havevar(t̂1) ≥A2σ2N−1∑n=dt1/T ee−2(t−t1)/τoff(1τoff−(1τoff+ 1τon)e−(t−t1)/τon)2−1 . (4.57)In order that the CRB holds for an arbitrarily placed spike, we remove the dependencyon the true spike time by averaging the result over several values of t1. DenotingCRB(t1) as the right hand side of Eq. (4.57), we compute σ2CRB =1M∑Mm=1 CRB(tm1 ),where tm1 are evenly placed in the interval (nT, (n+ 1)T ) for a fixed n.In Fig. 4.5, we plot σCRB as the sampling rate and peak signal-to-noise ratio (PSNR) ofthe data vary. The PSNR is computed as A2peak/σ2, where Apeak is the peak amplitudeof the fluorescence signal (the maximum) and σ is the standard deviation of the noise.For this example, we use τon = 32ms and τoff = 314ms; the parameters for a Cal-520pulse [133]. We see that the CRB decreases as either the scan rate or the PSNR of thedata increases.844.5. Experimental methods4.4.4 Pulse widthThe CRB can be used as a benchmark for temporal precision of any unbiased estimator.As such, we set the pulse width to ensure that, on average, an estimate at the precisionof the CRB achieves a relatively high score. We set the benchmark metric score at0.8, as this represents a relatively high value in the range of the metric, which isbetween 0 and 1. The importance of this score is not the particular benchmark value— there are a range of values that give similar performance — but rather that it isa reproducible number with a clear interpretation. In this chapter, we characterisethe discrimination performance of CosMIC with a benchmark value of 0.8, so that itsscores can be interpreted when applied to spike inference algorithms on real data. Thebenchmark value was set lower than the metric’s maximum value, 1, so that the scoredoes not saturate when the model assumptions are not ideally satisfied. On real data,the noise may not be stationary (σ may vary), and so it may be possible for algorithmsto exceed the theoretical CRB. A benchmark score of 0.8 allows leeway to be exceededin this scenario.We consider a true spike at t1 and an estimate, U , normally distributed around it atthe precision of the CRB, such that U ∼ N (t1, σ2CRB). Then, we fix the pulse widthso that, on average, E [M(t1, U)] = 0.8. In Appendix B, we show that this conditionis satisfied when0.4 = (Φ(1/β)− 0.5)(β2 + 1)+β√2π(exp(−1/2β2)− 2), (4.58)where β = σCRB/w, w is the pulse width and Φ denotes the cumulative distributionfunction of the standard normal distribution. We observe empirically that the pulsewidth that solves this equation is approximately equal to 7.3σCRB (Fig. 4.5B).4.5 Experimental methodsIn Section 4.6, we compare CosMIC with the two most commonly used metrics inthe spike inference literature, the success rate and STC, which we define in Sections4.5.1 and 4.5.2, respectively. We fix the parameters of both metrics with respect toCosMIC’s pulse width, to ensure that the scores are comparable (Fig. 4.6). Themetrics are assessed on both simulated and real data. The latter dataset, which was85Chapter 4. CosMIC: A consistent metric for spike inferenceFigure 4.6: Illustration of common metrics. In this chapter, we compare thescores of three metrics: CosMIC, the spike train correlation (STC) and the successrate (SR). None of the metrics compute scores directly from the true and estimatedspike trains, x(t) and x̂(t), shown in A. Rather, CosMIC initially convolves the spiketrains with a triangular pulse (B). The STC first discretises the temporal interval andutilises the counts of spikes in each time bin, the bin edges and counts are plotted inC. The SR uses a bin centered around each true spike — an estimate in that bin isdeemed a true detection (D). In order that the metric scores are comparable, we fixthe STC and SR bin widths to be equal to CosMIC’s pulse width.collected by Dr T. Abrahamsson and Dr P. J. Sjöström, our collaborators at McGillUniversity, is described in Section 4.5.3.4.5.1 Success rateThe success rate, which is defined as a function of the true and false positive ratesor, alternatively, as a function of precision and recall, appears in various forms in theliterature. Spike inference performance has been assessed using true and false positiverates [106], precision and recall analysis [107] and using the complement of the successrate, the error rate [30]. We study this class of metrics under the umbrella of thesuccess rate, which we define here.A ground truth spike is deemed to have been ‘detected’ if there is an estimate withinδ1/2 (s) of that spike, where δ1 is a free parameter. Only one estimate can be deemedto detect one ground truth spike. The recall is the percentage of ground truth spikesthat were detected. The precision is the percentage of estimates that detect a ground864.5. Experimental methodstruth spike. Then, the success rate is the harmonic mean of the precision and recall,such thatsuccess rate = 2precision ∗ recallprecision + recall. (4.59)From this formula and Eq. (4.19), the relationship between the success rate and theSDC is apparent. A binary true detection region centred around each ground truthspike is analogous to an implementation of CosMIC with a box function pulse. Toensure that the success rate ‘pulse’ has the same width as CosMIC’s pulse, we setδ1 = w, where w is CosMIC’s pulse width.4.5.2 Spike train correlationThe first step in the calculation of the STC is the discretisation of the temporal in-terval into bins of width δ2. Two vectors of spike counts, c and ĉ, are subsequentlyproduced, whose ith elements equal the number of spikes in the ith time bin for thetrue and estimated spike trains, respectively. The STC is the Pearson product-momentcorrelation coefficient of the resulting vectors, i.e.STC =〈c−m(c), ĉ−m(ĉ)〉√v(c)√v(ĉ), (4.60)where 〈·, ·〉, m(·) and v(·), represent the inner product, sample mean and sample vari-ance, respectively. To remain consistent with the success rate, in all numerical exper-iments, we take δ2 = δ1 = w.The STC takes values in the range [-1, 1]. In practice, however, it is rare for aspike inference algorithm to produce an estimate that is negatively correlated with theground truth [7]. Moreover, an estimate with maximal negative correlation is equallyas informative as one with maximal positive correlation. We utilise the normalisedSTC, the absolute value of the STC. This ensures that the range of each metric thatwe analyse is equivalent (and equal to [0,1]) and that, as a consequence, the distributionof metric values are comparable.4.5.3 Two-photon calcium imaging of cortical slicesThis data was collected by Dr T. Abrahamsson and Dr P. J. Sjöström, our collaboratorsat McGill University. P11-P15 mice were anesthetized with isoflurane, decapitated, and87Chapter 4. CosMIC: A consistent metric for spike inferencethe brain was rapidly dissected in 4◦C external solution consisting of 125 mM NaCl,2.5 mM KCl, 1 mM MgCl2, 1.25 mM NaH2PO4, 2 mM CaCl2, 26 mM NaHCO3, and25 mM dextrose, bubbled with 95% O2/5% CO2 for oxygenation and pH. Quadruplewhole-cell recordings in acute visual cortex slices were carried out at 32◦C-34◦C withinternal solution consisting of 5 mM KCl, 115 mM K-gluconate, 10 mM K-HEPES, 4mM MgATP, 0.3 mM NaGTP, 10 mM Na-phosphocreatine, and 0.1% w/v biocytin,adjusted with KOH to pH 7.2-7.4. On the day of the experiment, 20 µM Alexa Fluor594 and 180 µM Fluo-5F pentapotassium salt (Invitrogen) were added to the internalsolution. Electrophysiology amplifier (Dagan Corporation BVC-700A) signals wererecorded with a National Instruments PCI-6229 board, using in-house software runningin Igor Pro 6 (WaveMetrics). The two-photon imaging workstation was custom builtas previously described [16]. Briefly, two-photon excitation was achieved by raster-scanning a Spectraphysics MaiTai BB Ti:Sa laser tuned to 820 nm across the sampleusing an Olympus 40x objective and galvanometric mirrors (Cambridge Technologies6215H, 3 mm, 1 ms/line, 256 lines). Substage photomultiplier tube signals (R3896,Hamatsu) were acquired with a National Instruments PCI-6110 board using ScanImage3.7 running in MATLAB (MathWorks). Layer-5 pyramidal cells were identified by theirprominent apical dendrites using infrared video Dodt contrast.4.6 ResultsTo investigate metric properties, we simulated estimated and ground truth spike trainsand analysed the metric scores. To mimic the temporal error in spike time estimation,unless otherwise stated, estimates were normally distributed about the true spike times.In the following, we refer to the standard deviation of the normal distribution as thejitter of the estimates.4.6.1 CosMIC penalises overestimationAs opposed to the STC, CosMIC and the success rate penalised overestimation ofspikes (Fig. 4.7). We simulated spike train estimates that were normally distributedabout the true spike times. When the number of detected spikes (K̂) was less than thenumber of true spikes (K), the locations about which the estimates were distributedwere chosen without replacement. When K̂ > K, the set of locations included allthe true spikes plus a subset of extras chosen with replacement. The overestimation884.6. ResultsFigure 4.7: Results on spike trains with varying overestimation ratio. Incontrast to the spike train correlation, CosMIC and the success rate were maximisedwhen the correct number of spikes were detected. We display the distribution of metricscores as the number of estimated spikes (K̂) varies with respect to the number of truespikes (K). The true spike train, which was identical throughout, consisted of 200spikes simulated from a Poisson process with spike rate 1Hz. Estimated spikes werenormally distributed about the true spikes, with jitter σCRB (A), 2 σCRB (B) and3 σCRB (C), respectively, where σCRB = 20ms. When the number of estimated spikeswas greater than the number of true spikes, estimates were distributed around a set oflocations including all true spikes plus an extra subset chosen with replacement. Foreach metric we plot the mean (darker central line) and standard deviation (edges ofshaded region) of metric scores on a set of 100 spike train estimates generated at eachoverestimation and jitter combination.ratio (K̂/K) reflects the degree of accuracy to which an estimate matches the rate ofa ground truth spike train. We observed that, rather than penalising overestimation,the STC increased with the overestimation ratio. In contrast, CosMIC and the successrate were maximised when the correct number of spikes were detected. This behaviourwas consistent as the jitter of the estimated spikes varied; in this example, the jitterwas σCRB (Fig. 4.7A), 2 σCRB (Fig. 4.7B) and 3 σCRB (Fig. 4.7C), respectively.It is the type of normalisation used by the STC that causes it to be insensitive tooverestimation. Scaling factors present in the spike count vectors cancel out in thenumerator and denominator, see Eq. (4.60), rendering the STC invariant under scalartransformations of the inputs. Even when the STC was adapted to the continuous-time89Chapter 4. CosMIC: A consistent metric for spike inferenceassessment of spike train similarity, by first convolving spike trains with a smoothingpulse, this flaw persisted [95].4.6.2 CosMIC rewards high temporal precisionCosMIC is more sensitive to temporal precision than the STC or success rate (Fig.4.8). First, we investigated this characteristic at the level of estimates of a singlespike, t1. CosMIC depends only on the absolute difference between the estimate, t̂1,and the true spike — the further the distance, the smaller the score. The relationshipbetween CosMIC and the temporal error, δ = t1 − t̂1, isM({t1}, {t̂1}) =(|δ|w − 1)2if |δ| < w0 otherwise,(4.61)where w is the width of the pulse. The derivation of this result is given in AppendixB. The success rate, on the other hand, does not reward increasing temporal precisionabove the bin width; an estimate is assigned a score of 1 or 0, when its precision isabove or below the bin width, respectively. Moreover, the STC is asymmetric in thetemporal error; estimates the same distance from the true spike are not guaranteed tobe awarded the same score, see Fig. 4.8A. This asymmetry stems from the correlation’stemporal discretisation. The temporal interval is first discretised into time bins andthe number of spikes in each bin are counted (Fig. 4.6). It follows that estimatedspikes that are the same absolute distance from a true spike can fall into different timebins, thus achieving a different score.On simulated data, we investigated the effect of these properties when spike trainestimates, rather than single spikes, were evaluated. In particular, we analysed themetric scores when spike train estimates contained the correct number of spikes buttheir temporal precision varied. We simulated the ground truth spike train as a Poissonprocess with rate 1Hz over 200s. The corresponding calcium transient signal wasgenerated assuming a Cal-520 pulse shape (see Table 2.1) and a sampling rate of 30Hz. White Gaussian noise was added to the calcium transient signal to generate twofluorescence signals, one with low and the other with relatively high noise (Fig. 4.8B).The corresponding metric pulse widths, as calculated from the CRB, were 33ms and78ms, or 1 and 2.3 sample widths, respectively. Spike train estimates were normallydistributed about the true spikes with varying jitter. The metric scores were then904.6. ResultsFigure 4.8: Results on spike train estimates with varying temporal preci-sion. CosMIC was more sensitive to the temporal precision of estimates than thespike train correlation (STC) or success rate (SR). Unlike the STC, CosMIC awardsestimated spikes (t̂1) with the same proximity to the true spike (t1) the same score (A).In contrast to both the STC and SR, CosMIC rewards increasing precision above thepulse width (2ε) with strictly increasing scores. In (C) and (D), we plot the distribu-tion of scores awarded to estimates that detect the correct number of spikes at varyingtemporal precision, in a low and high noise setting, respectively. In (B), a sample ofeach of the following signals are plotted: the ground truth spike train, simulated asa Poisson process at rate 1Hz over 200s; the corresponding calcium transient signal,sampled with interval T =1/30s; the low and high noise fluorescence signal and thecorresponding pulse widths. At each noise and jitter level, 100 realisations of spiketrain estimates normally distributed about the true spike times were generated. Inboth the low (C) and high noise (D) settings, the STC exhibited a relatively largevariation in the scores awarded to estimates of the same jitter. CosMIC and the SRwere roughly linearly related (E). CosMIC was boosted with respect to the success ratewhen temporal error, represented by the root mean square error (RMSE) of estimatesas a fraction of the pulse width, was low (F). Conversely, CosMIC was relatively lowwith respect to the SR when temporal error was relatively high. The colour map in (E)and (F) is thresholded at the 1st and 99th percentiles of the RMSE for visual clarity.91Chapter 4. CosMIC: A consistent metric for spike inferencecalculated for 100 realisations of spike train estimates at each jitter level in both thelow and high noise settings (Fig. 4.8C and D, respectively).As the correct number of spikes were always estimated, the level of jitter representedthe quality of a spike train estimate in this setting. Ideally, a metric would reliablyreward spike train estimates of the same quality with the same score. The STC,however, took a relatively large range of values for estimates of the same jitter (Fig.4.8C and D), despite having the same range as CosMIC and the success rate. Thisinconsistency is a consequence of the edge effects introduced by binning. Here, we usethe term consistency in line with its semantic rather than mathematical definition.We observed a roughly linear trend in the scores of CosMIC and the success rate (Fig.4.8E). As expected, CosMIC was boosted with respect to the success rate when theroot mean square error (RMSE) of detected spikes was relatively low when measured asa fraction of the pulse width. In each case, the RMSE was computed empirically fromthe estimated spikes within the precision of the pulse width. Conversely, CosMIC wasrelatively low with respect to the success rate when the RMSE was relatively high. Thistrend is visible in the Bland-Altman plot [3], in which the mean of the two methods isplotted against the difference (Fig. 4.8F). We conclude that CosMIC is more sensitiveto the temporal precision of detected spikes, as, unlike the success rate, it discriminatesprecision above the bin width.4.6.3 Application to real imaging dataOn imaging data of the mouse visual cortex at a frame rate of 13 Hz, CosMIC wasmore sensitive than the success rate to the temporal precision of detected spikes. For adetailed description of the imaging data, see Section 4.5.3. As detailed in Section 4.4,the metric’s pulse width was set with respect to the CRB. On this dataset, the pulsewidths were concentrated between 1 and 3 sample widths — this range encompassed92% of the data, see (Fig. 4.9F). As the noise level of the data increases, so does thepulse width. Consequently, the tolerance of the metric with regards to the temporalprecision of estimates also increases. As a result, estimates on noisier data (Fig. 4.9B)were scored with more lenience than those on less noisy data (Fig. 4.9A). Spikeswere detected from this data using the version of the FRI spike inference algorithmpresented in Chapter 5 of this thesis.As was found on simulated data in Section 4.6.2, there was a linear trend between924.6. ResultsFigure 4.9: Results on real imaging data. On mouse in vitro imaging data,CosMIC was more sensitive than the success rate (SR) to the temporal precision ofdetected spikes. Spikes were detected using the FRI algorithm from 83 traces sampledfrom visual cortex slices at 13Hz. In A and B, we display from top to bottom: anexample fluorescence trace (∆F/F ), ground truth and detected spike trains, and thecorresponding pulse trains. There was an approximately linear relationship betweenCosMIC and the SR (C). CosMIC was relatively high with respect to the SR whentemporal error, represented by root mean square error (RMSE) as a fraction of the pulsewidth, was relatively low. Conversely, CosMIC was low with respect to the SR whentemporal error was relatively high. This pattern was conserved in the relationshipbetween the precision and CosMIC’s analogous ancestor metric, PCosMIC, (D) andbetween the recall and RCosMIC (E). The range of pulse widths as computed from theCramér-Rao bound (F) and the range of spike train correlation (STC) scores (G) arealso shown.93Chapter 4. CosMIC: A consistent metric for spike inferencethe scores of CosMIC and the success rate (Fig. 4.9C). CosMIC was relatively highwith respect to the success rate when the temporal error represented by RMSE asa fraction of the pulse width, was relatively low. Conversely, CosMIC was low withrespect to the success rate when the temporal error was relatively high. The RMSEwas computed empirically from the estimated spikes within the precision of the pulsewidth. This pattern was conserved when CosMIC’s ancestor metrics, PCosMIC andRCosMIC (Section 4.3.1), were compared to the precision and recall (Fig. 4.9D and E).The average RMSE over all traces was 27ms, or 0.37 sample widths. As CosMIC isable to discriminate precision above the pulse width, it is more able to reward thissuper-resolution performance than the success rate or STC.4.6.4 CosMIC discriminates precision and fall-out rate of spike trainsBy construction, CosMIC bears a strong resemblance to the Sørensen-Dice coefficientand, consequently, the success rate. Both metrics can be written as the harmonicmean of the precision and recall, two intuitive metrics which represent the proportionof estimates that detect a ground truth spike and the proportion of true spikes detected,respectively. In this section, we demonstrate that CosMIC can accurately discriminateboth the precision and recall of spike train estimates.When a spike train estimate detects exactly a subset of the true spikes, plus no re-mainders, CosMIC and the success rate depend only on the percentage of true spikesdetected (the recall) and not the location of that subset, see Fig. 4.10A and 4.10D.Denoting the size of the subset of true detections as K − R, with K the number oftrue spikes and 0 < R ≤ K, we haveM(S, Ŝ)= 1− 12K/R− 1, (4.62)see Appendix B for a proof. Thus, CosMIC depends only on the proportion of ‘missing’spikes, R/K, not their location. In contrast, the STC exhibits significant variation ateach level of recall. This is illustrated in Fig. 4.10A, in which we plot the distributionof CosMIC, success rate and correlation scores over 100 realizations of spike trainestimates at each level of recall. It can be seen that, in this setting, CosMIC and thesuccess rate are fixed with the recall of the spike train estimates.When all the true spikes were exactly detected together with R > 0 surplus spikes,CosMIC and the success rate depend only on the level of precision not the location of944.7. Summarythe surplus spikes, see Fig. 4.10B and 4.10E. We haveM(S, Ŝ) =11 +R/2K, (4.63)where K is the number of true spikes, see Appendix B for a proof. The fall-out rate,which is the complement of the precision, is the proportion of estimates that were notdeemed to have detected a ground truth spike. It is apparent from Eq. (4.63) that,in this setting, CosMIC depends only on the fall-out rate, R/K. The correlation, onthe other hand, varied with the location of the surplus spikes. In Fig. 4.10E, we plotthe distribution of the correlation scores for 100 realizations of spike train estimatesat each level of precision. CosMIC and the success rate, which were constant (andidentical) at a given precision, in this scenario, are also shown.4.7 SummaryIn this chapter, we presented a new metric, CosMIC, with which to assess the similarityof spikes inferred from calcium imaging data to the ground truth spiking activity.Rather than operating on the true and estimated spike trains directly, the proposedmetric assesses the similarity of the pulse trains obtained from convolution of the spiketrains with a smoothing pulse. The pulse width is derived from the Cramér-Rao bound,which, in turn, is computed from the statistics of the imaging data. As a result, themetric quantifies the accuracy of the inferred spikes with respect to the limitations ofthe data. The final metric score is the size of the commonalities of the pulse trains asa fraction of their average size.We demonstrated the advantage of CosMIC over the two most commonly used metrics,the STC and the success rate. Unlike the STC, which appears to reward overestima-tion, CosMIC is maximised when the correct number of spikes have been detected.Moreover, CosMIC is more sensitive to the temporal precision of estimated spikesthan the success rate. We consider a single summary score to be practical for userswho do not have the time or desire to analyse multi-dimensional trade-offs. Alterna-tively, CosMIC’s ancestor metrics, RCosMIC and PCosMIC, can be used to determine theextent to which errors stem from undetected or falsely-detected spikes.95Chapter 4. CosMIC: A consistent metric for spike inferenceFigure 4.10: Results on spike train estimates with varying precision andfall-out rate. CosMIC scored estimated spike trains of the same recall and fall-outrate consistently, unlike the spike train correlation (STC). When a spike train estimatedetected precisely the location of a subset of spikes from a true spike train, the scoresof CosMIC and the success rate depended only on the percentage of spikes detected(the recall), not the location of the detected spikes (A, D). In contrast, the STC variedwith the subset of spikes that were detected. When a spike train estimate detectedall the true spikes precisely plus a number of surplus spikes, the STC varied withthe placement of the surplus spikes (B, E). In contrast, the success rate and CosMICdepended only on the percentage of estimated spikes that did not correspond to groundtruth spikes (the fall-out rate, also known as the false positive rate). The distributionof correlation scores plotted in D and E stem from 100 realizations of estimated spiketrains at each recall and fall-out rate. In C, we plot an example of a true spike train.In D and E, we plot estimated spike trains, with a recall and fall-out rate of 50% and33%, respectively, along with the corresponding metric scores. The spikes with a black‘x’ marker in E indicate the surplus spikes.96Chapter 5.Spike inference using FRI theoryIn two-photon calcium imaging, when a neuron fires a spike, it also emits a character-istic pulse of fluorescence. The fluorescence signal from one neuron over time is thusmodelled as a stream of pulses, each initiated by a spike. Neuronal network activityis analysed using spike times inferred from the imaging data. Due to the challengesinherent in imaging through living tissue, this data is typically noisy. Furthermore,sampling rates are relatively low (on the order of 10Hz) with respect to the temporalprecision required for network connectivity analysis. Robust signal processing meth-ods that are able to infer spike times from noisy data with precision above the samplewidth are thus required.In previous work, it was identified that a neuron’s fluorescence signal belongs to theclass of signals with finite rate of innovation (FRI). Using FRI theory and a simplifiedmodel of the pulse shape, spike times were inferred with high accuracy and temporalprecision. The simplified pulse model is, however, not ideal for a subset of fluorescentindicators, including some of the most promising, genetically-encoded indicators. Inthis chapter, we extend the FRI framework to encompass a wider class of pulse shapes.We also introduce strategies to increase the robustness of the algorithm, includinga procedure to whiten the statistics of the noise and an updated approach to theestimation of the number of spikes in a window. We demonstrate on simulated datathat the updated signal model eliminates estimation bias due to model mismatch.Additionally, we compare the performance of the updated FRI algorithm and a stateof the art deconvolution algorithm on mouse in vitro imaging data.The work presented in this chapter led to the following publication [109].97Chapter 5. Spike inference using Finite Rate of Innovation theory5.1 IntroductionThe goal of spike inference is to recover a neuron’s spiking activity from its fluorescencesignal (Fig. 5.1). To this end, several authors proposed supervised learning algorithmsthat were trained on datasets with ground truth spiking activity [114, 135, 127]. De-spite promising results on familiar data, it has been shown that state of the art su-pervised learning approaches impose spurious auto-correlation patterns onto test data[93]. This type of error, which occurs when patterns of spiking activity in the trainingdata are implicitly learned by an algorithm, could precipitate misleading conclusionsabout neuronal activity. Until the nature of errors of supervised learning algorithmsare fully characterised, it is likely that model-based algorithms will be predominatelyused.In contrast to supervised learning approaches, model-based algorithms take advan-tage of the known structure of the fluorescence signal to perform spike inference[153, 148, 47, 87, 103, 30, 106, 93]. These algorithms can be divided into two groups,according to whether their model is in discrete or continuous time. The former groupmodel the signal as an auto-regressive process, whereby fluorescence intensity at a givensample depends on past fluorescence values and the (unknown) number of spikes thatoccurred in that sample [148, 147, 103, 105, 44]. With the exception of a Bayesianalgorithm [103], these approaches do not infer the locations of spikes in continuoustime. Rather, the output is a vector of equal length to the fluorescence signal, whereeach entry reflects the number of spikes fired in that sample. Vogelstein et al. arguedthat identification of the most likely vector of spiking activity is intractable when itsentries are integer-valued [147]. Instead, they relaxed their problem to infer a vectorwith non-negative, real entries, whereby each entry loosely represents the probabilitythat a spike is fired in that sample. The work of Vogelstein et al. inspired a family ofnon-negative deconvolution methods that infer vectors of spiking probabilities undercertain sparsity constraints [105, 44]. When the SNR of traces is low and, consequently,temporal precision is not likely to surpass the sample width, non-negative deconvolu-tion approaches can perform well. They are computationally relatively efficient andconvey uncertainty in the spike time estimates. However, as temporal resolution isimplicitly limited in this framework, they are not ideal when the end goal of an ex-periment is, for example, to analyse temporal coding hypotheses and causal networkactivity.985.1. IntroductionFigure 5.1: Example of a slow-rise fluorescence signal. The signal from oneneuron is modelled as a stream of pulses, each initiated by a spike. This figure showsdata sampled at 60Hz generated using the genetically encoded indicator GCaMP6s[132, 26]. In B, we plot a single calcium transient from the signal in A. The rise of thisindicator’s pulse is relatively slow — the time from baseline to peak is approximately300ms.Continuous-time models used for spike inference vary in the degree of detail with whichthey model calcium dynamics. Rahmati et al. formulated a detailed biophysical modelrelating the fluorescence signal to neuronal membrane potential, considering propertiesof ion channels and sub-threshold voltage fluctuations [106]. This approach, whilepotentially cumbersome for bulk image analysis, allows inference of a cell’s biophysicalproperties as well as spiking activity. As such, it could be fruitful when used incombination with pharmacological interventions, in which the aim is to alter thoseproperties [114]. Typically, continuous-time models use a simpler framework, in whichthe fluorescence signal is a convolution of a spike train with a characteristic pulse[47, 87, 30]. A variety of methods have been used to infer spike times in this setting,including template matching [47] and maximum likelihood estimation [30]. Of thesetwo, the latter was found to have the highest accuracy on a wide array of real datawith electrophysiological ground truth [30].Oñativia et al. inferred the timing of spikes by exploiting the fact that the fluorescencesignal belongs to the class of signals with finite rate of innovation [87]. A signalbelongs to this class if it admits a parametric representation with a finite number offree parameters per unit of time [145]. In the case of neuronal fluorescence signals,the free parameters are the unknown spike times and amplitudes. By mapping theproblem of spike inference to the FRI problem of reconstructing a stream of Diracs,Oñativia et al. were able to utilise well-established methods to detect spike locations.The mapping incorporated assumptions about the fluorescence pulse, which Oñativiaet al. modelled with an instantaneous rise and exponential decay. This simplification99Chapter 5. Spike inference using Finite Rate of Innovation theoryis common [147, 30] and well-suited to situations in which the pulse rise time is notsignificantly longer than the sampling period. However, scan rates are consistentlyincreasing [117] and some of the most promising fluorescent indicators have relativelyslow rise times [26]. For example, the genetically-encoded indicator GCaMP6s, whichhas a high dynamic range and brightness, rises over a period of 200-300ms (Fig. 5.1).In these cases, there is a risk of estimation bias due to model mismatch.In this chapter, we extend the FRI framework for spike inference from calcium imagingdata to encompass a slow-rise pulse shape. In the FRI framework, the sampling processof the continuous-time spike train is modelled explicitly, so we start by introducingrelevant aspects of sampling theory in Section 5.2. In order to increase the algorithm’srobustness to noise, we introduce a pre-whitening step [141], so that the statistics ofthe noise are suited to the parameter estimation techniques used. Furthermore, weinclude a model-based procedure to estimate the number of spikes in a window. Wepresent these contributions in Section 5.3. Finally, we assess the performance of theupdated algorithm on real and simulated data.5.2 Sampling theoryLet x(t) be a continuous-time signal, which is sampled by a generic acquisition device.The modifications that occur to the signal during the acquisition process are modelledby convolution of x(t) with a filter, h(t). In digital cameras, for example, this representsthe blurring of the original image by a low-quality lens. The actual signal from whichsamples are obtained is thusy(t) = x(t) ∗ h(t) =∫Rx(τ)h(t− τ) dτ. (5.1)Discrete samples, which are obtained with uniform time resolution, T , are denotedy[n] = y(t)|t=nT =∫Rx(τ)h(nT − τ) dτ, (5.2)for n ∈ {0, 1, ..., N − 1}. Defining h(t) = ϕ(− tT ) and using the inner product notation,the samples are expressed asy[n] =∫Rx(τ)ϕ( τT − n) dτ =〈x(t), ϕ(tT − n)〉, (5.3)1005.2. Sampling theorywhere we refer to ϕ as the sampling kernel.Whereas in some applications, such as standard digital photography, it can be suffi-cient to retain the digital samples, it is typically desirable to use the samples to recoveraspects of the original, continuous-time signal, x(t). Indeed, under certain conditionson the signal and filter, it is possible to perfectly reconstruct the original signal fromdigital samples. In the mid 20th century, it was established that a bandlimited sig-nal can be completely determined by sufficiently finely-spaced samples [119, 152, 67].Whilst an important development, the corresponding reconstruction formula is onlyapproximately applicable to real-world signals, which are limited in time and thusnot truly bandlimited [121]. In 2002, Vetterli et al. established conditions for per-fect reconstruction of a wider class of signals [145]. In the following, we introducethis framework, which will later be used to infer the locations of spikes from calciumimaging data.5.2.1 Signals with finite rate of innovationWe consider a class of signals that can be represented as a linear combination of shiftsof known functions, {gr(t)}Rr=0, such thatx(t) =∑k∈ZR∑r=1ar,kgr(t− tk), (5.4)where tk and ar,k are free parameters for r ∈ {1, ..., R} and k ∈ Z. In this case, theproblem of reconstructing x(t) is equivalent to that of retrieving the unknown param-eters. When the parameters are known, x(t) is defined everywhere. The complexity ofthe reconstruction problem is thus reflected by the number of free parameters per unitof time. To formalise this concept, a counting function Cx(ta, tb) is used that countsthe number of free parameters in x(t) over the temporal interval [ta, tb]. Then, the rateof innovation of x(t) is defined asρ = limτ→∞1τCx(− τ2 ,τ2 ). (5.5)A signal has finite rate of innovation if it admits a parametric representation as inEq. (5.4) with a finite ρ. This class of signals, which includes bandlimited signals,encompasses a wide range of physical processes [139, 35, 84]. FRI methods exploitthe knowledge of these signals’ parametric structure in order to retrieve the unknown101Chapter 5. Spike inference using Finite Rate of Innovation theoryFigure 5.2: Example of signals with finite rate of innovation. Both signalsare linear combinations of a function with varying shifts and amplitudes. In A, theconstructive function is smoothly varying whereas in B it is a point impulse.parameters and, hence, reconstruct the signal.An example of an FRI signal, which is a linear combination of shifts of a single function,g(t), is illustrated in Figure 5.2A. In this case, x(t) consists of a finite number, K, ofshifts of g(t) and can be writtenx(t) =K∑k=1akg(t− tk). (5.6)A further example, which is used to model neuronal spike trains, is a stream of pointimpulses (Fig. 5.2B). A point impulse is expressed mathematically using the DiracDelta function, δ(t), which is zero at all locations on the real line except at t = 0. TheDelta function is defined with respect to its effect on the integral of any continuousfunction, s(t), such that [75] ∫Rs(t)δ(t) dt = s(0). (5.7)A stream of Diracs at locations {tk}Kk=1 with amplitudes {ak}Kk=1 is writtenx(t) =K∑k=1akδ(t− tk). (5.8)Due to its connection with our neuroscience problem, in the following, we focus onhow FRI theory can be used to sample and reconstruct streams of Diracs.5.2.2 Sampling kernelsThe reconstruction scheme is tailored to the sampling kernel used in the acquisitionprocess. Among others, reconstruction schemes have been proposed for sinc and Gaus-sian kernels [145], kernels with rational Fourier transform [38] and polynomial and1025.2. Sampling theoryFigure 5.3: Exponential reproduction with an E-spline of order 2. In A,we plot the E-spline in the time domain. In B and C, the exponentials eα0t and eα1tare reproduced, respectively, over the interval [1, 10].exponential reproducing kernels [38]. Here, we focus on the latter class of kernels.An exponential reproducing kernel is a function that, when linearly combined with itsuniformly shifted versions, can reproduce exponentials. More formally, it is a function,ϕ(t), for which there exists a set of coefficients, cm,n ∈ C, and exponents, αm ∈ C,such that ∑n∈Zcm,nϕ(t− n) = eαmt for m ∈ {0, ..., P}, (5.9)almost everywhere. A function reproduces exponentials if and only if it satisfies thegeneralised Strang-Fix conditions [149]:L{ϕ} (αm) 6= 0 and L{ϕ} (αm + 2πil) = 0, (5.10)for m ∈ {0, ..., P} and l ∈ Z \ 0, where L{ϕ} represents the double-sided Laplacetransform of ϕ(t).A family of functions that satisfy the exponential reproduction property are the ex-ponential B-splines [140], otherwise known as E-splines, whose name refers to theirrelation to polynomial-reproducing B-splines. A zero order E-spline is defined as fol-lows:βα(t) =eαt 0 ≤ t ≤ 1,0 otherwise,F−→ β̂α(ω) =1− eα−iωiω − α. (5.11)In the case of the zero-order E-spline, the parameter α indicates the exponent of the103Chapter 5. Spike inference using Finite Rate of Innovation theoryexponential that can be reproduced. More generally, an E-spline of order P ≥ 0 isassociated with a vector of exponents α = {α0, α1, ..., αP } and is obtained throughconvolution of zero order E-splines, such thatβα(t) = βα0(t) ∗ βα1(t) ∗ ... ∗ βαP (t)F−→ β̂α(ω) =P∏m=01− eαm−iωiω − αm. (5.12)This E-spline can reproduce functions in the subspace spanned by {eα0t, eα1t, ..., eαP t},as long as no two exponents are separated by a multiple of 2πi, which would violate theStrang-Fix conditions. In Fig. 5.3, we illustrate the exponential reproduction propertyfor an E-spline of order two. To reproduce the exponentials for all t ∈ R would requirethe sum in Eq. (5.9) to be evaluated over all n ∈ Z. In this example, we truncate thesum and reproduce the exponentials over the interval [1, 10].5.2.3 Reconstruction of a stream of DiracsWe now consider the situation in which a stream of Diracs, x(t), has been sampledwith an exponential reproducing kernel, h(t) = ϕ(− tT ), at time resolution T to producesamples {y[n]}N−1n=0 . To retrieve the unknown spike times and amplitudes, the samplesare linearly combined using the exponential reproduction coefficients of the samplingkernel, cm,n, such thats[m] =∑n∈Zcm,ny[n] =∑n∈Zcm,n〈x(t), ϕ(tT − n) 〉=〈x(t),∑n∈Zcm,nϕ(tT − n) 〉, (5.13)for m ∈ {0, 1, ..., P}, where the last step follows from the linearity of the inner product.We refer to s[m] as the sample moments. Due to the exponential reproduction property,see Eq. (5.9), it follows thats[m] =〈x(t), eαmt/T〉=∫Rx(t)eαmt/T dt =K∑k=1akeαmtk/T , (5.14)where the last step follows from Eq. (5.7). When the exponents are evenly spaced,that is αm = α0 + λm, the sample moments can be written in power-sum form, suchthats[m] =K∑k=1bkumk , (5.15)with bk = akeα0tk/T and uk = eλtk/T .1045.2. Sampling theoryFigure 5.4: Example of sampling and reconstruction of a stream of Diracs.A stream of 3 Diracs, x(t), is filtered with an exponential reproducing kernel, h(t) =ϕ(−t/T ), to produce y(t) = x(t) ∗ h(t). Samples y[n] are obtained at time resolutionT . The spike times and amplitudes are recovered from the samples y[n].The formation of the sample moments can also be written in matrix vector form, suchthat s = Cy, where y and s represent the vectors of samples and sample moments, re-spectively, and the (i, j)th entry of C is ci,j . We select αm to ensure that the coefficientmatrix C is well-conditioned. Urigüen et al. demonstrated that E-splines withαm =πP+1 (2m− P ) i and cm,n = eαmn (ϕ̂(αm))−1 , (5.16)enabled relatively stable reconstruction of streams of Diracs [141]. This is the formu-lation we use in the rest of this chapter. In Fig. 5.4, we display an example of suchan E-spline. We note that the kernel itself is real, as the exponents occur in complexconjugate pairs.In order to retrieve the spike times and amplitudes, one needs to recover {bk}Kk=1 and{uk}Kk=1 from the sample moments. This is a common problem in spectral estimationtheory, for which there are many methods [129]. For example, Prony’s method can beused to identify the filter h[m] that annihilates s[m], such thath[m] ∗ s[m] =K∑l=0h[l]s[m− l] = 0. (5.17)One can show that the filter h[m] that satisfies Eq. (5.17) has roots that correspondto the parameters uk [38]. Thus, the parameters uk are obtained as the roots of theannihilating filter. Using Prony’s method it is possible to perfectly retrieve uk, andthus tk, using the minimum number of measurements, P + 1 = 2K, see Fig. 5.4.In practice, the samples and, consequently, the sample moments are corrupted by105Chapter 5. Spike inference using Finite Rate of Innovation theorynoise. In this case, it is preferable to use P ≥ 2K moments in order to increase therobustness of the parameter estimation procedure. One can also use Prony’s methodwith denoising strategies, such as Cadzow denoising [19]. In this thesis, we use therobustified matrix pencil method [53, 54], which is as effective as Cadzow denoisingwith Prony’s method but has the advantage that it is not iterative.Matrix pencil methodThe matrix pencil method identifies uk using the Toeplitz matrix of sample moments:S =s[M ] s[M − 1] ... s[0]s[M + 1] s[M ] ... s[1]....... . ....s[P ] s[P − 1] ... s[P −M ] , (5.18)where M = bP/2c. First, the singular value decomposition of this matrix is computed,such that S = UΣVH . Here, Σ is diagonal and U and V contain the left and rightsingular vectors of S, respectively. A truncated matrix UK is subsequently formedfrom the left singular vectors corresponding to the K largest singular values. Whenthe noise level is sufficiently low, these singular vectors form a basis of the signalsubspace. Defining S0 and S1 as the matrix UK with the first and last row removed,respectively, the parameters {uk}Kk=1 are identified as the eigenvalues that solve thefollowing generalised eigenvalue problem(S†1S0 − µI)v = 0, (5.19)where S†1 denotes the Moore-Penrose psuedoinverse of S1.5.3 Spike inference from calcium imaging dataIn the following, we outline how the problem of spike inference from calcium imagingdata is mapped to the FRI problem of reconstructing a stream of Diracs. We start byre-introducing the mathematical model of a neuronal fluorescence signal. For a moredetailed discussion on the origins of the signal model, we refer the reader to Section2.3.1065.3. Spike inference from calcium imaging dataWhen a spike is fired from a neuron, a characteristic pulse of fluorescence is emitted,such thatp(t) = c(1− e−t/τon)e−t/τoff1t>0, (5.20)where c is a normalisation constant that ensures that max {p(t)} = 1 and τon, τoff areparameters defining the speed of the rise and decay, respectively. The parameters τonand τoff are largely determined by the calcium indicator that was used to generate thedata (see Section 2.3) and are straightforward to estimate [148, 103, 105, 30]. In thefollowing, therefore, we consider these to be known parameters. A neuron’s fluorescencesignal is modelled as a stream of pulses at spike times {tk}Kk=1 with amplitudes {ak}Kk=1,such thatf(t) =K∑k=1akp(t− tk), (5.21)where we assume the baseline fluorescence has been subtracted. There are manyprocedures for removing the baseline component, typically a parametric model is fit tothe signal and subsequently subtracted [30, 106, 44, 104].For the moment, we consider the number of spikes (the model order) to be known. Wediscuss the estimation of this value in Section 5.4.2. The fluorescence signal, therefore,has 2K unknown parameters — the spike times and amplitudes — and also belongsto the class of FRI signals. Denoting x(t) =∑Kk=1 akδ(t− tk), we havef(t) =K∑k=1akp(t− tk) = p(t) ∗K∑k=1akδ(t− tk) = p(t) ∗ x(t). (5.22)We assume that the fluorescence signal has been filtered with an exponential repro-ducing kernel, h(t) = ϕ(− tT ) and sampled with time resolution T , such that we havesamplesy[n] =〈f(t), ϕ(tT − n)〉=〈x(t) ∗ p(t), ϕ(tT − n)〉, (5.23)107Chapter 5. Spike inference using Finite Rate of Innovation theoryfor n ∈ {0, ..., N}. Then, expanding the integrals,y[n] =∫R[∫Rx(τ) p(t− τ) dτ]ϕ( tT − n) dt (5.24)=∫R[∫Rϕ( tT − n) p(t− τ) dt]x(τ) dτ (5.25)=∫Rx(τ)(p(−τ) ∗ ϕ( τT − n))dτ (5.26)=〈x(t), p(−t) ∗ ϕ(tT − n)〉, (5.27)where, in Eq. (5.25), we are able to exchange the order of integration as y[n] is finite(the integrand is bounded and has compact support).Practical implementationSo far, we have assumed that the fluorescence signal was sampled with an exponentialreproducing kernel. While not impossible, two-photon microscopes do not currentlyimplement such a sampling kernel. In practice, we have access to f [n] = f(nT )|t=nT .We compute the samples y[n] through discrete convolution of f [n] and discretisedsampling kernel, such thaty[n] =∑m∈Zf [n]h[n−m]. (5.28)To enable the theoretical analysis, in the following we use the continuous-time for-mulation. The discrete implementation of the sampling kernel has been shown to beeffective in a range of setups [88], and has been tested on real data from a range ofindicators with sampling rate ranging from 8Hz [87] to 60Hz [109].5.3.1 Mapping samples to the FRI settingIn the following proposition, we establish the filtering operations that transform sam-ples of slow-rise calcium transients into samples of the underlying spike train filteredwith an exponential reproducing kernel. We define a discrete filter, g[n], such that,when convolved with y[n], we obtainw[n] = y[n] ∗ g[n] =〈x(t), ψ( tT − n)〉. (5.29)1085.3. Spike inference from calcium imaging dataFigure 5.5: Example of spike inference with the FRI algorithm. A fluo-rescence signal, f(t), is sampled with kernel ϕ(−t/T ) to produce samples y[n]. Thesamples are subsequently convolved with a discrete-time filter to produce w[n], fromwhich the unknown spike times can be retrieved.The analogous filter was identified under the assumption of pulses with an instanta-neous rise and exponential decay in [87]. In that case, g[n] had two non-zero coefficientswhose value depended on the decay parameter of the pulse. The resulting kernel, ψ(t),was a convolution of the original sampling kernel, ϕ(t), and a zero order E-spline.In our setting, assuming the samples are as in Eq. (5.27), g[n] has three non-zerocoefficients andψ(t) = βατ (− tT ) ∗ ϕ(t), (5.30)where ατ = {−T/τoff, −T/τoff − T/τon}. The coefficients of both the filter and theE-spline depend on the pulse parameters, τon and τoff. As the exponential reproductionproperty is conserved through convolution [140], ψ(t) is also an exponential reproducingkernel and, in particular, reproduces the exponentials that are reproduced by ϕ(t).Although the exponents are conserved, the reproduction coefficients are, in general,different. We compute the sample moments using the reproduction coefficients of ψ(t),such thats[m] =∑n∈Zdm,nw[n]. (5.31)Then, the methods of Section 5.2.3 can be used to reconstruct the unknown spike timesand amplitudes (see Fig. 5.5).Proposition 5.1. Samples y[n] are convolved with filter g[n] = τonT g1[n]∗g2[n], whereg1[n] =δ[n]− exp(− Tτoff)δ[n− 1], (5.32)g2[n] =δ[n]− exp(− Tτoff −Tτon)δ[n− 1], (5.33)109Chapter 5. Spike inference using Finite Rate of Innovation theoryand δ[n] indicates the discrete delta function. This producesw[n] = y[n] ∗ g[n] =〈x(t), βατ (−t/T ) ∗ ϕ( tT − n)〉, (5.34)with ατ = {−T/τoff, −T/τoff − T/τon}.Proof. From the definition of w[n], we havew[n] =∑m∈Zg[m]y[n−m]. (5.35)Substituting the value of y[n] from Eq. (5.27) and using the linearity of the innerproduct and convolution, w[n] becomesw[n] =〈x(t), p(−t) ∗∑m∈Zg[m]ϕ(tT − (n−m))〉. (5.36)Representing the filter coefficients in continuous time with g(t) =∑m∈Z g[m]δ(t−mT ),we haveg(−t) ∗ ϕ( tT − n) =∫Rg(−τ)ϕ( t−τT − n) dt (5.37)=∫R∑m∈Zg[m]δ(−τ −mT ) ϕ( t−τT − n) dt (5.38)=∑m∈Zg[m]ϕ( tT − (n−m)). (5.39)It follows thatw[n] =〈x(t), p(−t) ∗ g(−t) ∗ ϕ( tT − n)〉. (5.40)Decomposing the continuous-time Fourier transform of g(−t) into the discrete-timeFourier transforms of filters g1[n] and g2[n], we haveF {g(−t)} (ω) =∑m∈Zg[m]eiωmT =τonTG1(e−iωT )G2(e−iωT ). (5.41)1105.3. Spike inference from calcium imaging dataWriting c1 = − 1τoff and c2 = −1τoff− 1τon , this becomesF {g(−t)} = 1T (c1 − c2)(1− ec1T eiωT) (1− ec2T eiωT)(5.42)=1T (c1 − c2)(1− ec1T+iωT) (1− ec2T+iωT). (5.43)Noting thatF {p(−t)} (ω) = c1 − c2(c1 + iω)(c2 + iω), (5.44)and denoting ατ = {α1, α2} = {c1T, c2T}, we haveF {p(−t) ∗ g(t)} = F {p(−t)}F {g(t)} (5.45)=1T2∏m=11− ecmT+iωTcm + iω(5.46)= T2∏m=11− eαm+iωT−αm − iωT(5.47)= F {βατ (−t/T )} . (5.48)Finally, it follows thatw[n] =〈x(t), βατ (− tT ) ∗ ϕ(tT − n)〉. (5.49)5.3.2 Amplitude estimationThe fluorescence signal has 2K unknown parameters: K spikes times and K ampli-tudes. Accurate computation of the amplitudes, which are of limited practical interest,is essential for the procedure that we use to estimate the model order (Section 5.4.2).To estimate the amplitudes, we use the following matrix-vector formulation of thefluorescence signal model:f [0]f [1]...f [N − 1] =p1[0] p2[0] . . . pK [0]p1[1] p2[1] . . . pK [1]....... . ....p1[N − 1] p2[N − 1] . . . pK [N − 1]a1a2...aK , (5.50)111Chapter 5. Spike inference using Finite Rate of Innovation theorywhere pj [n] = p(t − tj)|t=nT . Writing Eq. (5.50) more compactly as f = Pa, weestimate the amplitude vector asâ = mina∥∥Pa− f∥∥2 where 12A ≤ â ≤ 32A. (5.51)Here, we constrain the entries of â to be within a certain range of physiologicallyplausible amplitudes, where A is the expected amplitude of a calcium transient.5.4 Spike inference from noisy dataIn practice, sampling a signal with any acquisition device generates measurement noise.We model the samples as corrupted by additive noise, such thatỹ[n] = y[n] + εy[n], (5.52)where εy[n] are i.i.d. noise samples with distribution N (0, σ2). As discussed in Section2.4.4, noise generated in fluorescence microscopy is predominantly photon shot noise,which is typically modelled in this way [147, 105, 106, 44, 127]. As noted earlier, inpractice, we obtain ỹ[n] from the raw signal samples by applying the convolution withan exponential reproducing kernel in post-processing. The model of the noise in Eq.(5.52), therefore, only approximately holds.The filtering operations that are used to map the 2PCI spike inference problem to thatof reconstructing a stream of Diracs cause the noise to become correlated. We havew̃[n] =∑m∈Zg[m]ỹ[n−m] (5.53)=∑m∈Zg[m]y[n−m] +∑m∈Zg[m]εy[n−m] (5.54)= w[n] + εw[n], (5.55)where the noise samples, εw[n], are still Gaussian but are no longer independent. Whenwe compute the sample moments, the noise distribution is further altered, such thats̃[m] =∑n∈Zdm,nw̃[n] =∑n∈Zdm,nw[n] +∑n∈Zdm,nεw[n] = sm + εs[m], (5.56)1125.4. Spike inference from noisy datawhereεs[m] =∑n∈Z∑k∈Zdm,ng[k]εy[n− k]. (5.57)Finally, the Toeplitz matrix of noisy sample moments, which is analogous to Eq. (5.18),is denoted S̃ = S + B, where B contains the contributions from the noise.When the data is not corrupted by noise, the sample moment matrix has rank K.Determining the number of innovations is therefore straightforward. In the noisysetting, S̃ is no longer rank-deficient and the estimation of the number of innovationsis more challenging. In Section 5.4.2, we outline the method we use to estimate thisparameter. Additionally, the performance of the methods that infer spike times fromthe noisy sample moment matrix are negatively affected by correlated noise. In Section5.4.1, we discuss a method to whiten the statistics of the noise.5.4.1 Pre-whiteningThe matrix pencil method, which we use to estimate the spike times from the samplemoment matrix, aims to separate the signal subspace from the noise subspace. To doso, it assumes that the singular vectors corresponding to the largest K singular valuescharacterise the signal subspace. When the noise is correlated, however, it is likelythat some of the dominant singular vectors characterise trends in the noise rather thanthe signal. In this case, some singular vectors corresponding to the signal subspacewill be discarded and the unknown parameters will not all be retrieved correctly.To mitigate this problem, a pre-whitening transform can be applied to S̃ [141], suchthatS̃′ = S̃W = (S + B) W = SW + BW, (5.58)where W is designed to ensure that BW contains white noise. Specifically, we requirethatE[(BW)H (BW)]= I, (5.59)where I is the identity matrix. This whitening condition is satisfied by W = R†/2B ,where RB = E[BHB][41, 142] and the superscript †/2 denotes the square root of thepseudo inverse. When W is full-rank, the rank of SW is equal to K. It can be shownthat the same parameters that satisfy the generalised eigenvalue problem in Eq. (5.19)113Chapter 5. Spike inference using Finite Rate of Innovation theoryFigure 5.6: Model order estimation procedure. For a range of prospectivemodel orders, the unknown spike times and amplitudes are estimated and the estima-tion error is computed between the data samples and the resynthesised fluorescencesignal. The selected model order, Ksharp, is the one that causes the sharpest decreasein estimation error and is in the vicinity of Kmin, the order that achieves minimumerror. In A, these values are indicated on the estimation error curve along with thetrue model order, K. In B, we plot the raw data and the resynthesised fluorescencesignals for the estimates with Kmin and Ksharp spikes.for S also do so for SW [142]. The pre-whitened matrix, S̃′, can therefore be used toretrieve the unknown parameters.To calculate W, we note that εs = DGεy, withεs =d0,0 d0,1 ... d0,N−3d1,0 d1,1 ... d1,N−3....... . ....dP,0 dP,1 ... dP,N−3g[0] g[1] g[2] ... 00 g[0] g[1] ... 0.......... . ....0 0 0 ... g[2]εy[0]εy[1]...εy[N − 1] . (5.60)Writing A = DG, the (k, l)th entry of RB is [142][RB]k,l = σ2〈A[M−k:P−k,:], A[M−l:P−l,:],〉F, (5.61)where A[a:b,:] denotes the submatrix of A consisting of rows a to b and 〈·, ·〉F denotesthe Frobenius inner product.1145.4. Spike inference from noisy data5.4.2 Model order estimationTo estimate the model order, K, we adapt the method of Doğan et al. [36] to oursetting. For a range of prospective model orders, J ∈ {1, ..., Jmax}, we estimate theunknown parameters of the spike train, {t̂k, âk}Jk=1. We set Jmax = min(32 λNT,P2),where λ(Hz) is the expected neuronal firing rate. We then compute the estimationerror between the fluorescence signal samples, f [n], and the samples resynthesisedfrom the parameter estimates, f̂J [n] = f̂J(t)|t=nT wherefJ(t) =J∑k=1âkp(t− t̂k). (5.62)The estimation error is computed asET (J) =∥∥f̂J [n]− f [n]∥∥2. (5.63)The estimation error exhibits a U-curve shape, which indicates the transition fromunderfitting (too few spikes) to overfitting (too many), see Fig. 5.6A. The optimalfitting region lies between these two regions. The fitting level is subsequently definedas the minimum estimation error,σ̂2 = minJET (J). (5.64)Prospective estimates of the model order are narrowed down to those in the vicinityof the model order that minimises the estimation error, such thatN(σ̂2) ={J : 12 σ̂2 ≤ ET (J) ≤ 32 σ̂2}. (5.65)Ideally, this region encompasses those in between the underfitting and overfitting re-gions. Finally, the estimate of the model order, K̂, is the order in N(σ̂2) that causesthe largest drop in the estimation error. This is a more conservative estimate sincethe model order that yields the minimum estimation error tends to overfit the data,adding spikes to compensate for noise (Fig. 5.6B).Whereas Doğan et al. compute the estimation error from the sample moments, weuse the raw signal samples. This is because we prefer to perform the model fitting inthe native domain of the spike inference problem, rather than in the FRI domain. Wenote that the amplitude fitting procedure (Section 5.3.2) includes an upper and lower115Chapter 5. Spike inference using Finite Rate of Innovation theoryAlgorithm 1: Procedure by which spike times and amplitudes are estimatedfrom a noisy fluorescence signal.Input: ỹ[n] for n ∈ {0, 1, ..., N − 1}Output: t̂k and âk for k ∈ {1, 2, ..., K̂}1 Apply discrete filter to samples: w̃[n] = ỹ[n] ∗ g[n].2 Compute sample moments: s̃[m] =∑n dm,nw̃[n].3 Compute matrix S̃ from sample moments.4 Perform pre-whitening: S̃′ = S̃W.5 for J ∈ {1, ..., Jmax} do6 Use matrix pencil method to estimate {t̂j}Jj=1 from S̃′.7 Estimate {âj}Jj=1 by least squares from samples and signal model.8 Compute estimation error.9 end10 Select model order, K̂ from {1, ..., Jmax}, by assessing estimation error.11 Retrieve corresponding estimates, {t̂k, âk}K̂k=1, from Steps 6 and 7 with J = K̂.bound on the amplitudes that are obtained from the calcium transient model. Withouta lower bound, the estimation error would be likely to be minimised by overfitting witha large number of spikes with relatively small amplitudes. The amplitude lower boundimplicitly acts as a penalty on the number of spikes, meaning that spikes are nottypically inserted where there is a high degree of noise or fluctuations in the baselinefluorescence.5.4.3 Spike inference algorithmIn Algorithm 1, we summarise the procedure with which spikes are detected fromcalcium imaging data. When we detect spikes from long streams of data, we use thesliding window approach of Oñativia et al. [87], which we describe briefly here. In thisapproach, the fluorescence signal has length N samples but only Nwin < N consecutivesamples are analysed at one time. Firstly, the unknown parameters are estimated inthe stretch of fluorescence signal defined by samples {y[0], y[1], ..., y[Nwin]}. Then, thesliding window is advanced one sample, so that the parameters are estimated for thestretch corresponding to {y[1], y[N ], ..., y[Nwin + 1]}. This process is repeated untilthe sliding window reaches the end of the data. A histogram is then computed of theestimates from all the sliding windows. Typically, estimates that correspond to truespikes are detected consistently in numerous sliding windows. The final spike locationsare then estimated from the peaks of the histogram.1165.5. Experimental methodsFigure 5.7: Example of simulated data. We simulate calcium transients withpulse shape corresponding to the calcium indicator Cal-520 (A). In B, we plot examplesof fluorescence signals produced from the same underlying spike train at SNRs equalto 5 (dB), 10 (dB) and 15 (dB).5.5 Experimental methodsIn Section 5.6, we demonstrate the performance of the spike inference algorithm on realand simulated data. The real data, which was described in detail in Section 4.5.3, wascollected by Dr T. Abrahamsson and Dr P. J. Sjöström. In the following, we presentthe procedures used to generate the simulated data.5.5.1 Simulated dataThe spike trains are generated from a Poisson distribution with varying rate. Wegenerate the fluorescence pulse with parameters corresponding to the commonly-usedsynthetic indicator Cal-520 [133], see Fig. 5.7. As is typical for synthetic indicators,the amplitude of Cal-520 pulses decrease as the local spike rate increases. In thesimulated data, we compute the amplitude of a spike based on the number of spikes inthe previous 250ms time bin, see Table 5.1. This procedure was derived from trendsobserved in real data by Tada et al. [133]. After the spike train has been convolvedwith the pulse to generate the fluorescence signal, white Gaussian noise is added tothe discrete samples, such that f̃ [n] = f [n] + ε[n], where f [n] is the fluorescence signalsampled at rate 16Hz and ε[n] are samples of white Gaussian noise.We vary the SNR of the fluorescence signals in order to fully characterise algorithmperformance. We use the decibel format of SNR, where SNR = 10 log10 Psig/Pnoiseand Psig and Pnoise are the powers of the noiseless fluorescence signal samples, f [n],117Chapter 5. Spike inference using Finite Rate of Innovation theoryNumber of spikes 0 1 2 3 4Amplitude 0.27 0.18 0.18 0.14 0.1Table 5.1: Amplitude of simulated calcium transients. The amplitude ofsimulated spikes varies with the local spike rate. On real data, it has been observedthat amplitudes of calcium transients generated by synthetic indicators decrease asthe local spike rate increases [133]. In our simulations, we compute the amplitude ofa spike based on the number of spikes in the previous 250ms time bin.and the noise samples, ε[n], respectively. We calibrate the standard deviation of thenoise that attains a given SNR on a fluorescence signal containing one spike at t1 = 0observed over the range [0, 1]s (as in Fig. 5.7A). This standard deviation is then usedto attain that SNR in all simulations. This prevents noise power at a given SNR fromvarying with the spike rate. In Fig. 5.7B, we show an example of fluorescence signalsat different noise levels.5.6 Results5.6.1 Comparison with instantaneous-rise algorithmOn simulated data, we examined the temporal error between the true spike timesand spike times inferred by the FRI algorithm. The fluorescence signals, which lasted10s and contained 7 spikes, were generated at three noise levels: 5 dB, 10 dB and15 dB. Fluorescence signals were generated over 1000 realisations of noise at eachSNR. Spikes were detected using the FRI algorithm with both the instantaneous-rise and slow-rise models. The only difference between our implementation of theseversions of the algorithm is that the coefficients of the filter, g[n], and the entries of thepre-whitening matrix, W, differ in the slow-rise and instantaneous-rise case. As thesimulated fluorescence signals stem from a slow-rise pulse model, the instantaneous-riseversion of the algorithm suffers from model mismatch in this setting. In this example,as the model order estimation procedure of both implementations are identical, weassume that K is known. This allows us to restrict the comparison to the temporalerror of each implementation of the algorithm.In Fig. 5.8, we plot the distribution of the temporal error between estimates and truespikes as a fraction of the sampling period, which is T = 1/16. The instantaneous-1185.6. ResultsFigure 5.8: Comparison of the temporal error of the FRI algorithm withinstantaneous-rise and slow-rise models. The fluorescence signals, which aresimulated over 10s and contain 7 spikes, were generated at three noise levels: 5 dB,10 dB and 15 dB. We plot the temporal error between estimates and true spikes (eachestimate is paired with only one true spike) as a fraction of the sample width, T .The distribution of temporal error is plotted over 1000 realisations of noise at eachSNR. The instantaneous-rise model exhibits an average bias of 0.33 samples, due tothe mismatch between the algorithm model and the data.rise algorithm exhibits estimation bias due to mismatch between the algorithm modeland the data. For data with SNR equal to 15 dB, 10 dB and 5 dB, the averageestimation bias is 14s, 21ms and 22ms, respectively. In contrast, using the correct pulsemodel enables average estimation bias of 1ms in all cases. Using the slow-rise model,99.3%, 89.2% and 62.7% of estimates, respectively, have super-resolution precision, i.e.precision above the sample width. While using the correct pulse model seems to haverectified the estimation bias, it does not seem to have affected the variability in theestimates. The standard deviations of the temporal error are 10ms, 20ms, and 40ms,respectively, for both the slow and instantaneous-rise algorithm. We conclude that,using the slow-rise pulse model, which is a more accurate representation of a calciumtransient, enables detection of spikes with higher temporal precision.5.6.2 Model order estimationWe tested the accuracy of the model order estimation procedure on simulated data.Specifically, we simulated fluorescence signals over 8s with varying numbers of spikesand SNRs. We used the procedure presented in Section 5.4.2 to estimate the number ofspikes in the trace. In Fig. 5.9, we display the distribution of the discrepancy betweenthe true number of spikes, K, and the estimated number, K̂. At each noise level and119Chapter 5. Spike inference using Finite Rate of Innovation theoryFigure 5.9: Accuracy of the model order estimation procedure on simulateddata. A fluorescence signal was simulated over 8s with varying numbers of spikes,K, and SNRs. At each combination of K and SNR, 100 realisations of fluorescencesignals were generated. In A, B and C, we display the distribution of the model orderestimation error, K̂ − K, under each set of conditions. In D, E and F, we plot anexample true spike train, x(t), and estimated spike train, x̂(t), at SNR equal to 10(dB) for each value of K.number of true spikes, the procedure is repeated over 100 realisations of noise andspike trains.The estimated model order is selected as the one that causes the sharpest decrease inthe estimation error, whilst being in the vicinity of the model order that minimisesthe estimation error. The estimated value is largely accurate; when SNR is sufficientlyhigh, the estimated model order is predominantly equal to or one less than the truevalue. However, the estimate is generally overly conservative — the number of spikesis more commonly underestimated than overestimated. Furthermore, as the number ofspikes in the time period increases, the accuracy decreases. This can be attributed totwo factors. Firstly, as spike rates increase, the amplitudes of spikes decrease. Theselower amplitude spikes are harder to detect above the noise. Secondly, as the spikerate increases, it is more likely to have closely neighbouring spikes. As is illustrated inFig. 5.9, these are sometimes identified by the algorithm as a single, higher-amplitudespike.1205.6. Results5.6.3 Comparison with state of the artOn imaging data of the mouse visual cortex acquired at a frame rate of 13 Hz, wecompare the performance of the FRI algorithm and a state of the art, non-negativedeconvolution algorithm [44]. Whereas the FRI algorithm uses a continuous timeformulation of the fluorescence signal, the deconvolution algorithm models the signalas an auto-regressive process. At each sample, the fluorescence intensity is assumedto depend upon previous intensity values and the (unknown) number of spikes fired inthat sample. From the data and signal model, a spiking activity vector is inferred, inwhich each entry represents the probability that a spike was fired in a given sample.Spike inference is performed by solving an isotonic regression problem under sparsityconstraints on the spiking activity vector. By subsequently applying a threshold, thespiking probability is converted into the number of spikes per sample.The dataset contains fluorescence signals from 8 cells with ground truth spiking activ-ity, some of which is plotted in Fig. 4.9. We used one training example per cell withwhich to fit the parameters for each algorithm. In the case of the FRI algorithm, wefitted the expected amplitude, A, and rise and decay parameters, τon and τoff. For thedeconvolution algorithm, we identified the most suitable threshold, which is relatedto the calcium transient amplitude, and the optimal lag of the estimates, which isrelated to the rise-time. Then, we applied the algorithms to the remaining 75 tracesin the dataset, see Fig. 5.10. We observed that, on average, the FRI algorithm out-performed the deconvolution algorithm; the algorithms achieved mean CosMIC scoresof 0.53 and 0.47, respectively. We also investigated the source of errors by computingCosMIC’s ancestor metrics, PCosMIC and RCosMIC, which are analogous to the preci-sion and recall. On average, the FRI algorithm was awarded higher scores by both ofthese metrics (Fig. 5.10B). The lower scores of the deconvolution algorithm are likelydue to its quantisation of the temporal interval. As spikes are estimated per sample,the temporal precision of the algorithm is limited by the bin width. Deconvolutionalgorithms may perform relatively better on higher sampling rate or noisier data, onwhich precision above the bin width is unlikely.121Chapter 5. Spike inference using Finite Rate of Innovation theoryFigure 5.10: Algorithm comparison on real data. On a mouse in vitro imagingdataset with ground truth spiking activity, we compare the performance of the FRIalgorithm and a state of the art non-negative deconvolution algorithm. In A, we displaythe distribution of each algorithm’s respective CosMIC scores over all 75 traces. In B,we display the distribution of CosMIC’s ancestor metric scores (PCosMIC and RCosMIC)over all 75 traces.5.7 SummaryWe have presented an extension to the FRI framework for spike inference from calciumimaging data. In the original framework, calcium transients were modelled with aninstantaneous-rise and exponential decay. While this is suitable for calcium indicatorswith fast rises and low sampling rate data, in other cases it can lead to estimationbias due to model mismatch. We have extended the framework to encompass calciumindicators that produce slow-rise pulses of fluorescence. On simulated data, we demon-strated that this increases the temporal precision with which spikes are detected. Tomap the fluorescence signal to the FRI setting, we apply a discrete filter to the samples.This, in turn, alters the statistics of the noise. As the parameter estimation methodswe use are better suited to white noise, in the modified algorithm, we implement aprocedure to whiten the statistics of the noise. On real data from slices of the mousevisual cortex, we demonstrate that the updated algorithm outperforms a state of theart deconvolution algorithm.122Chapter 6.Conclusion6.1 SummaryIn this thesis, we considered the problem of detecting cells and cellular activity fromtwo-photon calcium imaging data. Firstly, we presented an algorithm for detecting thelocations of cells within a video. Next, we considered the problem of spike inference,in which the timing of spikes are inferred from a cell’s fluorescence signal. Our firstcontribution in this area was to develop a spike inference metric, which compares thesimilarity between inferred and ground truth spike trains. Finally, we presented a spikeinference algorithm, whose accuracy we assessed with the aforementioned metric.In Chapter 3, we presented an algorithm, which we refer to as ABLE, for detecting thelocation of cells within an imaging video. In our framework, multiple coupled activecontours evolve, guided by a model-based cost function, to identify cell boundaries.An active contour seeks to partition a local region into two subregions, a cell interiorand exterior, in which all pixels have maximally ‘similar’ time courses. This simple,local model allows contours to be evolved predominantly independently. ABLE isa flexible method: we include no priors on a region’s morphology or stereotypicaltemporal activity. Due to this versatility, ABLE segmented cells with varying size andshape from a mouse in vivo dataset and both active and inactive cells from brain slices.Furthermore, we demonstrated its ability to detect the true boundaries of overlappingcells on both real and simulated data.In Chapter 4, we presented a metric, referred to as CosMIC, for comparing the similar-ity between inferred spikes and ground truth spiking activity. Rather than operating123Chapter 6. Conclusionon the true and estimated spike trains directly, the proposed metric assesses the simi-larity of the pulse trains obtained from convolution of the spike trains with a smoothingpulse. The pulse width is derived from the Cramér-Rao bound (CRB), a lower boundon the variance of any unbiased estimator, which, in turn, is computed from the statis-tics of the imaging data. As such, spike trains receive a score that is implicitly adjustedfor the difficulty of the spike inference problem. The final metric score is the size of thecommonalities of the pulse trains as a fraction of their average size. Viewed throughthe lens of fuzzy set theory, CosMIC resembles a continuous Sørensen-Dice coefficient– an index commonly used to assess the similarity of discrete, presence/absence data.Unlike the spike train correlation, which appears to reward overfitting, the proposedmetric score is maximised when the correct number of spikes have been detected. Fur-thermore, we show that CosMIC is more sensitive to the temporal precision of estimatesthan either the success rate or the spike train correlation.The problem of spike inference can be considered to be one of recovering a finite num-ber of unknown parameters from a noisy signal with known parametric structure. Thisis because neuronal fluorescence signals belong to the class of signals with finite rateof innovation (FRI). In Chapter 5, we extended the FRI framework for spike infer-ence from calcium imaging data. In the original framework, calcium transients weremodelled with an instantaneous-rise and exponential decay. While this is suitable forcalcium indicators with fast rises and low sampling rate data, in other cases it can leadto estimation bias due to model mismatch. We extended the framework to encompasscalcium indicators that produce slow-rise pulses of fluorescence. On simulated data,we demonstrated that this increases the temporal precision with which spikes are de-tected. We also introduced procedures to improve the robustness of the algorithm tonoise. On real data from slices of the mouse visual cortex, we demonstrated that theupdated algorithm outperforms a state of the art deconvolution algorithm.6.2 Future researchWe conclude this chapter with a discussion of potential directions for future research. In our proposed segmentation algorithm, an active contour moves toward a pixelif that pixel’s time course is more similar to the average time course of thecell interior than the narrowband. As techniques for in vivo imaging advance1246.2. Future research(for a review, see [98]), the length of imaging videos, and thus the length ofthe corresponding activity vectors, are consistently increasing. While we havedemonstrated that ABLE’s runtime is virtually unaffected by video length, thisdoes not hold when a video is too large to load into the workstation. It maytherefore become necessary to consider techniques by which the temporal dimen-sion of the video can be reduced prior to segmentation. For example, Pachitariuet al. perform a singular value decomposition (SVD) [94] and Spaen et al. [126]represent each pixel in ‘correlation space’ by a vector of correlation coefficients.Rather than using a generic dimensionality reduction technique such as the SVD,it may be beneficial to develop a bespoke technique, which considers both theknown parametric structure of the data and the manner in which temporal activ-ity vectors are compared in our segmentation algorithm. This could both extendthe range of videos to which ABLE can be applied and, potentially, increase thesegmentation accuracy. An active contour in our proposed segmentation algorithm is advanced using thelevel set method. One of the virtues of the level set method is that it is relativelystraightforward to extend to higher spatial dimensions. For example, Dufour etal. developed a segmentation algorithm for biological imaging data with threespatial dimensions and one temporal dimension [40], which is closely related toChan and Vese’s active contours without edges [24]. With the advent of widefieldmethods such as lightsheet microscopy [1], it is now possible to collect (3+1)-Dcalcium imaging data. Accordingly, future work may focus on the extension ofABLE to higher spatial dimensions. When extending the framework, it maybe useful to consider the impact of photon scattering, which affects widefieldmicroscopy to a greater extent than two-photon microscopy, due to the less-localised excitation of fluorescence. Due to the manner in which calcium ions interact with calcium indicators, theamplitude of calcium transients from one neuron are not constant. Rather, theyvary with the local spike rate of the neuron. At high spike rates, calcium transientamplitudes can change substantially (see, for example [26]), which can negativelyaffect spike detection performance. In future, it could thus be beneficial to incor-porate a model of the variation in spike amplitudes into the FRI spike detectionalgorithm. In the FRI algorithm, locations of spikes are inferred sequentiallyusing a sliding window that advances through the data. It may therefore bebeneficial to use the knowledge of the number of spikes estimated in previous125Chapter 6. Conclusionwindows, to obtain a better estimate of calcium transient amplitudes and, con-sequently, the number and location of spikes.126Appendix A.ABLE: supplemental analytical resultsA.1 External velocityIn Section 3.4.2, we presented the external energy for M ≥ 1 cells. In the following,we compute the external velocity of the mth cell. We recall that Ω̄ represents thelocal spatial domain, which is equal to the union of all contour interiors and theirnarrowbands. Given two sets A ⊂ Ω̄ and B ⊂ Ω̄, we can decompose this asΩ̄ = (A ∩B) ∪ (Ac ∩B) ∪ (A ∩Bc) ∪ (Ac ∩Bc) , (A.1)where Ac denotes the complement of A with respect to Ω̄. In particular, when Arepresents the interior of the mth cell, A = Ωin,m, and B represents the union of theinteriors of all other cells, B =⋃j 6=m Ωin,j , we haveΩ̄ =Ωin,m ∩ ⋃j 6=mΩin,j ∪(Ωin,m)c ∩ ⋃j 6=mΩin,j (A.2)∪Ωin,m ∩ ⋂j 6=m(Ωin,j)c ∪(Ωin,m)c ∩ ⋂j 6=m(Ωin,j)c (A.3)= S1 ∪ S2 ∪ S3 ∪ S4, (A.4)where we use De Morgan’s laws to convert the complement of the union of sets intothe intersection of complements and, in the final line, we introduce a compact notationfor the sets. Writing Cm(x) as the function which reports the cell interiors that arepresent at x ∈ Ω excluding the mth cell, we can decompose the integral in the external127Appendix A. Analytical resultsenergy over the four disjoint sets, such thatEext(φ1, φ2, ..., φM ) =∫outsideD(I(x), fout)dx +∫insideD(I(x),∑m∈C(x)f in,m)dx (A.5)=∫S1D(I(x), f in,m +∑j∈Cm(x)f in,j)dx (A.6)+∫S2D(I(x),∑j∈Cm(x)f in,j)dx (A.7)+∫S3D(I(x), f in,m)dx +∫S4D(I(x), fout)dx (A.8)=∫S1c1(x)dx +∫S2c2(x)dx +∫S3c3(x)dx +∫S4c4(x)dx, (A.9)where, in the final step, we introduce a compact notation for the integrands. Usingthe Heaviside function of the LSFs as an implicit indicator function, we haveEext(φ1, φ2, ..., φM ) =∫Ω̄c1(x)H(φm)∏j∈Cm(x)H(φj)∏k 6∈Cm(x)(1−H(φk)) (A.10)+c2(x)(1−H(φm))∏j∈Cm(x)H(φj)∏k 6∈Cm(x)(1−H(φk)) (A.11)+c3(x)H(φm)∏j 6=m(1−H(φj)) (A.12)+c4(x)(1−H(φm))∏j 6=m(1−H(φj))dx. (A.13)As in Section 3.3.1, to obtain the derivative we substitute the regularised version of theHeaviside function and differentiate with respect to φm. Then, the external velocity is∂Eext∂φm(x) =δε(φm(x)) (c1(x)− c2(x)))∏j∈Cm(x)H(φj)∏k 6∈Cm(x)(1−H(φk)) (A.14)+δε(φm(x)) (c3(x)− c4(x)))∏j 6=m(1−H(φj)). (A.15)Substituting back in the values of ci(x), the external velocity is equal to∂Eext∂φm(x) = δε(φm(x))(D(I(x), f in,m)−D(I(x), fout)), (A.16)128A.2. Euclidean dissimilarity metricwhen x is not in another cell and∂Eext∂φm(x) = δε(φm(x))D(I(x), f in,m + ∑j∈Cm(x)f in,j)−D(I(x),∑j∈Cm(x)f in,j)(A.17)otherwise.A.2 Euclidean dissimilarity metricIn this section, we consider the external velocity of an active contour driven by theEuclidean dissimilarity metric. From Eq. (3.12), we have∂Eext∂φ(x) =δε(φ(x))m(‖I(x)− f in‖2 − ‖I(x)− fout‖2), (A.18)where I(x) is the time course of pixel x, and f in and fout are the average time coursesfrom the contour interior and exterior, respectively. We model the time courses as thesum of a stationary baseline component, the resting fluorescence, and a time-varyingactivity component, which is zero in the absence of activity, such thatfout = bout + aout, f in = bin + ain and I(x) = b + a. (A.19)We have‖I(x)− f in‖2 = ‖b + a−(bin + ain)‖2 (A.20)= ‖(b− bin)+(a− ain)‖2 (A.21)= ‖b− bin‖2 − 2〈b− bin, a− ain〉+ ‖a− ain‖2. (A.22)Similarly, we have‖I(x)− fout‖2 = ‖b− bout‖2 − 2〈b− bout,a− aout〉+ ‖a− aout‖2. (A.23)When cellular activity is correlated with the background and the cell’s baseline fluo-rescence is distinguishable above the background, the elements of a− ain and a− aoutare small compared to b−bin and b−bout. Then, we can approximate the Euclidean129Appendix A. Analytical resultsdistance as‖I(x)− f in‖2 ≈ ‖b− bin‖2 and ‖I(x)− fout‖2 ≈ ‖b− bout‖2. (A.24)In this case, the Euclidean distance in Eq. (A.18) is driven by the discrepancies inbaseline intensity rather than the time-varying activity, such that∂Eext∂φ(x) ≈ δε(φ(x))m(‖b− bin‖2 − ‖b− bout‖2). (A.25)130Appendix B.CosMIC: proof of analytical resultsIn this section we provide proofs of some analytical results presented in Section 4.Proposition B.1 presents a closed-form expression for the score of an estimate of asingle spike. We see that, for estimates within the precision of the pulse width, CosMICmonotonically decreases as estimates are further from the true spike. The result of B.1is then used to establish the pulse width that, on average, awards a score of 0.8 to anestimate at the precision of the CRB. This result and the corresponding derivation arepresented in Proposition B.2. To improve the readability of the formulae, we denotethe metric score of sets of single spikes S = {t1} and Ŝ = {t̂1} as m(t1, t̂1) := M(S, Ŝ).Finally, in Propositions B.3 and B.4, we derive an alternative expression for the metricscore when estimates exactly detect the location of true spikes, but the correct numberof spikes is not detected.Proposition B.1. Let there be a single ground truth spike at t1 and an estimate,t̂1 = t1 + u. The metric score ism (t1, t1 + u) =(|u|2ε − 1)2if |u| < 2ε0 otherwise,(B.1)where ε is half the width of the pulse.Proof. Without loss of generality, we let the true spike location be at t1 = 0, as themetric score depends on the relative rather than absolute locations of the estimated131Appendix B. Analytical resultsand ground truth spikes. From Eq. (4.25), we havem(0, u) = 2‖min(qε(t), qε(t− u))‖1‖qε(t)‖1 + ‖qε(t− u)‖1. (B.2)When |u| > 2ε, the pulses do not overlap and, consequently, the numerator is equal to0. Therefore, the metric score is zero for all |u| > 2ε. For |u| ≤ 2ε, we writem(0, u) =1ε(∫Aqε(t) dt+∫Bqε(t− u) dt), (B.3)which follows from ‖qε‖ = ε, A = {t ∈ R : qε(t) < qε(t− u)} and B = {t ∈ R : qε(t) ≥qε(t − u)}. From the change of variables v = t + u, we see that m(0, u) = m(0,−u).As M is even in the second argument, we must only evaluate m(0, u) for 0 < u < 2ε.To identify the support of A and B, we must identify the value of t at which qε(t) =qε(t− u). We haveqε(t) = qε(t− u)⇔ 1−|t|ε= 1− |t− u|ε⇔ |t| = |t− u|. (B.4)For 0 < u < 2ε, the relevant intersection point is at t = u/2. Omitting the range ofintegration at which the integrand is zero, Eq. (B.3) becomesm(0, u) =1ε(∫ εu/2qε(t) dt+∫ u/2u−εqε(t− u) dt)(B.5)=1ε(∫ εu/2qε(t) dt+∫ −u/2−εqε(v) dv)(B.6)=2ε∫ εu/2qε(t) dt, (B.7)which follows from the change of variables v = t+ u and the symmetry of qε(t) about0. Evaluating the integral, we obtain m(0, u) = (|u|/2ε− 1)2 for |u| < 2ε.Proposition B.2. Let t1 denote the location of the true spike. The spike estimate,modelled with the random variable U ∼ N(t1, σ2CRB), is normally distributed aboutthe true spike at the precision of the CRB. We denote β = σCRB/w, where w = 2ε isthe pulse width. Then, the expectation of the metric score is 0.8 if β satisfies0.4 = (Φ(1/β)− 0.5)(β2 + 1)+β√2π(exp(−1/2β2)− 2), (B.8)132where Φ denotes the cumulative distribution function of the standard normal distribu-tion.Proof. Without loss of generality, we consider the case where t1 = 0. Denoting fU (u)as the PDF of U , we haveEU [m(0, U)] =∫Rm(0, u)fU (u)du (B.9)= 2∫ w0( uw− 1)2fU (u)du (B.10)=2w2∫ w0u2fU (u)du−4w∫ w0ufU (u)du+ 2∫ w0fU (u)du (B.11)=2w2I1 −4wI2 + 2I3, (B.12)where Eq. (B.10) follows from Proposition B.1. Applying integration by parts to I1,we obtainI1 = σ2CRBP (U ∈ [0, w))− σ2CRBfU (w)w. (B.13)The remaining integrals are I2 = −σ2CRB(fU (w) − fU (0)) and I3 = P (U ∈ [0, w)),respectively. Putting the integrals together yieldsE [m(0, U)] = 2[P (U ∈ [0, w))(σ2CRBw2+ 1)+σ2CRBw(fU (w)− 2fU (0))]. (B.14)Writing β = σCRB/w and expanding the distribution function of U , we haveEU [m(0, U)] = 2((Φ(1/β)− 0.5)(β2 + 1)+β√2π(exp(−1/2β2)− 2)). (B.15)Proposition B.3. Let S be a set of K true spikes and Ŝ a set of K̂ estimates. IfŜ contains a subset of the ground truth spike times with the exception of R missingspikes and no extras, such that K̂ = K −R with 0 < R ≤ K, thenM(S, Ŝ) = 1− 12K/R− 1. (B.16)Proof. We denote with xr(t) and r(t) the spike train and pulse train, respectively, ofthe spikes present in S but missing from Ŝ. Due to the distributivity of the convolution133Appendix B. Analytical resultsoperation, we haveµ̂(t) = x̂(t) ∗ qε(t) = (x(t)− xr(t)) ∗ qε(t) = µ(t)− r(t).From the metric expression in Eq. (4.32), the metric score becomesM(S, Ŝ) = 1− ‖µ− µ̂‖1‖µ‖1 + ‖µ̂‖1(B.17)= 1− ‖µ− (µ− r)‖1‖µ‖1 + ‖µ− r‖1(B.18)= 1− R‖qε‖1K‖qε‖1 + (K −R) ‖qε‖1(B.19)= 1− 12K/R− 1. (B.20)Proposition B.4. Let S be set of K spikes and Ŝ a set of K̂ estimates containing allthe ground truth spikes together with R > 0 extras, such that K̂ = K +R. ThenM(S, Ŝ) =11 +R/2K. (B.21)Proof. We denote with xr(t) and r(t) the spike train and pulse train, respectively, ofthe spikes present in Ŝ but missing from S. Due to the distributivity of the convolutionoperator, the estimated pulse train can be writtenµ̂(t) = x̂(t) ∗ qε(t) = (x(t) + xr(t)) ∗ qε(t) = µ(t) + r(t).The metric score isM(S, Ŝ) =2‖min(µ, µ̂)‖1‖µ‖1 + ‖µ̂‖1(B.22)=2‖min(µ, µ+ r)‖1‖µ‖1 + ‖µ+ r‖1(B.23)=2‖µ‖12‖µ‖1 + ‖r‖1(B.24)=11 + ‖r‖1/2‖µ‖1, (B.25)where the penultimate line follows from the non-negativity of µ and r. As ‖µ‖1 =134K‖qε‖1 and ‖r‖1 = R‖qε‖1, it follows thatM(S, Ŝ) =11 +R/2K. (B.26)135Appendix B. Analytical results136Bibliography[1] M. B. Ahrens, M. B. Orger, D. N. Robson, J. M. Li, and P. J. Keller. Whole-brainfunctional imaging at cellular resolution using light-sheet microscopy. Naturemethods, 10:413–420, 2013.[2] J. Akerboom, T.-W. Chen, T. J. Wardill, L. Tian, J. S. Marvin, S. Mutlu, N. C.Calderón, F. Esposti, B. G. Borghuis, X. R. Sun, A. Gordus, M. B. Orger,R. Portugues, F. Engert, J. J. Macklin, A. Filosa, A. Aggarwal, R. A. Kerr,R. Takagi, S. Kracun, E. Shigetomi, B. S. Khakh, H. Baier, L. Lagnado, S. S.-H.Wang, C. I. Bargmann, B. E. Kimmel, V. Jayaraman, K. Svoboda, D. S. Kim,E. R. Schreiter, and L. L. Looger. Optimization of a GCaMP calcium indicatorfor neural activity imaging. Journal of Neuroscience, 32(40):13819–13840, 2012.[3] D. G. Altman and J. M. Bland. Measurement in medicine: the analysis ofmethod comparison studies. Journal of the Royal Statistical Society. Series D(The Statistician), 32(3):307–317, 1983.[4] N. J. Apthorpe, A. J. Riordan, R. E. Aguilar, J. Homann, Y. Gu, D. W. Tank,and H. S. Seung. Automatic neuron detection in calcium imaging data usingconvolutional networks. In Proceedings of the 29th International Conference onNeural Information Processing Systems, NIPS’16, pages 3278–3286, USA, 2016.Curran Associates Inc.[5] G. Aubert and P. Kornprobst. Mathematical Problems in Image Processing:Partial Differential Equations and the Calculus of Variations. Springer-VerlagNew York, second edition, 2006.[6] R. Bellman and M. Giertz. On the analytic formalism of the theory of fuzzy sets.Information Sciences, 5:149–156, 1973.[7] P. Berens, J. Freeman, T. Deneux, N. Chenkov, T. McColgan, A. Speiser, J. H.Macke, S. C. Turaga, P. Mineault, P. Rupprecht, S. Gerhard, R. W. Friedrich,137BibliographyJ. Friedrich, L. Paninski, M. Pachitariu, K. D. Harris, B. Bolte, T. A. Machado,D. Ringach, J. Stone, L. E. Rogerson, N. J. Sofroniew, J. Reimer, E. Froudarakis,T. Euler, M. Romn Rosn, L. Theis, A. S. Tolias, and M. Bethge. Community-based benchmarking improves spike rate inference from two-photon calciumimaging data. PLOS Computational Biology, 14(5):1–13, 2018.[8] J. Berent and P. L. Dragotti. Plenoptic manifolds. IEEE Signal ProcessingMagazine, 24(6):34–44, 2007.[9] M. J. Berridge, P. Lipp, and M. D. Bootman. The versatility and universality ofcalcium signalling. Nature Reviews Molecular Cell Biology, 1:11 – 21, 2000.[10] V. Bonin, M. H. Histed, S. Yurgenson, and R. C. Reid. Local diversity andfine-scale organization of receptive fields in mouse visual cortex. Journal ofNeuroscience, 31(50):18506–18521, 2011.[11] J. G. G. Borst and B. Sakmann. Calcium current during a single action potentialin a large presynaptic terminal of the rat brainstem. The Journal of Physiology,506(Pt 1):143–157, 1998.[12] M. B. Bouchard, V. Voleti, C. S. Mendes, C. Lacefield, W. B. Grueber, R. S.Mann, R. M. Bruno, and M. C. Hillman. Swept confocally-aligned planar ex-citation (SCAPE) microscopy for high-speed volumetric imaging of behavingorganisms. Nature Photonics, 9:113–119, 2015.[13] J. R. Bray and J. T. Curtis. An ordination of the upland forest communities ofsouthern Wisconsin. Ecological Monographs, 27(4):325–349, 1957.[14] R. Brette and A. Destexhe. Handbook of Neural Activity Measurement. Cam-bridge University Press, 2012.[15] P. Brigger, J. Hoeg, and M. Unser. B-spline snakes: a flexible tool for parametriccontour detection. IEEE Transactions on Image Processing, 9(9):1484–1496,2000.[16] K. A. Buchanan, A. V. Blackman, A. W. Moreau, D. Elgar, R. P. Costa,T. Lalanne, A. A. Tudor Jones, J. Oyrer, and P. J. Sjöström. Target-specificexpression of presynaptic NMDA receptors in neocortical microcircuits. Neuron,75(3):451–466, 2012.138Bibliography[17] E. Bullmore and O. Sporns. Complex brain networks: graph theoretical analysisof structural and functional systems. Nature Reviews Neuroscience, 10:186–198,2009.[18] G. Buzsáki. Large-scale recording of neuronal ensembles. Nature Neuroscience,7:446–451, 2004.[19] J. A. Cadzow. Signal enhancement-a composite property mapping algorithm.IEEE Transactions on Acoustics, Speech, and Signal Processing, 36(1):49–62,1988.[20] V. Caselles, F. Catté, T. Coll, and F. Dibos. A geometric model for activecontours in image processing. Numerische Mathematik, 66(1):1–31, 1993.[21] V. Caselles, R. Kimmel, and G. Sapiro. Geodesic active contours. InternationalJournal of Computer Vision, 22(1):61–79, 1997.[22] T. Chan and L. Vese. An active contour model without edges. In M. Nielsen,P. Johansen, O. F. Olsen, and J. Weickert, editors, Scale-Space Theories inComputer Vision, pages 141–151, Berlin, Heidelberg, 1999. Springer Berlin Hei-delberg.[23] T. F. Chan, B. Sandberg, and L. A. Vese. Active contours without edges forvector-valued images. Journal of Visual Communication and Image Representa-tion, 11(2):130–141, 2000.[24] T. F. Chan and L. A. Vese. Active contours without edges. IEEE Transactionson Image Processing, 10(2):266–277, 2001.[25] T. F. Chan and L. A. Vese. A level set algorithm for minimizing the Mumford-Shah functional in image processing. In Proceedings IEEE Workshop on Vari-ational and Level Set Methods in Computer Vision, pages 161–168, Vancouver,BC, 2001. IEEE.[26] T.-W. Chen, T. J. Wardill, Y. Sun, S. R. Pulver, S. L. Renninger, A. Baohan,E. R. Schreiter, R. A. Kerr, M. B. Orger, V. Jayaraman, L. L. Looger, K. Svo-boda, and D. S. Kim. Ultrasensitive fluorescent proteins for imaging neuronalactivity. Nature, 499:295–300, 2013.[27] A. Danielyan, Y.-W. Wu, P.-Y. Shih, Y. Dembitskaya, and A. Semyanov. Denois-ing of two-photon fluorescence images with block-matching 3D filtering. Methods,68(2):308–316, 2014.139Bibliography[28] R. Delgado-Gonzalo, V. Uhlmann, D. Schmitter, and M. Unser. Snakes on aplane: A perfect snap for bioimage analysis. IEEE Signal Processing Magazine,32(1):41–48, 2015.[29] R. Delgado-Gonzalo and M. Unser. Spline-based framework for interactive seg-mentation in biomedical imaging. IRBM, 34(3):235–243, 2013.[30] T. Deneux, A. Kaszas, G. Szalay, K. Katona, T. Lakner, A. Grinvald, B. Rózsa,and I. Vanzetta. Accurate spike estimation from noisy calcium signals for ul-trafast three-dimensional imaging of large neuronal populations in vivo. NatureCommunications, 7:12190, 2016.[31] W. Denk, J. H. Strickler, and W. W. Webb. Two-photon laser scanning fluores-cence microscopy. Science, 248(4951):73–76, 1990.[32] L. R. Dice. Measures of the amount of ecologic association between species.Ecology, 26(3):297–302, 1945.[33] A. S. Dickey, A. Suminski, Y. Amit, and N. G. Hatsopoulos. Single-unit stabilityusing chronically implanted multielectrode arrays. Journal of Neurophysiology,102(2):1331–1339, 2009.[34] F. Diego Andilla and F. A. Hamprecht. Sparse space-time deconvolution forcalcium image analysis. In Proceedings of the 27th International Conferenceon Neural Information Processing Systems, NIPS’14, pages 64–72, USA, 2014.Curran Associates Inc.[35] Z. Doǧan, T. Blu, and D. Van De Ville. Detecting spontaneous brain activity infunctional magnetic resonance imaging using finite rate of innovation. In 11thInternational Symposium on Biomedical Imaging, ISBI, pages 1047–1050. IEEE,2014.[36] Z. Doǧan, C. Gilliam, T. Blu, and D. Van De Ville. Reconstruction of finite rateof innovation signals with model-fitting approach. IEEE Transactions on SignalProcessing, 63(22):6024–6036, 2015.[37] R. J. Douglas and K. A. C. Martin. Neuronal circuits of the neocortex. AnnualReview of Neuroscience, 27(1):419–451, 2004.[38] P. L. Dragotti, M. Vetterli, and T. Blu. Sampling moments and reconstructingsignals of finite rate of innovation: Shannon meets Strang-Fix. IEEE Transac-tions on Signal Processing, 55(5):1741–1757, 2007.140Bibliography[39] D. Dubois and H. Prade. Fuzzy sets and systems: theory and applications. Aca-demic Press, New York, 1980.[40] A. Dufour, V. Shinin, S. Tajbakhsh, N. Guillen-Aghion, J. C. Olivo-Marin, andC. Zimmer. Segmenting and tracking fluorescent cells in dynamic 3-D microscopywith coupled active surfaces. IEEE Transactions on Image Processing, 14(9):1396–1410, 2005.[41] Y. C. Eldar and A. V. Oppenheim. MMSE whitening and subspace whitening.IEEE Transactions on Information Theory, 49(7):1846–1851, 2003.[42] M. A. T. Figueiredo, J. M. N. Leitao, and A. K. Jain. Unsupervised contourrepresentation and estimation using B-splines and a minimum description lengthcriterion. IEEE Transactions on Image Processing, 9(6):1075–1087, 2000.[43] A. Foi, M. Trimeche, V. Katkovnik, and K. Egiazarian. Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data. IEEE Transac-tions on Image Processing, 17(10):1737–1754, 2008.[44] J. Friedrich, P. Zhou, and L. Paninski. Fast online deconvolution of calciumimaging data. PLOS Computational Biology, 13(3):1–26, 2017.[45] A. Giovannucci, J. Friedrich, M. Kaufman, A. Churchland, D. Chklovskii,L. Paninski, and E. A. Pnevmatikakis. OnACID: Online analysis of calciumimaging data in real time. In Proceedings of the 30th International Conferenceon Neural Information Processing Systems, NIPS’17, pages 2378–2388, USA,2017. Curran Associates Inc.[46] W. Göbel and F. Helmchen. In vivo calcium imaging of neural network function.Physiology, 22(6):358–365, 2007.[47] B. F. Grewe, D. Langer, H. Kasper, B. M. Kampa, and F. Helmchen. High-speedin vivo calcium imaging reveals neuronal network activity with near-millisecondprecision. Nature Methods, 7:399–405, 2010.[48] C. Grienberger and A. Konnerth. Imaging calcium in neurons. Neuron, 73(5):862–885, 2012.[49] Z. V. Guo, N. Li, D. Huber, E. Ophir, D. Gutnisky, J. T. Ting, G. Feng, andK. Svoboda. Flow of cortical activity underlying a tactile decision in mice.Neuron, 81(1):179–194, 2014.141Bibliography[50] O. P. Hamill, A. Marty, E. Neher, B. Sakmann, and F. J. Sigworth. Improvedpatch-clamp techniques for high-resolution current recording from cells and cell-free membrane patches. Pflügers Archiv, 391(2):85–100, 1981.[51] K. D. Harris, R. Q. Quiroga, J. Freeman, and S. Smith. Improving data qualityin neuronal population recordings. Nature Neuroscience, 19:1165–1174, 2016.[52] F. Helmchen and W. Denk. Deep tissue two-photon microscopy. Nature Methods,2:932–940, 2005.[53] Y. Hua and T. K. Sarkar. Matrix pencil method for estimating parametersof exponentially damped/undamped sinusoids in noise. IEEE Transactions onAcoustics, Speech, and Signal Processing, 38(5):814–824, 1990.[54] Y. Hua and T. K. Sarkar. On SVD for estimating generalized eigenvalues ofsingular matrix pencil in noise. In IEEE International Sympoisum on Circuitsand Systems, volume 5, pages 2780–2783, 1991.[55] D. Huber, D. A. Gutnisky, S. Peron, D. H. O’Connor, J. S. Wiegert, L. Tian,T. G. Oertner, L. L. Looger, and K. Svoboda. Multiple dynamic representationsin the motor cortex during sensorimotor learning. Nature, 484:473–478, 2012.[56] H. Inan, M. A. Erdogdu, and M. Schnitzer. Robust estimation of neural signals incalcium imaging. In Proceedings of the 30th International Conference on NeuralInformation Processing Systems, NIPS’17, pages 2901–2910, USA, 2017. CurranAssociates Inc.[57] M. Jacob, T. Blu, and M. Unser. Efficient energies and algorithms for parametricsnakes. IEEE Transactions on Image Processing, 13(9):1231–1244, 2004.[58] S. Jehan-Besson, M. Barlaud, and G. Aubert. DREAM2S: Deformable regionsdriven by an Eulerian accurate minimization method for image and video seg-mentation. International Journal of Computer Vision, 53(1):45–70, 2003.[59] N. Ji, J. Freeman, and S. L. Smith. Technologies for imaging neural activity inlarge volumes. Nature Neuroscience, 19:1154–1164, 2016.[60] P. Kaifosh, J. D. Zaremba, N. B. Danielson, and A. Losonczy. SIMA: Pythonsoftware for analysis of dynamic fluorescence imaging data. Frontiers in Neu-roinformatics, 8:80, 2014.142Bibliography[61] M. Kass, A. Witkin, and D. Terzopoulos. Snakes: Active contour models. Inter-national Journal of Computer Vision, 1(4):321–331, 1988.[62] S. M. Kay. Fundamentals of Statistical Signal Processing: Estimation Theory.Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1993.[63] A. M. Kerlin, M. L. Andermann, V. K. Berezovskii, and R. C. Reid. Broadlytuned response properties of diverse inhibitory neuron subtypes in mouse visualcortex. Neuron, 67(5):858–871, 2010.[64] J. N. D. Kerr, D. Greenberg, and F. Helmchen. Imaging input and output ofneocortical networks in vivo. Proceedings of the National Academy of Sciences,102(39):14063–14068, 2005.[65] A. Klibisz, D. Rose, M. Eicholtz, J. Blundon, and S. Zakharenko. Fast, simplecalcium imaging segmentation with fully convolutional networks. In M. J. Car-doso et al., editors, Deep Learning in Medical Image Analysis and MultimodalLearning for Clinical Decision Support, DLMIA 2017, ML-CDS 2017. LectureNotes in Computer Science, vol 10553. Springer, Cham, 2017.[66] J. J. Knierim. Chapter 19 - Information processing in neural networks. InJ. H. Byrne, R. Heidelberger, and M. N. Waxham, editors, From Molecules toNetworks, pages 563 – 589. Academic Press, Boston, 2014.[67] V. A. Kotelnikov. On the transmission capacity of the ether and wire in electro-communications. Izd. Red. Upr. Svyazzi RKKA (Moscow), 1933.[68] C. Li, R. Huang, Z. Ding, J. C. Gatenby, D. N. Metaxas, and J. C. Gore. A levelset method for image segmentation in the presence of intensity inhomogeneitieswith application to MRI. IEEE Transactions on Image Processing, 20(7):2007–2016, 2011.[69] C. Li, C. Xu, C. Gui, and M. D. Fox. Distance regularized level set evolution andits application to image segmentation. IEEE Transactions on Image Processing,19(12):3243–3254, 2010.[70] A. Longstaff. Neuroscience. BIOS instant notes series. Taylor and Francis,London, third edition, 2011.[71] L. L. Looger and O. Griesbeck. Genetically encoded neural activity indicators.Current Opinion in Neurobiology, 22(1):18–23, 2012.143Bibliography[72] Z. Lu, G. Carneiro, and A. P. Bradley. Automated nucleus and cytoplasmsegmentation of overlapping cervical cells. In Medical Image Computing andComputer-Assisted Intervention, MICCAI 2013, pages 452–460, Berlin, 2013.Springer.[73] H. Lütcke, F. Gerhard, F. Zenke, W. Gerstner, and F. Helmchen. Inferenceof neuronal network spike dynamics and topology from calcium imaging data.Frontiers in Neural Circuits, 7:201, 2013.[74] R. Malladi, J. A. Sethian, and B. C. Vemuri. Shape modeling with front propaga-tion: a level set approach. IEEE Transactions on Pattern Analysis and MachineIntelligence, 17(2):158–175, 1995.[75] S. Mallat. A wavelet tour of signal processing: the sparse way. Academic press,third edition, 2009.[76] R. Maruyama, K. Maeda, H. Moroda, I. Kato, M. Inoue, H. Miyakawa, andT. Aonishi. Detecting cells using non-negative matrix factorization on calciumimaging data. Neural Networks, 55:11 –19, 2014.[77] R. McAulay and E. Hofstetter. Barankin bounds on parameter estimation. IEEETransactions on Information Theory, 17(6):669–676, 1971.[78] A. Miri, K. Daie, R. D. Burdine, E. Aksay, and D. W. Tank. Regression-basedidentification of behavior-encoding neurons during large-scale optical imaging ofneural activity at cellular resolution. Journal of Neurophysiology, 105(2):964–980, 2011.[79] A. Mitchie, R. Feghali, and A. Mansouri. Motion tracking as a spatio-temporalmotion boundary detection. Robotics and autonomous systems, 43:39–50, 2003.[80] A. Miyawaki, J. Llopis, R. Heim, J. M. McCaffery, J. A. Adams, M. Ikura, andR. Y. Tsien. Fluorescent indicators for Ca2+ based on green fluorescent proteinsand calmodulin. Nature, 388:882–887, 1997.[81] V. B. Mountcastle. An organizing principle for cerebral function: The unit modeland the distributed system. In G. M. Edelman and V. V. Mountcastle, editors,The Mindful Brain, pages 7–50. MIT Press, Cambridge, MA, 1978.[82] E. A. Mukamel, A. Nimmerjahn, and M. J. Schnitzer. Automated analysis ofcellular signals from large-scale calcium imaging data. Neuron, 63(6):747–760,2009.144Bibliography[83] D. Mumford and J. Shah. Optimal approximations by piecewise smooth func-tions and associated variational problems. Communications on Pure and AppliedMathematics, 42(5):577–685, 1989.[84] J. Murray-Bruce and P. L. Dragotti. Estimating localized sources of diffusionfields using spatiotemporal sensor measurements. IEEE Transactions on SignalProcessing, 63(12):3018–3031, 2015.[85] http://neurofinder.codeneuro.org/. Accessed: 2018-01-04.[86] K. Ohki, S. Chung, P. Kara, M. H ubener, T. Bonhoeffer, and R. C. Reid. Highlyordered arrangement of single neurons in orientation pinwheels. Nature, 442:925–928, 2006.[87] J. Oñativia, S. R. Schultz, and P. L. Dragotti. A finite rate of innovation algo-rithm for fast and accurate spike detection from two-photon calcium imaging.Journal of Neural Engineering, 10(4):046017, 2013.[88] J. Oñativia Bravo. Sampling and reconstruction of finite rate of innovationssignals with applications in neuroscience and sparse representation. PhD thesis,Imperial College London, 2015.[89] S. Osher and R. Fedkiw. Level set methods and dynamic implicit surfaces, volume153 of Applied Mathematical Sciences. Springer-Verlag New York, 2003.[90] S. Osher and N. Paragios. Geometric level set methods in imaging, vision, andgraphics. Springer-Verlag New York, 2003.[91] S. Osher and J. A. Sethian. Fronts propagating with curvature-dependent speed:algorithms based on Hamilton-Jacobi formulations. Journal of computationalphysics, 79(1):12–49, 1988.[92] M. Pachitariu, A. M. Packer, N. Pettit, H. Dalgleish, M. Häusser, and M. Sahani.Extracting regions of interest from biological images with convolutional sparseblock coding. In Proceedings of the 26th International Conference on NeuralInformation Processing Systems, NIPS’13, pages 1745–1753, USA, 2013. CurranAssociates Inc.[93] M. Pachitariu, C. Stringer, and K. D. Harris. Robustness of spike deconvolutionfor calcium imaging of neural spiking. bioRxiv, 2017. doi:10.1101/156786.145Bibliography[94] M. Pachitariu, C. Stringer, S. Schröder, M. Dipoppa, L. F. Rossi, M. Carandini,and K. D. Harris. Suite2p: beyond 10,000 neurons with standard two-photonmicroscopy. bioRxiv, 2016. doi:10.1101/061507.[95] A. R. C. Paiva, I. Park, and J. C. Pŕıncipe. A comparison of binless spike trainmeasures. Neural Computing and Applications, 19(3):405–419, 2010.[96] C. P. Pappis and N. I. Karacapilidis. A comparative assessment of measures ofsimilarity of fuzzy values. Fuzzy Sets and Systems, 56(2):171–174, 1993.[97] J. B. Pawley. Handbook of Biological Confocal Microscopy. Springer US, thirdedition, 2006.[98] S. Peron, T.-W. Chen, and K. Svoboda. Comprehensive imaging of corticalnetworks. Current Opinion in Neurobiology, 32:115–123, 2015.[99] S. Peron, J. Freeman, V. Iyer, C. Guo, and K. Svoboda. Volumetric calciumimaging data recorded during performance of a single whisker object localizationtask, sampling activity in the majority of the relevant superficial barrel cortexneurons (75%, 12,000 neurons per mouse). CRCNS.org. doi:10.6080/K0TB14TN,2015.[100] S. P. Peron, J. Freeman, V. Iyer, C. Guo, and K. Svoboda. A cellular resolutionmap of barrel cortex activity during tactile behavior. Neuron, 86(3):783–799,2015.[101] D. S. Peterka, H. Takahashi, and R. Yuste. Imaging voltage in neurons. Neuron,69(1):9–21, 2011.[102] E. C. Pielou. The interpretation of ecological data: A primer on classificationand ordination. John Wiley and Sons, Inc., 1984.[103] E. A. Pnevmatikakis, J. Merel, A. Pakman, and L. Paninski. Bayesian spikeinference from calcium imaging data. In 2013 Asilomar Conference on Signals,Systems and Computers, pages 349–353. IEEE, 2013.[104] E. A. Pnevmatikakis and L. Paninski. Sparse nonnegative deconvolution forcompressive calcium imaging: algorithms and phase transitions. In Proceedingsof the 26th International Conference on Neural Information Processing Systems,NIPS’13, pages 1250–1258, USA, 2013. Curran Associates Inc.146Bibliography[105] E. A. Pnevmatikakis, D. Soudry, Y. Gao, T. A. Machado, J. Merel, D. Pfau,T. Reardon, Y. Mu, C. Lacefield, W. Yang, et al. Simultaneous denoising, de-convolution, and demixing of calcium imaging data. Neuron, 89(2):285–299,2016.[106] V. Rahmati, K. Kirmse, D. Markovi, K. Holthoff, and S. J. Kiebel. Inferringneuronal dynamics from calcium imaging data using biophysical models andbayesian inference. PLOS Computational Biology, 12(2):1–42, 2016.[107] S. Reynolds, T. Abrahamsson, R. Schuck, P. J. Sjöström, S. R. Schultz, andP. L. Dragotti. ABLE: An activity-based level set segmentation algorithm fortwo-photon calcium imaging data. eNeuro, 4(5), 2017.[108] S. Reynolds, T. Abrahamsson, P. J. Sjöström, S. R. Schultz, and P. L. Dragotti.CosMIC: A consistent metric for spike inference from calcium imaging. NeuralComputation, In press.[109] S. Reynolds, C. S. Copeland, S. R. Schultz, and P. L. Dragotti. An extensionof the FRI framework for calcium transient detection. In 13th InternationalSymposium on Biomedical Imaging, ISBI, pages 676–679. IEEE, 2016.[110] M. Ristivojevic and J. Konrad. Space-time image sequence analysis: objecttunnels and occlusion volumes. IEEE Transactions on Image Processing, 15(2):364–376, 2006.[111] D. W. Roberts. Ordination on the basis of fuzzy set theory. Vegetatio, 66(3):123–131, 1986.[112] N. L. Rochefort, H. Jia, and A. Konnerth. Calcium imaging in the living brain:prospects for molecular medicine. Trends in Molecular Medicine, 14(9):389–399,2008.[113] H. L. Royden. Real analysis. Prentice Hall, Englewood Cliffs, third edition, 1988.[114] T. Sasaki, N. Takahashi, N. Matsuki, and Y. Ikegaya. Fast and accurate detectionof action potentials from somatic calcium fluctuations. Journal of Neurophysiol-ogy, 100(3):1668–1676, 2008.[115] M. Scanziani and M. Häusser. Electrophysiology in the age of light. Nature, 461:930–939, 2009.147Bibliography[116] S. Schreiber, J. Fellous, D. Whitmer, P. Tiesinga, and T. Sejnowski. A newcorrelation-based measure of spike timing reliability. Neurocomputing, 52-54(Supplement C):925–931, 2003.[117] R. Schuck, M. A. Go, S. Garasto, S. Reynolds, P. L. Dragotti, and S. R. Schultz.Multiphoton minimal inertia scanning for fast acquisition of neural activity sig-nals. Journal of Neural Engineering, 15(2):025003, 2018.[118] S. R. Schultz, K. Kitamura, A. Post-Uiterweer, J. Krupic, and M. Häusser.Spatial pattern coding of sensory information by climbing fiber-evoked calciumsignals in networks of neighboring cerebellar purkinje cells. Journal of Neuro-science, 29(25):8005–8015, 2009.[119] C. E. Shannon. Communication in the presence of noise. Proceedings of theIEEE, 37(1):10–21, 1949.[120] M. E. J. Sheffield and D. A. Dombeck. Calcium transient prevalence across thedendritic arbour predicts place field properties. Nature, 517:200–204, 2015.[121] D. Slepian. On bandwidth. Proceedings of the IEEE, 64(3):292–300, 1976.[122] S. L. Smith and M. Häusser. Parallel processing of visual space by neighboringneurons in mouse visual cortex. Nature Neuroscience, 13:1144–1149, 2010.[123] N. J. Sofroniew, J. D. Cohen, A. K. Lee, and K. Svoboda. Natural whisker-guidedbehavior by head-fixed mice in tactile virtual reality. Journal of Neuroscience,34(29):9537–9550, 2014.[124] N. J. Sofroniew, D. Flickinger, J. King, and K. Svoboda. A large field of viewtwo-photon mesoscope with subcellular resolution for in vivo imaging. eLife, 5:e14472, 2016.[125] T. Sørensen. A method of establishing groups of equal amplitude in plant sociol-ogy based on similarity of species and its application to analyses of the vegetationon danish commons. Biologiske Skrifter, 5:1–34, 1948.[126] Q. Spaen, D. S. Hochbaum, and R. Aśın-Achá. HNCcorr: A Novel CombinatorialApproach for Cell Identification in Calcium-Imaging Movies. ArXiv e-prints,2017. arXiv:1703.01999.[127] A. Speiser, J. Yan, E. W. Archer, L. Buesing, S. C. Turaga, and J. H. Macke. Fastamortized inference of neural activity from calcium imaging data with variational148Bibliographyautoencoders. In Proceedings of the 30th International Conference on NeuralInformation Processing Systems, NIPS’17, pages 4024–4034, USA, 2017. CurranAssociates Inc.[128] J. M. Squirrell, D. L. Wokosin, J. G. White, and B. D. Bavister. Long-termtwo-photon fluorescence imaging of mammalian embryos without compromisingviability. Nature Biotechnology, 17:763–767, 1999.[129] P. Stoica and R. Moses. Spectral Analysis of Signals. Prentice Hall, 2005.[130] C. Stosiek, O. Garaschuk, K. Holthoff, and A. Konnerth. In vivo two-photoncalcium imaging of neuronal networks. Proceedings of the National Academy ofSciences, 100(12):7319–7324, 2003.[131] K. Svoboda and S. M. Block. Biological applications of optical forces. AnnualReview of Biophysics and Biomolecular Structure, 23:247–285, 1994.[132] K. Svoboda (contact). GENIE project, Janelia Farm Campus, HHMI. (2015).Simultaneous imaging and loose-seal cell-attached electrical recordings from neu-rons expressing a variety of genetically encoded calcium indicators. CRCNS.org.doi:10.6080/k02r3pmn.[133] M. Tada, A. Takeuchi, M. Hashizume, K. Kitamura, and M. Kano. A highlysensitive fluorescent indicator dye for calcium imaging of neural activity in vitroand in vivo. European Journal of Neuroscience, 39(11):1720–1728, 2014.[134] A. E. Taylor and W. R. Mann. Advanced calculus. John Wiley and Sons, Inc.,USA, third edition, 1983.[135] L. Theis, P. Berens, E. Froudarakis, J. Reimer, M. R. Rosón, T. Baden, T. Euler,A. S. Tolias, and M. Bethge. Benchmarking spike rate inference in populationcalcium imaging. Neuron, 90(3):471–482, 2016.[136] L. Tian, S. A. Hires, T. Mao, D. Huber, M. E. Chiappe, S. H. Chalasani, L. Pe-treanu, J. Akerboom, S. A. McKinney, E. R. Schreiter, C. I. Bargmann, V. Ja-yaraman, K. Svoboda, and L. L. Looger. Imaging neural activity in worms, fliesand mice with improved GCaMP calcium indicators. Nature Methods, 6:875–881,2009.[137] Y. A. Tolias, S. M. Panas, and L. H. Tsoukalas. Generalized fuzzy indices forsimilarity matching. Fuzzy Sets and Systems, 120(2):255–270, 2001.149Bibliography[138] R. Y. Tsien. New calcium indicators and buffers with high selectivity againstmagnesium and protons: design, synthesis, and properties of prototype struc-tures. Biochemistry, 19(11):2396–2404, 1980.[139] R. Tur, Y. C. Eldar, and Z. Friedman. Innovation rate sampling of pulse streamswith application to ultrasound imaging. IEEE Transactions on Signal Processing,59(4):1827–1842, 2011.[140] M. Unser and T. Blu. Cardinal exponential splines: part I - theory and filteringalgorithms. IEEE Transactions on Signal Processing, 53(4):1425–1438, 2005.[141] J. A. Urigüen, T. Blu, and P. L. Dragotti. FRI sampling with arbitrary kernels.IEEE Transactions on Signal Processing, 61(21):5310–5323, 2013.[142] J. A. Urigüen Garaizabal. Exact and approximate Strang-Fix conditions to re-construct signals with finite rate of innovation from samples taken with arbitrarykernels. PhD thesis, Imperial College London, 2013.[143] I. Valmianski, A. Y. Shih, J. D. Driscoll, D. W. Matthews, Y. Freund, andD. Kleinfeld. Automatic identification of fluorescently labeled brain cells forrapid functional imaging. Journal of Neurophysiology, 104(3):1803–1811, 2010.[144] M. Van Rossum. A novel spike distance. Neural Computation, 13(4):751–763,2001.[145] M. Vetterli, P. Marziliano, and T. Blu. Sampling signals with finite rate ofinnovation. IEEE Transactions on Signal Processing, 50(6):1417–1428, 2002.[146] I. Vlahinić, E. Andò, G. Viggiani, and J. E. Andrade. Towards a more accu-rate characterization of granular media: extracting quantitative descriptors fromtomographic images. Granular Matter, 16(1):9–21, 2014.[147] J. T. Vogelstein, A. M. Packer, T. A. Machado, T. Sippy, B. Babadi, R. Yuste,and L. Paninski. Fast nonnegative deconvolution for spike train inference frompopulation calcium imaging. Journal of Neurophysiology, 104(6):3691–3704,2010.[148] J. T. Vogelstein, B. O. Watson, A. M. Packer, R. Yuste, B. Jedynak, andL. Paninski. Spike inference from calcium imaging using sequential monte carlomethods. Biophysical Journal, 97(2):636–655, 2009.150Bibliography[149] C. Vonesch, T. Blu, and M. Unser. Generalized Daubechies wavelet families.IEEE Transactions on Signal Processing, 55(9):4415–4429, 2007.[150] G. Wang, D. R. Wyskiel, W. Yang, Y. Wang, L. C. Milbern, T. Lalanne, X. Jiang,Y. Shen, Q.-Q. Sun, and J. J. Zhu. An optogenetics- and imaging-assisted si-multaneous multiple patch-clamp recording system for decoding complex neuralcircuits. Nature Protocols, 10(3):397–412, 2015.[151] B. Weber and F. Helmchen. Optical Imaging of Neocortical Dynamics. Neu-romethods. Humana Press, 2014.[152] J. M. Whittaker. The Fourier theory of the cardinal functions. Proceedings ofthe Edinburgh Mathematical Society, 1:169–176, 1928.[153] E. Yaksi and R. W. Friedrich. Reconstruction of firing rate changes across neu-ronal populations by temporally deconvolved ca2+ imaging. Nature Methods, 3:377–383, 2006.[154] W. Yang, J. K. Miller, L. Carrillo-Reid, E. Pnevmatikakis, L. Paninski, R. Yuste,and D. S. Peterka. Simultaneous multi-plane imaging of neural circuits. Neuron,89(2):269–284, 2016.[155] L. A. Zadeh. Fuzzy sets. Information and Control, 8(3):338–353, 1965.[156] C. Zimmer and J. C. Olivo-Marin. Coupled parametric active contours. IEEETransactions on Pattern Analysis and Machine Intelligence, 27(11):1838–1842,2005.[157] H.-J. Zimmermann. Fuzzy set theory. Wiley Interdisciplinary Reviews: Compu-tational Statistics, 2(3):317–332, 2010.[158] J. Ziv and M. Zakai. Some lower bounds on signal parameter estimation. IEEETransactions on Information Theory, 15(3):386–391, 1969.[159] K. H. Zou, S. K. Warfield, A. Bharatha, C. M. Tempany, M. R. Kaus, S. J.Haker, W. M. Wells, F. A. Jolesz, and R. Kikinis. Statistical validation ofimage segmentation quality based on a spatial overlap index: Scientific Reports.Academic Radiology, 11(2):178–189, 2004.151

Inferring Spike Trains from Calcium Imaging

Abstract

Extracted data

Inferring Spike Trains from Calcium Imaging

Abstract

Extracted data

Related items

Related items