# R/findEpisodes.R In sjemea: Multi electrode analysis

#### Defines functions extract.burst.times.with.idburstmatrix.to.episodefindEpisodes

```## SJE version of finding episodes by overlapping bursts.
## 2011-06-09

findEpisodes <- function(s, max.burst.length=40) {
## Find the episodes of overlapping bursting activity.
##
## MAX.BURST.LENGTH is the longest allowed burst.
## Any burst longer than this is ignored. Set this value to "Inf" if
## you do not want any filtering of burst times.

## First, perform burst analysis (could check if not already done.)
if (is.null(s\$allb)) {
s\$allb <- lapply(s\$spikes, mi.find.bursts)
}
allb <- s\$allb

## Find the number of bursts in each channel; if a channel has no bursts, just
## let it be of length 1.
nbursts <- sapply(s\$allb, function(m) {n <- nrow(m); ifelse(is.null(n), 1, n)})
rep(names(nbursts), times=nbursts)
burst.ids <- as.numeric(rep(names(nbursts), times=nbursts))

## Flatten the "allb" burst information into one big matrix, to
## create the BURSTS matrix.
## Each row of the burst matrix contains three columns (beg,end, id):
## BEG, END: Spike number of the first and last spike in the burst.
## ID: channel number from where this burst came from

allb.flat <- do.call(rbind, s\$allb)
allb.flat2 <- cbind(allb.flat[,1:2], burst.ids)

## Convert the spike number into a spike time, by looking up times
## in the spike trains for each channel.

all.ids <- mapply(extract.burst.times.with.id,
s\$spikes, allb, as.numeric(names(allb)))
flat <- do.call(rbind, all.ids)

## remove any electrodes which had no bursts.
zero.bursts <- which(is.na(flat[,1]))
if (any(zero.bursts))
flat <- flat[-zero.bursts,]
colnames(flat) <- c("beg", "end", "id")

## Remove any bursts which are too long.
long.bursts <- apply(flat, 1, function(x) { (x["end"] - x["beg"]) > max.burst.length})
if (any(long.bursts)) {
printf("%d long bursts removed\n", sum(long.bursts))
flat <- flat[!long.bursts,]
}

## Now that we have the BURSTS information in the right format, we
## can now find the episodes of overlapping bursts.

if (nrow(flat)==0 || is.null(flat)) {
## We have no bursts to make episodes from
## (e.g. if max.burst.length is too small).
episodes <- NULL
} else {
episodes <- burstmatrix.to.episode(flat)
}

episodes
}

burstmatrix.to.episode <- function(bursts) {
## Given a BURST matrix, sort it according to the burst start time
## and divide it up into smaller episodes, such that each block
## represents the bursts within one episode.

## An episode is defined as all the bursts that overlap in time;
## e.g. if burst A and B overlap and burst B and C overlap, even
## though A and C may not overlap, A,B,C form an episode.
##
## By definition, if one electrode fires and does not overlap with
## any other burst, that counts as an episode, with just one burst.
##
##

## For each episode, we then compute a few statistics, such as the
## start and end time of the episode, and the number of electrodes
## recruited.

get.stats <- function(be) {
burstinfo <- bursts[be[1]:be[2],,drop=FALSE]
beg <- burstinfo[1,"beg"]
end <- max(burstinfo[,"end"])
num.channels <- length(unique( burstinfo[,"id"]))
c(beg=beg, end=end, num.channels=num.channels)
}

## Sort the burst information according to the burst onset time.
## This makes episode detection straightforward: bursts i and i+1 do
## not overlap (and hence do not belong to the same episode) if the
## start time for burst i+1 is greater than the end time for burst i.

bursts <- bursts[order(bursts[,"beg"]), ]

B <- nrow(bursts)

new.episode <- TRUE
looking <- TRUE
i <- 1

## worst case: each burst is on its own, so max number
## of episodes is equal to number of bursts
breaks <- rep(NA, B)                    #store start of new episode
##breaks[1] <- 1
j <- 1
end <- bursts[1, "end"]

while (looking) {

if ( bursts[i+1, "beg"] <= end) {
## keep the chain going: possible extend end of episode.
end = max(end, bursts[i+1, "end"])
} else {
## come to the end of an episode.
breaks[j] <- i
j <- j + 1
end <- bursts[i+1, "end"]
}

i <- i + 1
if (i == B) {
## finished all the bursts.
breaks[j] <- i
looking <- FALSE
}
}

ends <- breaks[1:j]
starts <- c(1, 1+ends[-length(ends)])
episode.begend <- cbind(starts, ends)

episodes <- apply(episode.begend, 1, get.stats)
res <- list(beg=episodes[1,],
end=episodes[2,],
num.channels=episodes[3,])
res
}

extract.burst.times.with.id <- function(spikes, burstinfo, id) {
## Given the beg+end indexes of a burst, lookup the times of the burst
## and return this along with the channel name.
nbursts <- nrow(burstinfo)
if(!is.null(nbursts)) {
begt = spikes[burstinfo[,"beg"]]
endt = spikes[burstinfo[,"end"]]
id2 = rep(id, nbursts)
}
else {
## no bursts
begt = NA
endt = NA
id2 = id
}
cbind(begt, endt, id2)
}
```

## Try the sjemea package in your browser

Any scripts or data that you put into this service are public.

sjemea documentation built on May 31, 2017, 3:45 a.m.