GenerateSpikeTrain

GenerateSpikeTrain - Generate spike train using a leaky integrate and fire neuron model.

  USAGE

    [spikes,V] = GenerateSpikeTrain(I,threshold,<options>)

    I              input current <a href="matlab:help samples">samples</a>
    threshold      spiking threshold
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'delay'       delay in seconds (default = 0)
     'tau'         decay rate (default = 0.015)
     'Vr'          resting membrane potential (default = 0)
     'sigma'       noise standard deviation (default = 0, no noise)
     'R'           input resistance
    =========================================================================

  OUTPUT

    spikes         list of spike timestamps
    V              membrane potentials across time


  EXAMPLE

    To get spike trains of two neurons receiving correlated inputs, provide
      I = mu + sigma(sqrt(1-c)*E + sqrt(c)*C)
    where c is the correlation, E is the exclusive input (or noise) and C is
    the common signal.

0001 function [spikes,V] = GenerateSpikeTrain(I,threshold,varargin);
0002 
0003 %GenerateSpikeTrain - Generate spike train using a leaky integrate and fire neuron model.
0004 %
0005 %  USAGE
0006 %
0007 %    [spikes,V] = GenerateSpikeTrain(I,threshold,<options>)
0008 %
0009 %    I              input current <a href="matlab:help samples">samples</a>
0010 %    threshold      spiking threshold
0011 %    <options>      optional list of property-value pairs (see table below)
0012 %
0013 %    =========================================================================
0014 %     Properties    Values
0015 %    -------------------------------------------------------------------------
0016 %     'delay'       delay in seconds (default = 0)
0017 %     'tau'         decay rate (default = 0.015)
0018 %     'Vr'          resting membrane potential (default = 0)
0019 %     'sigma'       noise standard deviation (default = 0, no noise)
0020 %     'R'           input resistance
0021 %    =========================================================================
0022 %
0023 %  OUTPUT
0024 %
0025 %    spikes         list of spike timestamps
0026 %    V              membrane potentials across time
0027 %
0028 %
0029 %  EXAMPLE
0030 %
0031 %    To get spike trains of two neurons receiving correlated inputs, provide
0032 %      I = mu + sigma(sqrt(1-c)*E + sqrt(c)*C)
0033 %    where c is the correlation, E is the exclusive input (or noise) and C is
0034 %    the common signal.
0035 
0036 
0037 % Copyright (C) 2015-2016 by Ralitsa Todorova, Michaël Zugaro
0038 %
0039 % This program is free software; you can redistribute it and/or modify
0040 % it under the terms of the GNU General Public License as published by
0041 % the Free Software Foundation; either version 3 of the License, or
0042 % (at your option) any later version.
0043 
0044 
0045 % Defaults
0046 delay = 0;
0047 Vr = 0;
0048 tau = 0.015;
0049 sigma = 0;
0050 R = 1;
0051 
0052 % Check parameters
0053 if ~isdmatrix(I,'@2'),
0054     error('Incorrect input current (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0055 end
0056 if ~isdscalar(threshold),
0057     error('Incorrect spike threshold (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0058 end
0059 
0060 % Browse parameters
0061 for i = 1:2:length(varargin),
0062     if ~ischar(varargin{i}),
0063         error(['Optional parameter ' num2str(i) ' is not a valid property (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).']);
0064     end
0065     switch(lower(varargin{i})),
0066         case 'delay',
0067             delay = varargin{i+1};
0068             if ~isdscalar(delay,'>=0'),
0069                 error('Incorrect value for property ''delay'' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0070             end
0071         case 'tau',
0072             tau = varargin{i+1};
0073             if ~isdscalar(tau,'>=0'),
0074                 error('Incorrect value for property ''tau'' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0075             end
0076         case 'vr',
0077             Vr = varargin{i+1};
0078             if ~isdscalar(Vr),
0079                 error('Incorrect value for property ''Vr'' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0080             end
0081         case 'sigma',
0082             sigma = varargin{i+1};
0083             if ~isdscalar(sigma,'>=0'),
0084                 error('Incorrect value for property ''sigma'' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0085             end
0086         case 'r',
0087             R = varargin{i+1};
0088             if ~isdscalar(R),
0089                 error('Incorrect value for property ''R'' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).');
0090             end
0091         otherwise,
0092             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help GenerateSpikeTrain">GenerateSpikeTrain</a>'' for details).']);
0093     end
0094 end
0095 
0096 % Initialization
0097 Vmin = -threshold;
0098 if tau == 0,
0099     warning('The value of tau cannot be exactly 0 (division by zero) - using ''realmin'' instead');
0100     tau = realmin;
0101 end
0102 dt = median(diff(I(:,1)));
0103 m = min(I(:,1));
0104 M = max(I(:,1));
0105 t = (m:dt:M+dt/100)'; % +dt/100 added for underflow errors
0106 delay = round(delay/dt);
0107 
0108 % Create Gaussian noise
0109 randn('seed',sum(100*clock)); % set seed for random number generator
0110 n = length(t);
0111 noise = sigma*sqrt(1/dt)*randn(n,1);
0112 
0113 % Make sure I is evenly sampled (interpolate)
0114 I = Interpolate(I,t);
0115 
0116 % Start
0117 spikes = [];
0118 V = nan(n,1);
0119 V(1) = Vr;
0120 
0121 for i = 1:n-1,
0122     if V(i) > threshold,
0123         spikes = [spikes;t(i)];
0124         V(i) = Vr;
0125     elseif V(i) < Vmin,
0126         V(i) = Vmin;
0127     end
0128     if i+delay > length(I),
0129         V(i+1) = V(i) + dt*(-(V(i) - Vr)/tau);
0130     else
0131         V(i+1) = V(i) + dt*((R*I(i+delay,2)-(V(i) - Vr))/tau + noise(i+delay)/sqrt(tau));
0132     end
0133 end