How do I get a raster plot of nonhomogeneous spikes if I have a sequence of spike times over the total duration, in R?

Question

I am new to R (relatively) so this question might have a simple answer. I have been trying to plot non homogeneous Poisson Process in R. Through the code below, for one trial of 10 seconds, I have a sequence nhpp1 which consists of the time stamps where spikes occurred, in this particular trial. How do I take this sequence nhpp1 and get a raster plot of it. More importantly, I want to repeat (replicate?) this whole thing for 10 trials and get a raster plot which looks something like this: (please look below the code)

 nhpp <- function(lambda){
 set.seed(1)
 t_max = 10
 t = 0 
 s = 0
 Lambda <- function(tupper) integrate(f=lambda, lower =0, upper= 
           tupper)$value
 LambdaInv <- function(s) {
         v <- seq(0, t_max+1, length=1000)
         min(v[which(Vectorize(Lambda)(v) >= s)])
           }
 X = numeric(0)
 while(t <= t_max){
  u <- runif(1)
  s <- s-log(u)
  t <- LambdaInv(s)
  X <- c(X,t)

  }
 return(X)
}

lambda <- function(t) 100*(sin(pi*t)+1)
nhpp1 <- nhpp(lambda)

I have the spike time stamps already, I need help regarding finding a way to plot this one trial (with tiny bars occurring where spikes happen on the timeline) and how to replicate this process for 10 trials, then? Any help will be very appreciated.

Answer 1

I use another spike train generator based on the Poisson model.

arrival_time_v3 <- function(firing_rate,tMax,sampling_rate){
lambda <- firing_rate/sampling_rate

## find the number 'n' of exponential r.vs required by imposing that
## Pr{N(t) <= n} <= 1 - eps for a small 'eps'
n <- qpois(1 - 1e-8, lambda = lambda * tMax)

## simulate exponential interarrivals the
X <- rexp(n = 2*n, rate = lambda)
S <- c(0, cumsum(X))

arr_time <- S[which(S <= tMax)]
arr_time <- as.integer(arr_time)
arr_time <- arr_time[which(arr_time!=0)]
arr_time <- unique(arr_time)
return(arr_time)
}

num_of_spike_trains <- 10
firing_rate_arr <- matrix(c(10,20,30,40,50,60,70,80,90,100),1,num_of_spike_trains)
sampling_rate <- 10000
durartion_sample <- sampling_rate*10 ##10 sec

spike_arrival_time_mat_list <- list()
for(i in seq(1,num_of_spike_trains,1)){
spike_arrival_time_mat_list[[i]] <- t(as.matrix(arrival_time_v3(firing_rate_arr[i],durartion_sample,sampling_rate)))
}

After generating 10 spike trains (events trains), we can utilize barcode package in R:

#install.packages("barcode")

library(barcode)
barcode(spike_ariv_time_mat_list,xlab="Time",main="Spike Trains")