Methods in Computational Neuroscience Obidos, 2004 
Assembled, edited and written by Oren Shriki, Oren Farber and Dmitri Bibitchkov.
MATLAB is a generalpurpose mathematical software that is user friendly and very useful for data analysis, simulations, presentation of results, and more. The software runs on both UNIX and WINDOWS.
MATLAB is an interpreter,
which means that it performs every command line immediately after it is
entered. So you can either enter your commands line by line (prompt mode), or
prepare a script file first, and then run a program. A script file for MATLAB
is called an mfile and it must have the extension 'm' (i.e. something
like myprogram.m). To run an mfile you just enter the name of the file without the extension at the MATLAB prompt. Graphical output is available to supplement
numerical results.
For
those who prefer to find out everything themselves!
Online help is available from the Matlab prompt (a
double arrow), both generally (listing all available commands):
>> help
[a long list of help topics follows]
and for specific commands:
>> help fft
[a help message on the fft function
follows in the prompt window]
You can also get help on a function or a topic in a separate window using a command helpwin For example,
>> helpwin datafun
opens a desktop help window with information about functions for data analysis and Fourier transforms. At the end of each text concerning a function you can find links to related functions, indicated by 'See Also'. For general help, for example if you would like to get a general introduction into a specific Toolbox, use the command
>>
helpdesk
This
will display the start page for the online documentation in a MATLAB help
browser or in your web browser, depending on the settings of your system. The
online documentation is also available from http://www.mathworks.com/access/helpdesk/help/helpdesk.html
If you want to search for
commands related to some topic, say, complex numbers, you enter
>> lookfor complex
If you are a firsttime user,
you should invoke MATLAB now and follow along. In MATLAB you don't have to
define a variable type  you just enter it. The main object of the MATLAB is the matrix (The name MATLAB stands
for "matrix laboratory").
You can enter matrices into MATLAB in several different ways.
· Enter an explicit list of elements.
· Load matrices from external data files.
· Generate matrices using builtin functions.
· Create matrices with your own functions in Mfiles.
A matrix of one element is called a scalar and we can create one by entering:
>> s = 3
If you enter a variable this way, MATLAB immediately displays its value in the prompt window. To avoid this, put the ’;’ at the end of the line. You can use the ';' also to separate different commands on the same line like this
>> s = 3; a=10;
A matrix of one row or one column is called a vector. We can create a row vector by entering:
>> v = [8 3 5[
The result is a 1 by 3 matrix
To create a column vector we will use the ';' symbol to separate between the elements
>> u = [4.2; 5; 7]
The result is a 3 by 1 matrix
We can create a 3 by 3 matrix simply by entering
>> A
= [3 4 0; 1 4 6; 9 7 2]
We can now calculate
quantities like the transpose of A
>> A'
or the transpose of v which
will be a column vector:
>> v'
For example we can create another 3 by 3 matrix by entering
>> B=[v; u'; 1 0
4[
Note that if A has complex entries then A' means a complex conjugate transpose. A nonconjugate transpose is performed by
>> A.'
Check this by typing
[X, Y]=meshgrid(1:1);
Z=X+i*Y
Z.'Z
MATLAB has some builtin functions for the construction of special matrices.
>> ones(3,4)
will create a matrix of 1's with
3 rows and 4 columns. And you can figure out what
>> zeros(5,7(
does.
As you see, every Matlab
function (even those you will write) has input arguments and output
arguments . The input arguments are what you put between the brackets, for
example in ones(m,n) m and n are the inputs. The output is what you get back from
the function. For example in A = ones(5) the result of the ones
function is placed in A , which will be a 5by5
matrix of zeros. The variables you use to get the output are the glue between the
various commands. In order to find about the input and output of functions, simply
type the magic word help.
A command
>>rand(3,4)
creates a 3 by 4 matrix with random values chosen from a uniform distribution on the interval (0.0,1.0).
>>repmat(a,M,N)
creates an M by N matrix filled by values of a. It produces a large matrix consisting of M by N copies of a, if a is a matrix. To get an array of copies of a matrix A, try
>>B=repmat(A,[1 1 4])
and see what it does.
To create an identity matrix of order 5 you type
>> eye(5)
Another example matrix appears in the Renaissance
engraving Melancholia I by the German artist and amateur mathematician Albrecht
Dürer. This matrix is known as a magic square and has a property that it has equal row, column, and diagonal
sums. You can create such a matrix using a function magic
>> A=magic(4)
creates an 4by4 magic square.
If you type
>> x=
this will give a vector with all
the integer numbers from 1 to 10, and
>> y=0:0.1:8
will give a vector of all the numbers between 0 and 8 with a 0.1 step.
To summarize all we have learned, enter
>> s=3; y=[1 s ones(1,5); 5*rand(1,7); 1:7[;
and type y to see what you've got.
After we created a matrix we
might want to extract some elements from it. To extract the j'th element of a
vector you type
>> v(j)
Typing the command A(j,k) will extract the element of A that sits in the j'th row in the k'th column. For example:
>> A(2,3)
will extract the element in the 2^{nd} row in the 3^{rd} column, and
>> A(:,3:end(
will extract the last columns of
A beginning from the 3rd. The
first semicolon stands for all the indexes along the corresponding dimension
and end denotes
the last index in an indexing expression.
It is also possible to refer to
the elements of a matrix with a single subscript, A(k). This
is the usual way of referencing row and column vectors. But it can also apply
to a fully twodimensional matrix, in which case the array is regarded as one
long column vector formed from the columns of the original matrix. So, for our
magic square, A(8) is another way of referring to the value 15
stored in
A(4,2). If you try to use the value of an element outside of the matrix, it
is an error:
>>
t = A(4,5)
Index exceeds matrix dimensions.
On the other hand, if you store a value in an element outside of the matrix, the size increases to accommodate the newcomer:
>> X = A;
>> X(4,5) = 17
X =
16
3 2 13
0
5
10 11 8
0
9
6 7 12
0
4
15 14 1
17
If we have several matrices, we can add them
>> A+B
or multiply them
>> A*B
but we have to take care that the dimensions are appropriate: the number of columns in A should match the number of rows in B. The multiplication operation we just used obeys the rules of ordinary matrix multiplication. There is another kind of multiplication, which is multiplying each element in a matrix by the corresponding element of another matrix which has to be the same size (this is called array multiplication). To do this you enter
>> A.*B
This dot rule is general, so if you enter
>> A.^2
it will raise all the elements of A to the power of 2
(which is of
course different from A^2).
MATLAB uses a dot, or decimal point, as part of the notation for
elementbyelement operations.
The list of operators includes:
+ 

 

.* 

./ 

.\ 

.^ 

.' 

A command
>>sum(A)
performs a summation along columns. Suppose we want to know the sum of the 3^{rd} row. We will use the colon (:) and the function sum like this:
>>
sum(A(3,:))
for the sum of the third column we would type:
>>
sum(A(:,3))
The main diagonal of the matrix can be picked off by a function diag and the determinant is calculated by function det. The other diagonal, the socalled antidiagonal, is not so important mathematically, so MATLAB does not have a readymade function for it. But a function originally intended for use in graphics, fliplr flips a matrix from left to right.
>>sum(diag(fliplr(A)))
calculates the sum of the antidiagonal
Now, type the following command:
>>
A>3
What you got is a new matrix, with ones and zeros. Every element with one means that the corresponding element in A fulfills the condition "bigger than 3". Zeros mean that the opposite is true, that is, the number is smaller or equal to 3. We now go further by typing:
>> (A>3).* A
can you explain the result ?
hint: notice we used the element by element multiplier (.*)
Now try this and explain the result:
>> A.*A > 5*A
The logical vectors created from logical and relational operations
can be used to reference subarrays. Suppose X is an ordinary matrix and L is a
matrix of the same size that is the result of some logical operation. Then X(L)
specifies the elements of X
where the
elements of L are nonzero.
This kind of subscripting can be done in one step by specifying the logical operation as the subscripting expression. Suppose you have the following set of data.
x =
2.1 1.7 1.6 1.5
The
a questionnaire. To remove the missing data with logical indexing,
use finite(x), which is true for
all finite numerical values and false for
x = x(finite(x))
x =
2.1 1.7 1.6 1.5 1.9 1.8 1.5 5.1 1.8 1.4 2.2 1.6 1.8
Now there is one observation, 5.1, which seems to be very different from the others. It is an outlier. The following statement removes outliers, in this case those elements more than three standard deviations from the mean.
x = x(abs(xmean(x)) <=
3*std(x))
x =
2.1 1.7 1.6 1.5 1.9
1.8 1.5 1.8 1.4 2.2 1.6 1.8
The find
function determines the indices of array elements that meet a given
logical condition. In its simplest form, find
returns a column vector
of indices. Transpose that vector to obtain a row vector of indices. For
example
k = find(isprime(A))'
picks out the locations, using onedimensional indexing, of the primes in the magic square.
Display those primes, as a row vector in the order determined by k, with
A(k)
When you use k as a lefthandside index in an assignment statement, the matrix structure is preserved. For example,
A(k) =
will change previously prime entries of A into
While working with MATLAB, you can always list the variables in the workspace by entering who or whos ( whos shows also the size of the variables)
To log everything you type in the workspace (and the results), use the 'diary' command. To start, enter something like:
>> diary my_log.txt
and to stop saving, enter:
>> diary off
Type help diary to see exactly what it does.
You can store all your data currently defined in the workspace in a file :
>> save mydata
a file mydata.mat will be created with all the variables you currently have. The extension .mat indicates the MATLAB data format. Note, that if this file already exists, MATLAB will overwrite the old data. If you want to add some variable to an existing file mydata.mat type
>>save mydata A append
You can remove all the defined variables using the command clear and then load all the data again from file :
>> clear
>> load mydata
You can also load and store data using the MATLAB user interface. To save all variables in the workspace, choose a corresponding option in the File menu. If you want to save a certain variable, choose it in the Workspace window, click the right mouse button and select the appropriate option.
You can load a data file by selecting it in the Current Directory window. However, this does not allow you to load just selected variables from the file, which is possible using the command load:
>> load mydata A x
It is also possible to load the data stored in other formats rather than MATLAB. For this you need to Import data using the Import Wizard dialog. You can start a demo Importing Data from the helpdesk to get more information about this tool.
MATLAB uses columnoriented analysis for multivariate statistical data. It assumes that each column in a data set represents a variable and each row an observation. The (i,j)th element is the ith observation of the jth variable. For vector arguments, it does not matter whether the vectors are oriented in a row or column direction. For array arguments, the functions operate on columns. This means, for example, that if you apply max to an array, the result is a row vector containing the maximum values over each column.
Here we list few functions, which are very useful for data analysis:
max  maximum value
min  minimum value
mean  mean value
median  median value
std  standard deviation
sort – sorting
sum sum of the elements
prod  product of elements
cumsum  cumulative sum of elements
cumprod  cumulative product of elements
diff  differences and approximate derivatives
corrcoef  correlation coefficients
cov 
covariance matrix
Typing help on each of these functions will give you the full description of them.
Note that, although it is very easy to write most of these functions yourself,
it is nice to already have them and it makes the
scripts more compact.
As an example, consider a data set with three variables:
· Heart rate
· Weight
· Hours of exercise per week
For five observations, the resulting array might look like:
D =
72 134 3.2
81 201 3.5
69 156 7.1
82 148 2.4
75 170 1.2
To obtain the mean and standard deviation of each column:
>> mu = mean(D),
sigma = std(D)
For a list of the data analysis functions available in MATLAB, type
>> help datafun
If you have access to the Statistics Toolbox, type
>> help stats
To learn about the interpolation capabilities of MATLAB, try the following demo
>> census
Hypothesis testing for the difference in means of two samples.
Syntax
[h,significance,ci] = ttest2(x,y)
[h,significance,ci] =
ttest2(x,y,alpha)
[h,significance,ci] =
ttest2(x,y,alpha,tail)
Description
h = ttest2(x,y) performs a ttest to determine whether two samples from a normal distribution (in x and y) could have the same mean when the standard deviations are unknown but assumed equal.
h, the result, is 1 if you can reject the null hypothesis at the 0.05 significance level alpha and 0 otherwise.
significance is the pvalue associated with the Tstatistic.
significance is the probability that the observed value of T could be as large or larger under the null hypothesis that the mean of x is equal to the mean of y.
ci is a 95% confidence interval for the true difference in means.
[h,significance,ci]
= ttest2(x,y,alpha) gives control of the
significance level, alpha. For example if alpha = 0.01, and the result, h, is 1, you can reject the null
hypothesis at the significance
level 0.01. ci in this case is a
100(1alpha)% confidence interval for
the true difference in means.
ttest2(x,y,alpha,tail) allows specification of one or twotailed tests. tail is a flag that specifies one of three alternative hypotheses:
· tail = 0 (default) specifies the alternative, .
· tail = 1 specifies the alternative, .
· tail = 1 specifies the alternative, .
Examples
This example generates 100 normal random numbers with theoretical mean zero and standard deviation one. We then generate 100 more normal random numbers with theoretical mean one half and standard deviation one. The observed means and standard deviations are different from their theoretical values, of course. We test the hypothesis that there is no true difference between the two means. Notice that the true difference is only one half of the standard deviation of the individual observations, so we are trying to detect a signal that is only one half the size of the inherent noise in the process.
x = normrnd(0,1,100,1);
y = normrnd(0.5,1,100,1);
[h,significance,ci] = ttest2(x,y)
h =
1
significance =
0.0017
ci =
0.7352 0.1720
The result, h =1, means that we can reject the null hypothesis. The
significance is 0.0017, which means that by chance we would have
observed values of t more extreme than the one in this example in only 17
of 10,000 similar experiments! A 95% confidence interval on the mean is
[0.7352 0.1720], which includes the theoretical (and hypothesized) difference
of 0.5.
MATLAB is indeed a great tool for visualization of data of all sorts. The first command you should learn how to use is the plot. If
you want to plot vector y versus vector x just type
>> plot(x,y(
MATLAB will produce a figure window and will plot the graph.
For example, type the following:
>> x = 0:pi/8:pi*2;
>> m=sin(x);
>> plot(x, m)
To learn more about the plot command, try the help on it. To present several plots on the same figure, use the subplot command. To create another figure window, use the figure command.
Now suppose we have a normally distributed random variable y with the mean, which depends on the parameter x as sin(x), and has a constant standard deviation 0.1. Let's produce 100 'observations' of this variable per each x:
>>y=repmat(sin(x),100,1)+randn(100,length(x))*0.1;
Now, we plot the data:
>> plot(x,y,'.');
The 3rd argument specifies the symbol with which we plot the data. See help to find out what are the other options. Now, let's create a plot of the data with error bars:
>> errorbar(x,mean(y),std(y));
Now, we shall learn about a useful command: hist. With hist you get nice histograms very quickly. Try this:
>> x = rand(1000,1) ;
>> hist(x)
this will create a histogram of a random vector with 1000 elements.
>> max(x)
this is the greatest element
>> max(x) – min(x)
this is the difference between the greatest and lowest. A command range produces the same result. We can define the binning ourselves the following way:
>> dx=min(x):range(x)/100:max(x);
>> h=hist(x,dx);
>> bar(dx,h);
The command bar creates a bar graph like the one produced by the command hist.
The following commands will help you control the way your graphs
look: axis, xlabel, ylabel,
title, hold, grid, get, set, gca, gcf. You can also edit your graphs
manually and save figures in various formats using the
interactive menu of the window, in which your figure is displayed. An interactive
tool called Properties Editor
enables you to see and change object property values.
A 3D graph is plotted using commands surf and mesh. To illustrate this, consider the following example. The strength of connections between columns of neurons is often assumed to be dependent on a distance between these columns in form of a socalled Mexicanhat function, characterized by local excitation and longrange inhibition. Let us demonstrate such interactions of a neuron placed at the origin with its surrounding on an area of 3by3 mm
x=3:0.05:3;
[X, Y]=meshgrid(x);
R=X.^2+Y.^2;
surf(X,Y,R);
pause
surf(X,Y,exp(R/0.5)0.5*exp(R/3));
A command meshgrid returns the X and Y coordinates of the nodes of a grid built on the domain specified by vector x. The rows of the output array X and the columns of the output array Y are copies of the vector x. You can rotate the graph by clicking the corresponding button in the figure menu and then dragging the graph with the mouse. This way you can see the graph from any angle, especially you can observe any 2D projection.
The best way to learn about MATLAB visualization capabilities is to see the demo and then learn about the commands that are most useful for you. Note that the demos are written in MATLAB so you can "type" the files and just copy the commands you need.
Instead of entering all the commands in the command line you can write
them all in an Mfile. You create MATLAB scripts using a text editor, then use
them as you would any other MATLAB
command. Suppose you have created a script myscript.m , then you can simply run it from
the command prompt by typing
>> myscript
or use it in other scripts. For
this, you need to run both scripts from
the same directory or to specify the path to the directory where myscript.m is stored. Suppose you have stored your files in the directory C:\MATLAB6p5\work. To be able to access these files from other directories, type
>>addpath
C:\MATLAB6p5\work
or set the path to the directory using the File menu.
Scripts do not accept input arguments or
return output arguments. They operate on data in the workspace or can create
new data on which to operate. Any variables that they create remain in the
workspace, to be used in subsequent computations.
You can define your own functions which can accept input arguments and return output arguments.
Suppose, you have just discovered the formula that explains the universe and everything:
Now you can create a file called universe.m with the following content
function y=universe(x);
y=x.^2+42*sin(x);
Note, that the name of the function must be identical with the name of the file! The variable y defines the output variable of the function. This is what you get back from the function. Unlike in a script, the rest of the variables used by the function are internal and are not passed to the workspace from which the function was called. The input arguments are defined between the brackets. Your function can have several input arguments. For example, you can define an optional parameter a and set its default value to 42:
function y=universe1(x,a);
if nargin < 2
a=42;
end
y=x.^2+a*sin(x);
Inside the body of a userdefined
function, the variable nargin counts
the number of parameters actually used to call the function. Now, if you call
the function universe(x) it will not return an error but will calculate the value using a=42.
Another way to check if a certain variable is defined is to use a function exist:
if ~exist('a','var')
a=42;
end
Finally, you can declare a function with variable number of arguments using varargin.
function y=universe2(x,varargin)
varargin must be declared as the last input argument and collects all the inputs from that point onwards. Then we can use them the following way:
a=42;
b=1;
if nargin ==2
a=varargin{1};
elseif nargin>2
a=varargin{1};
b=varargin{2};
end
y=b*x.^2+a*sin(x)
Type: fplot(‘x^2+42*sin(x)’, [3*pi 3*pi]). As you guessed, this plots a function whose definition is specified by the first argument (a string). We could also enter 'universe' as a first argument, this would plot our function with its default values for a and b. The second argument is the x range of interest. Type help fplot to learn more.
Try to plot the function over larger and smaller ranges. Try to use the zoom tool to inspect the function in details. Add grids to the plot by : grid on (you can turn it of by grid off).
There is another way to define our function, which is especially suitable if we want to use the function as an argument of other functions of MATLAB:
f = inline('x^2 + 42*sin(x)')
You will get this answer:
f =
Inline function:
f(x) = x^2+42*sin(x)
Which means MATLAB now recognizes our function. In order to know how much the universe will last, you need to know for what x between –2 and –4 the function zeros. Type
t=fzero (f, [2 –4]))
If you want to test this, simply type: f(t). Use this function with different second arguments to find other intersections with zero. Now feel a strong urge to find the minimum points of f. This is simple.
Type: m=fmin(f,
8, 0)
This means that m will be assigned with the x that
gives a local minimum such that 8<m<0.
In order to perform analytical calculations of any kind, one can use a Symbolic Math Toolbox.
For example, we can calculate the derivative of f. Type:
s = diff('x^2 + 42*sin(x)');
The output will be a symbolic expression containing the derivative of the input, if the argument of the function is a string, or a symbolic expression (but not an inline function). You can define a symbolic variable by typing
>> syms x
Then we don't need to use '':
s = diff(x^2 + 42*sin(x));
If we want to plot the result, we need to use the string equivalent to s:
ss = char(s)
figure (this will create another figure)
fplot(ss, [3*pi 3*pi])
If we wanted to plot the integral f we would do
s = int(x^2 + 42*sin(x))
and repeat the other command as before.
Speed of simulations is usually considered as disadvantage of
MATLAB. It is especially slow if it has to operate with for  loops and allocate large amount
of memory. One of the possibilities to make programs work faster and still
enjoy the advantages of MATLAB interface is to write time and memoryconsuming
parts of the code in compiled languages like Fortran and C and use them as
functions in MATLAB scripts.
MATLAB gives this possibility via MEXfiles. MEXfiles
are dynamically linked subroutines produced from C or Fortran source code that
the MATLAB interpreter can automatically load and execute. They behave just
like Mfiles and builtin functions. While Mfiles have a platformindependent
extension, .m, MATLAB
identifies MEXfiles by platformspecific extensions. For Windows, for example,
it is .dll (dynamically
linked libraries).
Before you
can create MEXfiles on the Windows platform, you must define which compiler
you are going to use. For this, type
mex setup
You will get
a list of compilers installed on your platform.
MATLAB includes a C compiler, Lcc, that you can use to create C
MEXfiles if you don't have any other compiler installed.
Next, you need to create a Ccode.
The C code for a MEXfile differs from usual C programs, because it should be
able to convert MATLAB variables passed to the function into types
"digestable" by C, and transform the result of computations back into
a MATLAB array.
The source code for a MEXfile
consists of two distinct parts: A computational routine that contains
the code for performing the computations that you want implemented in the
MEXfile. Computations can be numerical computations as well as inputting and
outputting data. A gateway routine that interfaces the computational
routine with MATLAB by the entry point mexFunction (see example). Its
parameter right is an array of
input arguments, n_right is the number of input arguments, left is an array of output arguments,
and n_left
is the number of output arguments. The gateway calls the computational routine
as a subroutine.
Here is an example of
a function which acts like square root for positive inputs, but returns
negative values for negative inputs instead of imaginary. It is actually
equivalent to
sign(x).*sqrt(abs(x));
#include <stdlib.h>
#include <math.h>
#include
"mex.h"
void negsqrt(double *y, double *x, int n)
{
int i;
double dummy;
for
(i=0; i<n; ++i)
{
dummy
= x[i];
y[i]
= ((dummy < 0) ? sqrt(dummy) : sqrt(dummy));
}
return;
}
void mexFunction(int n_left, mxArray
*left[],
int
n_right, const mxArray *right[])
{
int n;
double
*source, *destination, dummy;
if
(n_left != 1)
mexErrMsgTxt("Please
enter the output array!\n");
if
(n_right != 1)
mexErrMsgTxt("Please
enter the input array!\n");
if
(!mxIsDouble(right[0]))
mexErrMsgTxt("Input
should be double!\n");
n
= mxGetNumberOfElements(right[0]);
left[0]
= mxCreateDoubleMatrix(1, n, mxREAL);
source
= mxGetPr(right[0]);
destination
= mxGetPr(left[0]);
/* do something... */
negsqrt(destination, source, n);
return;
}
1. Create a matrix of 100 random pairs in the range [1, 1]
2. Plot the numbers with red dots
3. Calculate the distances vector
4. Make a histogram of the distances vector
5. Find the maximal distance
6. Find the index of the maximal distance
7. Find the minimal distance
8. Find the index of the minimal distance
9. Plot the minimal distance and the maximal distance in blue (don’t forget to type hold on when adding new items to an existing figure)
1. Create a matrix with 10 rows and 4 columns. The numbers in the matrix should be random integers within 10 to 100.
2. Find the maximal values in each column.
3. Find the maximal values in each row.
4. Find the maximal entry in the matrix
5. Find the index of this entry (row and column)
6. Create a matrix with the same dimensions in which all entries equal the maximal value that you found in 4.
7. Calculate the difference between the matrix from 6 and the original matrix.
8. Find all the entries that are less than 40 and more than 20. Can you do it in one line?
1. Create a vector x with 100 random entries between 0 to 1
2. Create a vector y with 100 random entries between 0.1 to 1.1
3. Create a figure with two subplots. The upper one should show an histogram of x and the lower an histogram of y. Make sure that both histograms show the same range.
4. Calculate the mean, range and standard deviation for each vector.
5. Check the null hypothesis that the two vectors are from the same distribution (use ttest2 or ttest).
6. Find how many members of x are greater than their corresponding members in y.
7. Find the greatest difference between x and y.
1. Generate
spikes for 100 s using a Poisson spike generator with a constant rate of 100
Hz, and record their times of occurrence. Compute the coefficient of variation
of the interspike intervals, and the Fano factor for spike counts obtained over
counting intervals ranging from 1 to 100 ms. Plot the interspike interval
histogram.
Hints:
For a Poisson process with a firing rate r , the probability of an
interspike interval falling between τ and τ +t is
P[τ ≤ t_{i+1}
 t_{i} < τ +Δt] = rΔt exp(rτ). Compute spike times ti for
i = 1, 2,... iteratively by generating interspike intervals
from an exponential probability density. If xrand is uniformly
distributed over the range between zero and one, the negative of its logarithm
is exponentially distributed. Thus, we can generate spike times iteratively
from the
formula t_{i}_{+}_{1 }= t_{i}  ln(x_{rand})/r.
The coefficient of variation C_{V} is the ratio of the standard deviation to the
mean. The Fano factor for spike counts
is the ratio between the variance and mean of the spike count:
2. Add a
refractory period to the Poisson spike generator by allowing the firing rate to
depend on time. Initially, set the firing rate to a constant value r(t)=r_{0} . After every
spike, set r(t)=0 and then allow
it to recover exponentially back r_{0} with a time constant τ_{ref
}that controls the refractory recovery rate. In other words, have r(t) obey the equation
except
immediately after a spike, when it is set to 0. Plot the coefficient of
variation as a function of τ_{ref} over the range from 1ms
to 20ms and plot interspike interval histograms for a few different values of τ_{ref}.
Compute the Fano factor for spike counts obtained over counting intervals
ranging from 1 to 100 ms for the case τ_{ref} =10 ms.
Hints: The simplest way to simulate a differential
equation is according to Euler approximation. Divide the time into small
(compared to 1/r and to τ_{ref}) time intervals and calculate the rate iteratively:
However, for this simple equation, it is possible
to calculate the rate exactly, by integration from the time of the last spike t_{last}
:
The algorithm discussed in the exercise 1 only
works for constant firing rates. There are two ways to generate trains with a
variable firing rate. A simple procedure for generating spikes, which may be
more suitable for this exercise, is based on the fact that the estimated
probability of firing a spike during a short interval of duration Δt
is r(t)·t. The program progresses through time in small steps
of size Δt and
generates, at each time step, a random number x_{rand }chosen
uniformly in the range between zero and one. If r(t)·
Δt
> x_{rand}
at that time step, a spike is fired, otherwise it is not. For the second
method, see exercise 3.
3. Compute autocorrelation histograms of spike trains generated by Poisson generator with a constant firing rate
of 100 Hz, a constant firing rate of 100 Hz together with a refractory period
modeled as in exercise 2 with τ_{ref}=10ms, and a variable
firing rate r(t)=100(1+cos(2πt/50ms)). Plot the histograms
over a range from 0 to 100 ms. Approximate the firing rates from these spike
trains in the following ways:
a) by
binning time and counting spikes with Δt = 100 ms.
b) using a
sliding a rectangular window function along the spike train with Δt = 100 ms.
c) using a
Gaussian window function with σ_{t} = 100 ms.
d)
Approximate firing rate for an α function
window w(τ) = [α^{2}τ exp(ατ)]_{+}_{ }with 1/α = 100 ms. Hints:
The method used in exercise 1 can be extended to timedependent rates by using a procedure called the rejection sampling or spike thinning. The thinning technique requires a bound r_{max} on the estimated firing rate such that r ≤ r_{max} all times. We first generate a spike sequence corresponding to the constant rate r_{max} by iterating the rule t_{i}_{+}_{1 }= t_{i}  ln(x_{rand})/r. The spikes are then thinned by generating another y_{rand} for each i and removing the spike at time t_{i} from the train if r(t_{i})/r_{max} < y_{rand}. If r(t_{i})/r_{max }≥ y_{rand}, spike i is retained. Thinning corrects for the difference between the estimated timedependent rate and the maximum rate._{}
The spiketrain autocorrelation function is defined
as
, with
The
autocorrelation function is just a histogram of all the spike intervals that
are present in the data. It is a particular case of the crosscorelation
function. Here is a simple MATLAB function for crosscorrelation:
function crosscorr(t1,t2,N_bins);
M1 = length(t1);
M2 = length(t2);
D = ones(M2,1)*t1  t2'*ones(1,M1);
D = D(:);
hist(D,N_bins)
For
rate approximation, use function conv
which calculates a convolution of your spike train with a vector containing the
values of the window function on the appropriate time interval. The central
element of the vector should correspond the time delay τ=0.