function [t,s] = fastgammatrain(duration,meanrate,order)

% function [t,s] = fastgammatrain(duration,meanrate,order)

%

% returns spike train output of an integrate and fire neuron with

% a random spiking threshold which is distributed according to

% a gamma distribution of order n. n=1 gives a poisson train.

% increasing n gives increasingly regular spike trains.

%

% where

% duration = duration of spiketrain (s)

% meanrate = mean firing rate of the spike train (spikes/s)

% order = order of the gamma distribution for resetting threshold

%

% return parameters are:

% t = vector containing the time index for spike train (s)

% s = spike train (1’s represent spikes, 0’s represent no spikes)

delta_t = 0.001; % resolution of spike train = 1 ms

t = 0:delta_t:duration;

idur=length(t);

Vrandth = gamrnd(order,1/order,1,idur);

itonextspike =round(Vrandth/(meanrate*delta_t));

itonextspike = itonextspike(itonextspike > 2); % refractory period = 2 ms

ispikes = cumsum(itonextspike);

ispikes = ispikes(ispikes < (duration/delta_t));

nspikes = length(ispikes);

s = sparse(1,ispikes,1,1,idur,nspikes);

if (nargout == 0)

plot(t,s)

end

s = full(s);