SpikeOMatic Tutorial 2 SpikeOMatic goes MCMC

Matthieu Delescluse, Christophe Pouzat

Description

SpikeOMatic Tutorial 2
SpikeOMatic goes MCMC
Matthieu Delescluse and Christophe Pouzat
July 1, 2005Contents
1 Introduction 2
1.1 Why does SpikeOMatic go MCMC? . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Data generation model: the Dynamic Hidden Markov Model . . . . . . . . . . . . 2
1.2.1 Spike amplitude relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.2 Neuronal discharge statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 The MCMC approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.2 Con guration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3.3 Posterior density and the MCMC method . . . . . . . . . . . . . . . . . . 4
1.3.4 Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.5 The Replica Exchange Method . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3.6 Last check before we start . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Re-generate this tutorial 6
3 Start and load data 6
3.1 Loading the SpikeOMatic Package in R . . . . . . . . . . . . . . . . . . . . . . . . 6
3.2 Reading the data into R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.3 Detecting spikes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.4 Making spikes and noise sweeps . . . . . . . . . . . . . . . ...

Sujets

• Neurone
• Energy de l'Iowa
• Tutoriel
• Vecteur
• Amplitude
• Orbite cristallographique

Informations

• Publié par: Astlyo
• Nombre de lectures: 164
• Langue: English
• Poids de l'ouvrage: 2 Mo

Extrait

SpikeOMatic Tutorial 2 SpikeOMatic goes MCMC

Matthieu

Delescluse

July

and Christophe

1,

2005

Pouzat

Contents 1 Introduction 1.1 Why does SpikeOMatic goMCMC?. . . . . . . . . . . . . . .. . . . . . . . . . . 1.2 Data generation model: the Dynamic Hidden Markov Model. . . . . . . . . . . . 1.2.1 Spike amplitude relaxation. . . . . . . . . . . . . . .. . . . . . . . . . . . 1.2.2 Neuronal discharge statistics. . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 TheMCMCapproach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Model parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Configuration. . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . 1.3.3 Posterior density and theMCMCmethod. . . . . . . . . .. . . . . . . . 1.3.4 Energy. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . 1.3.5 TheReplica Exchange Method. . . . . . . . . . . . . . . . . . . . . . . . . 1.3.6 Last check before we start. . . . . . . . . . . . . . . . . . . . . . . . . . . Re-generate this tutorial Start and load data 3.1 Loading the SpikeOMatic Package in R. . . . . . . . . . . . .. . . . . . . . . . . 3.2 Reading the data into R. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Detecting spikes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Making spikes and noise sweeps. . . . . . . . . . . . . . . .. . . . . . . . . . . . 3.5 Whitening. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . 3.6 Gaussian mixture based spike-sorting. . . . . . . . . . . . . .. . . . . . . . . . . DHMM based spike-sorting 4.1 Building theMCMCinputs. . . . . . . . . . . . . . . . .. . . . . . . . . . . . . 4.2 Running theMCMCalgorithm. . . . . . . . . . . . . . . .. . . . . . . . . . . . 4.3 Output files. . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . 4.4 Reading the output files. . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . 4.4.1 Reading the energy file. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Reading the configuration file. . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 Reading the neurons file. . . . . . . . . . . . . . . .. . . . . . . . . . . .

2 3 4

1

2 2 2 2 3 3 3 4 4 5 5 6 6 6 6 7 7 8 9 9 12 12 14 16 16 16 19 24

1 Introduction This tutorial gives details on how to achieve a complete analysis of raw multi-unit data using SpikeOMaticand theSpikeOMatic Library you have not read Ifin the R environment.SpikeOMatic Tutorial 1[4] yet it is a good idea to stop and read it for we will only briefly comment on the stuff which has already been exposed there. Another point you should consider before going further is: “Do I need this Markov Chain Monte Carlo (MCMC) stuff and what is it anyway?”

1.1 Why does SpikeOMatic goMCMC? We give only a brief explanation here for we have already spent a fair amount of time (i.e., paper length) exposing the point in a recent paper [6] and book chapter [3 want to consider using]. You theMCMCextension ofSpikeOMaticwhen the neurons you record from exhibit one or both of the following features: •depends on the inter-spike interval (The spike amplitude of the neurons ISI). For instance, one commonly observes during a burst a reduction of the spike amplitude upon shortISI. •You want to include in the spike-sorting procedure information about the discharge statistics of the neurons, because, among other things, you can then perform good sorting even for neurons whose clusters overlap strongly onWilson plots(see [3]). The capability do deal with such situations, at least our solution to these situations, comes at a price. We need to use an explicit data generation model which is significantly more sophisticated than the “simple” Gaussian mixture model (GMM) we used inSpikeOMatic Tutorial 1[4]. We will dub this model, theDynamic Hidden Markov Model(DHMM). The first term (Dynamic) refers to the spike amplitude dynamics, the next two terms (Hidden Markov) refer to the neuron’s discharge statistics. Statistical inference on theGMMcould be easily done with theExpectation- Maximizationalgorithm (EM), whereas theDHMMrequires the rather heavy computational machinery of theMCMCapproach.

1.2 Data generation model: the Dynamic Hidden Markov Model The Dynamic Hidden Markov Model we present here is our current compromise between model complexity (the number of parameters of the model) and its flexibility to account “reasonably” well for actual neuronal data. It is not designed to provide a perfect fit, but to capture essential features of the data. These features are the spike amplitude dynamics and a possibly dynamic (e.g., bursty) inter-spike interval generation.

1.2.1 Spike amplitude relaxation In what follows we assume that the relation between the amplitudeaof a spike (with no recording noise) and the elapsed time since the last spike of the same neuron, the inter-spike interval (isi), is given by an exponential relaxation:

a(isi|P=p,Δ

a(isi|P max =p max ,Δ =δ,Λ =λ) =p max ∙(1−δ∙exp(−isi∙λ)),(1) whereP max is the maximal possible amplitude of the spikes (that is, the one we would observe following a long enoughisia modulation factor (dimensionless) and Λ is the inverse of), Δ is the relaxation time constant (measured ins −1 ). Remark that we use vectors to represent spikes’ amplitudes. This is because we typically use several amplitude samples to represent a spike (e.g., in the sequel we will use the peak amplitude of the spikes on each of four recording sites, oura andp max objects  Thiswill therefore be vectors in a 4-dimensional space). relaxation model means 2

that for each neuron we will have to estimaten s + 2amplitude parameters(n s is the number of amplitude samples used to describe a spike, 4 in the sequel). We will now assume that our measured amplitudes are corrupted by aGaussian white noise (with a standard deviation of 1) which is statistically independent of the signal (the spikes). In other words we will have towhiten the noisebefore running our spike-sorting procedure [5] & [4] as explained in Sect.3.5 the probability (density) to observe a spike,. Thens, following the previous one (from the same neuron) byisi(seconds) becomes:

π amp (S=s|ISI=isi,P max =p max ,Δ =δ,Λ =λ) = (2π) − n 2 s ∙exp (− 1 2 ks−a(isi)k 2 ),(2) wherea(isi) is given by Eq.1andkvkrepresents the Euclidean norm of vectorv. We will try to consistently use symbolπ something to designate aproperly normalizedprobability (density) function. We will use capital letters to designaterandom variablesand small cap letters to designate theirrealizations.

1.2.2 Neuronal discharge statistics We will extend here our previous model [6] by introducing multipledischarge states,D=d, with d∈ {1, . . . , N ds }, for each neuron. The goal is to be able to describe at the same timemulti-modal ISIdensities (i.e., densities with several peaks) and dependence between successiveISIssimilar to what goes on during a burst, where (single) “silent periods” (longISIs) are followed by many short ISIs Camproux et al [. Following1] we will assume that successiveisiare independentconditioned on the neuron’s discharge state,d the emission of each spike (. Afteri.e.,the generation of eachisi) the neuron can change itsdischarge states The inter discharge state dynamicsor keep the same. is given by aMarkov matrix,Q= (q i,j ). We will moreover assume that theisidistribution of each state islog-normal other words we assume that the neuron after its. Inmth spike is in stated and that theisibetween this spike and the next one is distributed as:

ISI|d∼log-normal(s d , f d )

(3)

Eq.3should be read as: π isi (ISI=isi|S=s d , F=f d ) = isi∙f d 1 ∙ √ 2π ∙exp  − g(1l 2 o ( isi) f − d log(s d ) ) 2  (4) whereSis ascaleparameter (measured in s) andFis ashapeparameter (dimensionless). These random variables do depend on the value taken by the random variableD. After theisi has been generated, the neuron can “move” to any of itsN ds states according to a probabilistic dynamics described by: P(D m+1 =k|D m =j) =q j,k (5) You see therefore that if we work with a neuron with 3 discharge states we have 12 independent isi parameters pairs ( 2to estimates:s, f) per state andN ds ∙(N ds −1) state transition parameters (don’t forget that matrix (q i,j ) isstochastictherefore its rows sum to 1).and More details about this model will hopefully be soon readily available to everyone [2].

1.3 TheCMCMoachrppa 1.3.1 Model parameters You have probably noticed and you’ve perhaps been surprised by the fact that in Eq.1,2&4 we treat our model parameters: 3

{P max ,Δ,Λ, S j , F j ,(q i,j )}(6) likerandom variables for you, Luckily is because we are using a Bayesian approach.. That you do not have to take side in the religious war opposingfrequentistsandBayesianto use our algorithm. If you’re deeply convinced that the likelihood approach is the beginning and the end of any respectable statistical analysis you can view our use of the Bayesian framework as an easy way to introduce constraints on our model parameters. Considering as well that the nature of our problem would require a Monte Carlo algorithm (that is anMCMCmethod) to tackle the associated maximum likelihood problem. It just turns out that going fully Bayesian represents a slight extra computational cost. Fine, so let’s say we consider for now a model withKneurons. Our task includes the estimation of the parameters of every neuron. We will call Θ this complete parameter vector:

Θ = (P max,1 ,Δ 1 ,Λ 1 , S d,1 , F d,1 ,(q i,j ) 1 , . . . ,P max,K ,Δ K ,Λ K , S d,K , F d,K ,(q i,j ) K ) (7)

1.3.2 Configuration We are of course directly interested (not to say mainly interested) in the neuron of origin of each detected spike. Or more precisely, the posterior probability for each neuron of the model to have generated each spike. In order to do that, we will associate to every detected spike,j∈ {1, . . . , N}, 2 random variables:L j andD j |l j .L j takes value in{1, ..., K}and its distribution is what we are fundamentally looking for: the probability for each neuron of the model to have generated spike j.D j |l j is the discharge state we introduced in Sect.1.2.2(we are just making explicit here the fact thatD j depends on the value taken byL j ). We will therefore have to estimate the posterior distribution of what we will call theconfiguration:

C= (( DL ,), . . . ,(L , D))

C= ((L 1 , D 1 ), . . . ,(L N , D N )) (8) The fundamental reason to use a simulation based inference (that is, a Monte Carlo method) is the number of elements of the space on whichC (takes value:KN ds ) N . This becomes astronomical for realistic values:K= 10,N ds = 3 andN= 1000.

1.3.3 Posterior density and theMCMCmethod MCMCalgorithms solve these inference problem by generating a Markov Chain on the space defined by the problem’s unknowns. In our case the space on which (Θ, C) are defined. This Markov Chain is built such that it converges to auniquedistribution which is precisely the posterior distribution of the unknowns given the data (see [6] and [3 our case: In] for details). π post (θ, c| D) = L(D, c|θ) Z ∙π prior (θ) ,(9) whereDrepresents the data (spike amplitudes and occurrence time),L(D, c|θ) is the likeli- hood of the dataandthe configuration given the model parameters,π prior (θ) is the prior density of the model parameters andZis the normalizing constant: Z= ∈ X Z θ dθL(D, c|θ)π prior (θ),(10) cC whereCis the space on whichCis defined. Formally our Markov Chain tells us how we move from the present “system” state, (θ, c), to a new one, (θ 0 , c 0 ), by specifying a transition kernel:T((θ, c)→(θ 0 , c 0 )). What we do on our computer is a Monte Carlo simulation of this Markov Chain (that’s what this fancy acronym MCMC will therefore end up with a sequence of states:means). We 4

{(θ 0 , c 0 ),(θ 1 , c 1 ), . . . ,(θ M , c M )}(11) where (θ 0 , c 0 ) is the initial state andMis the number of Monte Carlo (MC) steps performed. Then if we want to get the most likely neuron of origin for spike 100 we compute the following quantity:

j∈{ m 1,. a .., x K} M 1 −m i= X m M δ l 1 i 00 ,j (12) whereδ l i 100 ,j is the Kroenecker’s symbol andmis the step at which we start “taking measure- ments”. That the number of steps we judge sufficient to forget about the initial state (see [6] and [3]).

1.3.4 Energy Now, as a user ofMCMCmethods you should be aware that the convergence theorems hold for aninfinite Onnumber of steps. your computer you can of course perform only afinitenumber of steps. You therefore have to make sure that the behavior of your simulated Markov Chain iscompatible In this tutorial we will do thatwith a chain at “steady-state” or “equilibrium”. qualitatively by monitoring the evolution of theenergyof the Markov Chain that we define by: E(θ, c) =−log (L(D, c|θ)π prior (θ)) (13) You can qualitatively see E as aχ 2 it is, the better the fit. Of course,value, that is, the smaller if the Markov Chain has reached equilibrium, E should oscillate around a constant mean value.But be careful, for this is a necessary, not a sufficient condition for equilibrium.

1.3.5 TheReplica Exchange Method One more “theoretical” point before we can start playing with the algorithm: the energy landscape of our systems are typically rough and our Markov Chain dynamics is local. That means we can get traped in local energy minima for a long time. To circumvent this problem we have recourse to theReplica Exchange Method(REM) which consists in simulating multiple replica of our system at higher and higher “temperatures”, so that the replica at high temperature can escape from local minima (see [6] and [3] for details). The temperature parameter, or rather the inverse temperature parameter,β, will show up in the following way: A Markov Chain simulated at a single fixedβ≥0 will have its stationary density given by:

π post,β (θ, c| D  ( exp ) = − Z β (β ∙ ) E(θ, c)) ,(14) where E(θ, c) is defined by Eq.13. The symbolβand the formalism we use here come from statistical physics, where:β= (kT −1 ), with, k, the Boltzmann constant and T, the absolute temperature. What we will do in fact is simulating a Markov Chain on an extended space. If we callPthe space on which Θ is defined (we already usedCfor the space on whichCis defined) and if we use N β different temperature ((β 1 , . . . , β N β )), then our new Markov Chain is defined on: N β Y (P × C) (15) p=1

and its stationary density is given by:

5

N β Y π post,β p (16) p=1 To get the statistics of interest we then just need to use the sub-chain atβ= 1. 1.3.6 Last check before we start In order to perform anMCMCbased analysis we will have to call a stand alone C function som_mcmc. We first need to find out if the function is accessible from the directory where we are now running R. We can make this check as follows: > system("som_mcmc --version")

2 Re-generate this tutorial This tutorial was “auto-generated” from fileSpikeOMaticTutorial2.Rnwusing functionSweave. You can re-generate it on your own machine. That will allow you to test that everything works as it should and give you an idea of the time required to run an analysis. To re-generate the tutorial on your computer, enter the following commands at the R prompt 1 (i.e., within the R environment): > library(tools) > vigdir <- system.file(’doc’, package = ’SpikeOMatic’) > vig <- list_files_with_type(vigdir, ’vignette’) > Sweave(vig[2]) That will generate some activity, figures will pop up and a whole analysis will be performed. By the end of it you will find new files in the directory from which R was running. These files are mostlypdfandpsversions of the figures you just generated. There will be also a SpikeOMaticTutorial2.texthat you can process to get apdfversion by typing the following commandsat the shell level(not within R): pdflatex SpikeOMaticTutorial2.tex makeindex SpikeOMaticTutorial2.idx pdflatex SpikeOMaticTutorial2.tex

3 Start and load data 3.1 Loading the SpikeOMatic Package in R Once in R, if you want to use a specific package like SpikeOMatic, you must load it using: > library(SpikeOMatic) You can then read your data into the R environment. 1 The following commands are correct for R version 2 or greater, if you’re using R 1.9.x you should typelist- FilesWithTypeinstead oflist_files_with_type.

6

3.2 Reading the data into R The case we will illustrate in this tutorial comes from a linear probe recording (see [2] and Fig. 3 of [3 four Thesettings) recording from young rat cerebellar slice.] for info about the recording recording sites of the probe (which are 50µm apart) were positioned along the Purkinje cell layer. We will therefore be dealing with 4 data files:PK_1.fl.gz, PK_2.fl.gz, PK_3.fl.gz, PK_4.fl.gz /usr/local/lib/R/library/SpikeOMatic/doc.. On your machine they are located in folder: 58 seconds of recording are used. Data were sampled at 15 KHz and filtered between 300 and 5K Hz. The data are compressed and in float format. Loading the data can be done interactively by first calling functionselect.data.files: > data.names <- select.data.files() which makes a window pop up inviting you to select the data files that can then be read with read.data.files(see below). It is also possible to manually build a character vector like thedata.namesabove as follows (we are proceeding here in that way because we want this tutorial to be self generative): > vigdir <- system.file("doc", package = "SpikeOMatic") > data1.name <- paste(vigdir, "/PK_1.fl.gz", sep = "") > data2.name <- paste(vigdir, "/PK_2.fl.gz", sep = "") > data3.name <- paste(vigdir, "/PK_3.fl.gz", sep = "") > data4.name <- paste(vigdir, "/PK_4.fl.gz", sep = "") > data.names <- c(data1.name, data2.name, data3.name, data4.name) After selecting the data in one or the other way, we can read them by callingread.data.files: > data <- read.data.files(trace.names = data.names, trace.length = 3e+05, + size.of = 4) We can visualize the data in an interactive way withdisplay.raw.data.gui: > display.raw.data.gui(data)

3.3 Detecting spikes Spikes are detected with functionfind.spikes.with.template, which uses a self generated spike template orfind.spikes.no.templatewhich does not (see Tutorial 1 [4]). We will use here the template based detection with a fairly large detection threshold to avoid having too many detected events (you are free to try other detection settings and you should do it): > analyze <- find.spikes.with.template(trace.matrix = data, threshold = 5, + minimal.distance = 5) In that case, 979 spikes are detected. You should at that stage check detection quality with functiondisplay.detection.gui: > display.detection.gui(analyze.with.template) which makes a GUI pop up allowing you to compare the raw trace and events detection.

7

Figure 1: A quick display of the spike sweeps for spikes detected with template. 200 (that’s the default value of the argumentnb.to.display scale) of the 979 detected spikes are shown. The bar on the central panel (Mean event scale bar on the The) goes from 0 to 1 in noise SD units. lower panel (SD) goes from 1 to 2 in noise SD units.

3.4 Making spikes and noise sweeps We can now extract spikes sweeps with functionmake.sweeps: > make.sweeps(spikes = analyze, sweep.length = 45, peak.location = 20, + nb.noise.evt = 2000) The spike sweeps we just generated can be visualized with functiondisplay.spike.sweeps: > display.spike.sweeps(spikes = analyze, main.top = "Spike detected with template") Which generates Fig.1. Functionmake.sweepsalso extracts 2000 noise sweeps in between the spike sweeps. Estimates of the noise auto- and cross-correlation functions can be obtained and displayed with function display.noise.correlation: > display.noise.correlation(noise = analyze) The result is shown on Fig.2.

8

Figure 2: Noise auto- and cross-correlation functions. The auto-correlation functions on each recording site are displayed on the diagonal of the plot. The cross-correlation functions are dis- played in an (hopefully obvious way) on the upper diagonal part.

3.5 Whitening Now that spikes and noise sweeps have been obtained, the spikes sweeps can be reduced and whitened. We first need to specify which amplitude samples of each spike sweep will be kept, by creating a vector of indices. In this analysis, we keep the peak amplitude on each recording site, i.ethe 20th point of the spike sweep on each site (recall that the length of a spike sweep on one site is 45 points): > selected.points <- c(20, 65, 110, 155) We can now call functionreduce.and.whiten: > reduce.and.whiten(spike.list = analyze, automatic = selected.points, + nb.samples.per.site = 1) reduce.and.whitennot only reduces and whitens the spike sweeps but generates also diagnostic plots as shown on Fig.3(for details see the on-line help and Tutorial 1 [4]). A Wilson plot of the reduced and whitened sample can be generated withdisplay.wilsonas shown on Fig.4. > display.wilson(analyze, , , , main = "Wilson plot with noise whitening")

3.6 Gaussian mixture based spike-sorting At this stage, the spike-sorting can be performed with Expectation-Maximization algorithm for Gaussian mixture models,GMM, (as detailed in Tutorial 1 [4]) by calling successivelly functions get.modelandclassify.from.mixture. We will force hereget.modelto work with a model made of 6 neurons. You are of course encouraged to try out models with more or less neurons (as 9

Figure 3: Comparison between theoretical and empirical noise properties. The upper panel shows the cumulative distribution of theMahalanobisdistances of an actual noise sample (reduced and whitened), in black, and the expectedχ 2  The lower distribution, in red.panel is a tougher test of the difference between the two distributions. The presence of more largeMahalanobisdistances (upper bend of the black curve ) than expected is mainly due to the presence of sub- detection threshold spikes.

10