Modeling | MATLAB Workshop

Poisson Simulation of Neural Data: Building a Simple Model

Building models can help us simulate and predict the behavior of neurons. Here, we're going to build a non-homogenous Poisson model of a neuron, fit this model to a dataset, and generate more spike trains to augment the dataset. This is a relatively simple model that predicts neuronal firing based on the assumption that spiking is stochastic (random) and driven by some time-variant function r(t). This model was proposed by Siebert based on observations from the auditory system (paper here).

Poisson distributions describe the probability (P) of an event with a known rate happening during a given time period, with the assumption that each event is independent of any events before or after it (while this is obviously not fully true of a neuron, we're going to write a simple model assuming that it is). Poisson distributions are described by one parameter, λ, that describes the mean number of events during a given time period. For us, this will be r(t).

The data used for fitting were published as supplemental materials to MATLAB for Neuroscientists (Wallisch et al.).

* I want to caveat this section by saying that I am not a MATLAB modeler (or any type of modeler). Much of this is based on this tutorial, and you'll learn much more useful things about the math in Neuro 527 and Neuro 545. The focus here is still on learning to use MATLAB effectively, not really on generating a super useful model.

MATLAB Concepts Covered:

using random number generators, creating logical arrays, using other people's code, performing sanity checks on your data, using sinusoids, writing functions

Beginner Project: Non-Homogeneous Poisson Process to Predict Neural Firing

Because we're not starting with data, this project will follow a different pattern than the others. I'm not certain how well this format will generalize to more complicated model-building, it has just been the way I've found it the easiest to approach basic modeling...

1. Build and test a simple version of the model.

We're going to start with a very simple, time-invariant Poisson model that assumes a constant firing rate. Again, the point here is to learn MATLAB, so just trust these equations for the time being (look at the tutorial linked above if you want more on the math). Let's begin with a 10 second time frame, which we will assume has a constant firing rate of 50 spikes/second. We're going to "sample" this time frame at intervals of 1 millisecond and generate a spike train. The equation we want to implement is:

P{1 spike during δ} = rδt

(the probability of a spike during a given time interval equals the firing rate times the length of the time interval)

Based on the above, what are we plugging in as r and δt? (click for answer)

r = 50 spikes/sec, δt = 0.01 sec

We want a yes/no (1/0) answer to the question "did the neuron fire at this time point" for each time point that we sample from, with the result adhering to the probability computed above. Any ideas of how to could accomplish this? Brainstorm, then click to continue.

* We're going to generate a random number between 0 and 1 for each time point. If it is less than or equal to the probability value (P). we'll say the neuron spiked, and if it is greater than P, we'll say it didn't spike. (Make sure you understand how this relates to the probability).

What built-in function generates uniformly distributed random numbers between 0 and 1, and how many random numbers do we need to generate? (click for answer)

rand(), 10000

* Create a logical array to represent the spike train. Plot it using the plotRaster() function to visualize it.

How many spikes per second do you have on average? (click for answer)

sum your logical and divide by 10

* Now let's generate 10000 spike trains with a firing rate of 15 spikes/second, a δt of 1 msec, and a length of 1 second. Plot them all using plotRaster(). Do you see any patterns? Probably not, because the firing rate is invariant and entirely stochastic. Plot a histogram of the number of spikes per second over these 10000 trials.

What's special about this histogram? (click for answer)

It is randomly distributed around a mean of 15 (the firing rate)

* If you got the answer above, things are going correctly! Our goal was to create a model in which spike timing was randomly determined, but driven by a set firing rate. Test this out with a few other firing rates and lengths of time to make sure that the same patterns continue to emerge and you're confident in your understanding of the basic probability principles you're using here.

2. Identify inaccuracies in your model that will matter for you to achieve your goals.

I'll give you one for free: the time invariance is going to be a problem! Any other properties of neurons that you might foresee being an issue? Click for my list, but there are definitely more inaccuracies of this model!

Time invariant firing rate
No accounting for bursting behavior
No refractory period
The firing rate has an upper cap based on sampling frequency
Need to know firing rate of your neuron

We're going to address the first one or two items on the list above in the next step.

3. Address these inaccuracies.

First, we'll address the time-invariant firing rate. (click for more)

Instead of having a fixed r, we're going to make r into a function, r(t). Our model is going to assume that we know the time course of the firing rate, which will be most useful when you have existing data and would like to replicate that data with a model rather than when you're trying to simulate data based on, for instance, an input current.

* Let's start with a hypothetical neuron that has a resting firing rate of 50 spikes/second that increases immediately to 100 spikes/second in the presence of a stimulus. Create a variable, r_t, that represents the firing rate of this neuron over a 1 second time course with 500 ms of rest followed by 500 ms of stimulus, sampled every millisecond. Plot r_t to make sure it looks like a step function.

* Use the same protocol above to generate a spike train, choosing a number between 0 and 1 randomly for each time interval and comparing it to the probability of a spike during that interval, but with the new equation:

P{1 spike during δ} = r(t)δt

Click for a hint.

When the <= function is given two vectors of the same dimension, it will compare elements pairwise:

a <= b will produce a vector of [a(1) <= b(1), a(2) <= b(2), ..., a(n) <= b(n)]

* Calculate the firing rate in your spike train for the first 500 ms and the last 500 ms. How do they compare to your expectations? Then generate a total of 10000 spike trains and plot histograms of the average firing rates during those same time ranges. How do these plots compare to your expectations?

What's the average probability of a spike in a given time bin during the first 500 msec? During the second 500 msec? Why does this make sense? (click for answer)

Approximately 5% for the first half, 10% for the second half.

This makes sense because it is the given firing rate times dt

Click for a hint.

The mean of the first 500 msec along the trials dimension will give you the probability of firing for each of the first 500 msec. You want the average probability of firing over this entire time period.

* Now let's model a scenario where a neuron's firing rate varies sinusoidally. Make a variable r_t that represents 5 seconds of a firing rate time course sampled every millisecond. The firing rate should vary at a rate of 5 Hz, the amplitude should be 25, and the minimum firing rate should be 25 spikes/second.

Click for a hint.

The sin() function requires an input, so you're going to want to create a time course (in seconds) to feed into the function. Do some Googling to try to figure out how to make a vector that increments by .001 seconds (1 msec) and represents the full 5 seconds. If you don't want to have to remember precalc to figure out what else to feed into the function, click here (but makes sure you know what each part means!).

Keywords to Google: colon, vector

t = 0.001:0.001:5; % this is equivalent to t = (1:5000)/0.001: from 1 msec to 5 sec with

% an interval of 1 msec

r_t = 25*sin(10*pi*t) + 50; % plot to make sure it fits the criteria outlined above, and

% make sure you know where the 10*pi comes from

* Generate 10000 spike trains and use plotRaster() to visualize them. Does it look as you expected it to?

* If you're still unconvinced that you understand how r(t) relates to the final rasters or histograms, try out a couple of other r(t)s. What about a firing rate that increases linearly over the course of 2 seconds? Or a much shorter step function?

If it is of interest, you can also add bursting behavior (although this won't be used to answer our main question). (click for more)

In our current model, the interval between spikes is totally random, when in reality there might be an increased probability of having multiple events in close succession. We can fix this by adding "burstiness" to our neuron.

* Let's return to a time-invariant scenario to simplify things for now, where the firing rate is 50 spikes/second throughout a 1 second interval sampled every millisecond. Generate 1000 spike trains, as before, and make sure the rasters and histograms appear as you would expect. Note that this will be a smaller array than you had during part 1: generating bursts requires more processing power and takes more time, so we're reducing the number of points

* Now, instead of accepting each 1 in our logical array as an individual spike, we're going to re-interpret it as an event. On average, an event will consist of 1 spike, but it could stochastically have more or fewer spikes.

What type of function could represent this distribution? (click for answer)

another Poisson function!

* Find the command to plot a cumulative density function (CDF) for the type of distribution above. Plot values for 0-10 with a mean of 1. Make sure you understand (in general terms) what a CDF is.

* To plot the CDF, you gave it x values, which for us represent the number of events. What we want to do is to take a probability value, which we will randomly generate, and pull its corresponding x value.

What MATLAB function can do this? (click for answer)

poissinv() - the inverse Poisson CDF

* Create an array of random numbers that is the same size as your logical array of 1000 trials. Feed this into the function above (again using <= and > to determine if random numbers pass the probability) and find the mean number of events in each trial. Plot a histogram of these values. Does it look as you expect it to?

* Replace each 1 in your event array with the number of events, as determined above. Try to plot it with plotRaster() - you'll get an error.

Why are you getting this error? (click for answer)

your array now has numbers instead of being all 0s and 1s

* Instead, create a bar graph of one trial, comparing your initial spike train to your new "event train." Also compare the histograms of numbers of spikes per trial with and without bursting (you might need to adjust the edges variable of your histograms to make them directly comparable). Do these plots look as you'd expect?

What is similar and different between the histograms of average number of events per trial pre- and post-bursting modifications? (click for answer)

same mean, different variances (you canconfirm this computationally)

* If you'd like, try mixing this with a time invariant firing rate such as the ones demonstrated above. Do your results still align with your expectations?

4. Formalize the model and make an easy to use version.

While this model is still certainly flawed, this is as complex as we're going to make it for the time being. Hopefully you've been keeping track of your work in a script of some sort, but as of now you have to repeatedly run multiple lines of code every time you want to make a change to your model. To make the model easier to use, we're going to write a function that accepts key parameters, performs all necessary calculations, and returns event trains (look at this documentation if you need a function refresher). What parameters do you think your function should accept? Try to create a complete list. Click to see the parameters that my function takes.

t: a time course that defines each time point (i.e. .001:.001:1 for 1 second sampled every msec). You could also accept two different variables - the length of the trial and dt, but by using length() and diff() you can get the information you need from just t.
r_t: a time course of firing rate that is the same length as t.
numTrials: the number of trials you would like to simulate.
[meanEvents: the mean number of events per potential spike (in the above examples, this was always 1). - you'll want to make a function without bursting as well for use in part 5]

5. Test how well your model fits an actual dataset.

Now that we can easily generate event trains, we're going to try to fit our model to an existing (small) dataset.

* Open MT.mat, which includes a variable t that corresponds to the time course of the trials, and a logical array spks that includes multiple trials of spike trains recorded from a cell in visual area MT. During each trial, 4 seconds of rest was recorded, followed by .5 seconds in the presence of a visual stimulus.

How many trials are included here? (click for answer)

What are the average firing rates during rest and during stimulus presentation? (click for answer)

~41 and 85 spikes/sec

* Using the given information, generate 10000 more simulated trials. Compare the average firing rates to the actual dataset. They should be quite similar (because this is exactly what you fed the model).

* Compare the variances during each period. What do you notice and why do you think this is?

Finally, write a function to determine the average interspike interval, and use it on both your simulated data and the given data. What do you notice and why do you think this is?

Click for hint.

You will need the functions find(), unique(), and diff().

Click for my interpretation of what's going on here.

The variance decreases a lot because, in reality, the firing rate of neurons is not going to be that stable over time. The interspike interval slightly increases, which I would expect reflects the existence of some sort of bursting behavior in the recorded cell.

Overall, a model like this is going to have limited utility because we don't usually know raw information about firing rates (and even if we do, we can't usually make predictions for other states' firing rates). A more useful version of this model would not require us to input a firing rate, but perhaps to input a current that is then translated to the firing rate using a membrane voltage model. You'll make and use models more like that during your first year coursework!

Advanced Project Ideas

If you clicked on this, you likely have better ideas about advanced modeling projects that would be fun than I do! Maybe try playing with this tool? Make a network of Poisson firing neurons and see how they behave? Make a function that will calculate the firing rate based on inputted current rather than requiring a known firing rate? If you have ideas of what to add to this part of the site please let me know :)