dsCalcISI

CALCISI - Calculate the interspike interval.

 Usage:
   data = dsCalcISI(data,'option',value)

 Inputs:
   - data: DynaSim data structure (see dsCheckData)
   - options:
     'variable'         : name of field containing data on which to calculate
                          ISIs (default: *_spikes or first variable in data.labels)
     'threshold'        : scalar threshold value for detecting events (default: 0)
     'exclude_data_flag': whether to remove simulated data from result
                          structure (default: 0)
     'output_suffix'    : suffix to attach to output variable names (default: '')

 Outputs:
   - data: data structure with ISIs [ms] in .variable_ISI_SUA and .variable_ISI_MUA

 Notes:
 - "variable" can be specified as the name of a variable listed in
     data.labels, a cell array of string listing variable names, or as a regular
     expression pattern for identifying variables to process. See dsSelectVariables
     for more info on supported specifications.
 - DynaSim spike monitor returns spike data in variables *_spikes.
   - e.g., `data=dsSimulate('dv/dt=@current+10; {iNa,iK}; monitor v.spikes');`
     returns spikes in data.pop1_v_spikes (where 'pop1' is the default
     population name if not specified by the user).

 Examples:
   s.populations(1).name='E';
   s.populations(1).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65';
   s.populations(2).name='I';
   s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65';
   data=dsSimulate(s);
   data=dsCalcISI(data,'variable','*_v');
   data % contains ISIs for E and I pops in .E_v_ISI_SUA/MUA and .I_v_ISI_SUA/MUA

 See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables

0001 function data = dsCalcISI(data, varargin)
0002 %CALCISI - Calculate the interspike interval.
0003 %
0004 % Usage:
0005 %   data = dsCalcISI(data,'option',value)
0006 %
0007 % Inputs:
0008 %   - data: DynaSim data structure (see dsCheckData)
0009 %   - options:
0010 %     'variable'         : name of field containing data on which to calculate
0011 %                          ISIs (default: *_spikes or first variable in data.labels)
0012 %     'threshold'        : scalar threshold value for detecting events (default: 0)
0013 %     'exclude_data_flag': whether to remove simulated data from result
0014 %                          structure (default: 0)
0015 %     'output_suffix'    : suffix to attach to output variable names (default: '')
0016 %
0017 % Outputs:
0018 %   - data: data structure with ISIs [ms] in .variable_ISI_SUA and .variable_ISI_MUA
0019 %
0020 % Notes:
0021 % - "variable" can be specified as the name of a variable listed in
0022 %     data.labels, a cell array of string listing variable names, or as a regular
0023 %     expression pattern for identifying variables to process. See dsSelectVariables
0024 %     for more info on supported specifications.
0025 % - DynaSim spike monitor returns spike data in variables *_spikes.
0026 %   - e.g., `data=dsSimulate('dv/dt=@current+10; {iNa,iK}; monitor v.spikes');`
0027 %     returns spikes in data.pop1_v_spikes (where 'pop1' is the default
0028 %     population name if not specified by the user).
0029 %
0030 % Examples:
0031 %   s.populations(1).name='E';
0032 %   s.populations(1).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65';
0033 %   s.populations(2).name='I';
0034 %   s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65';
0035 %   data=dsSimulate(s);
0036 %   data=dsCalcISI(data,'variable','*_v');
0037 %   data % contains ISIs for E and I pops in .E_v_ISI_SUA/MUA and .I_v_ISI_SUA/MUA
0038 %
0039 % See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables
0040 
0041 %% 1.0 Check inputs
0042 options=dsCheckOptions(varargin,{...
0043   'variable',[],[],...
0044   'threshold',1e-5,[],... % slightly above zero in case variable is point process *_spikes {0,1}
0045   'exclude_data_flag',0,{0,1},...
0046   'output_suffix','',[],...
0047   'auto_gen_test_data_flag',0,{0,1},...
0048   },false);
0049 
0050 %% auto_gen_test_data_flag argin
0051 if options.auto_gen_test_data_flag
0052   varargs = varargin;
0053   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0054   varargs(end+1:end+2) = {'unit_test_flag',1};
0055   argin = [{data}, varargs]; % specific to this function
0056 end
0057 
0058 data = dsCheckData(data, varargin{:});
0059 % note: calling dsCheckData() at beginning enables analysis function to
0060 % accept data matrix [time x cells] in addition to DynaSim data structure.
0061 
0062 if numel(data)>1
0063   % use dsAnalyzeStudy to recursively call dsCalcISI on each data set
0064   data=dsAnalyzeStudy(data,@dsCalcISI,varargin{:});
0065   return;
0066 end
0067 
0068 % time info
0069 time = data.time;
0070 dt = time(2)-time(1);
0071 % ntime=length(time);
0072 
0073 % set defaults
0074 % default variable to process
0075 if isempty(options.variable)
0076   if any(~cellfun(@isempty,regexp(data.labels,'_spikes$')))
0077     % use results from DynaSim spike monitor
0078     options.variable=data.labels(~cellfun(@isempty,regexp(data.labels,'_spikes$')));
0079     if length(options.variable)==1 % store in string
0080       options.variable=options.variable{1};
0081     end
0082   else
0083     % use first state variable in model
0084 %     options.variable=data.labels{1};
0085   end
0086 end
0087 
0088 %% 2.0 set list of variables to process as cell array of strings
0089 options.variable=dsSelectVariables(data(1).labels,options.variable, varargin{:});
0090 
0091 %% 3.0 calculate ISIs for each variable
0092 if ~isfield(data,'results')
0093   data.results={};
0094 end
0095 
0096 % 3.2 loop over variables to process
0097 for v=1:length(options.variable)
0098   % extract this data set
0099   var=options.variable{v};
0100   dat=data.(var);
0101   % determine how many cells are in this data set
0102   ncells=size(dat,2);
0103   % loop over cells
0104   ISI_SUA=cell(1,ncells);
0105   spike_times=cell(1,ncells);
0106   for i=1:ncells
0107     % get spikes in this cell
0108     spike_inds=1+find((dat(2:end,i)>=options.threshold & dat(1:end-1,i)<options.threshold));
0109 %     spikes=zeros(ntime,1);
0110     spike_times{i}=time(spike_inds);
0111     if length(spike_inds)>1
0112 %       spikes(spike_inds)=1;
0113       % calculate ISIs
0114 %       ISI_SUA{i}=diff(find(spikes))*dt; %[ms]
0115       ISI_SUA{i}=diff(spike_times{i}); %[ms]
0116     end
0117   end
0118   
0119   ISI_MUA = vertcat(ISI_SUA{:});
0120     
0121   % add firing rates to data structure
0122   data.([var '_ISI_SUA' options.output_suffix])=ISI_SUA;
0123   data.([var '_ISI_MUA' options.output_suffix])=ISI_MUA;
0124   data.([var '_spike_times' options.output_suffix])=spike_times;
0125   if ~ismember([var '_ISI_SUA' options.output_suffix], data.results)
0126     data.results{end+1}=[var '_ISI_SUA' options.output_suffix];
0127   end
0128   if ~ismember([var '_ISI_MUA' options.output_suffix], data.results)
0129     data.results{end+1}=[var '_ISI_MUA' options.output_suffix];
0130   end
0131   if ~ismember([var '_spike_times' options.output_suffix], data.results)
0132     data.results{end+1}=[var '_spike_times' options.output_suffix];
0133   end
0134 end
0135 
0136 if options.exclude_data_flag
0137   for l=1:length(data.labels)
0138     data=rmfield(data,data.labels{l});
0139   end
0140 end
0141 
0142 %% auto_gen_test_data_flag argout
0143 if options.auto_gen_test_data_flag
0144   argout = {data}; % specific to this function
0145   
0146   dsUnitSaveAutoGenTestData(argin, argout);
0147 end