CALCFR - Calculate firing rage for DynaSim data
Usage:
data = dsCalcFR(data,'option',value)
Inputs:
- data: DynaSim data structure (see dsCheckData)
- options:
'variable' : name of field containing data on which to calculate
firing rates (default: *_spikes or first variable in data.labels)
'threshold' : scalar threshold value for detecting events (default: 0)
'bin_size' : size of temporal window over which to calculate rate
[ms or fraction of data set] (default: 5% of the data set)
'bin_shift' : how much to shift the bin before calculating rate
again [ms or fraction of data set] (default: 1% of the data set)
'exclude_data_flag': whether to remove simulated data from result
structure (default: 0)
Outputs:
- data: data structure with firing rates [Hz] in .variable_FR
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=dsCalcFR(data,'variable','*_v');
data % contains firing rates for E and I pops in .E_v_FR and .I_v_FR.
See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables
Author: Jason Sherfey, PhD <jssherfey@gmail.com>
Copyright (C) 2016 Jason Sherfey, Boston University, USA0001 function data = dsCalcFR(data, varargin) 0002 %CALCFR - Calculate firing rage for DynaSim data 0003 % 0004 % Usage: 0005 % data = dsCalcFR(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 % firing rates (default: *_spikes or first variable in data.labels) 0012 % 'threshold' : scalar threshold value for detecting events (default: 0) 0013 % 'bin_size' : size of temporal window over which to calculate rate 0014 % [ms or fraction of data set] (default: 5% of the data set) 0015 % 'bin_shift' : how much to shift the bin before calculating rate 0016 % again [ms or fraction of data set] (default: 1% of the data set) 0017 % 'exclude_data_flag': whether to remove simulated data from result 0018 % structure (default: 0) 0019 % 0020 % Outputs: 0021 % - data: data structure with firing rates [Hz] in .variable_FR 0022 % 0023 % Notes: 0024 % - "variable" can be specified as the name of a variable listed in 0025 % data.labels, a cell array of string listing variable names, or as a regular 0026 % expression pattern for identifying variables to process. See dsSelectVariables 0027 % for more info on supported specifications. 0028 % - DynaSim spike monitor returns spike data in variables *_spikes. 0029 % - e.g., `data=dsSimulate('dv/dt=@current+10; {iNa,iK}; monitor v.spikes');` 0030 % returns spikes in data.pop1_v_spikes (where 'pop1' is the default 0031 % population name if not specified by the user). 0032 % 0033 % Examples: 0034 % s.populations(1).name='E'; 0035 % s.populations(1).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65'; 0036 % s.populations(2).name='I'; 0037 % s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65'; 0038 % data=dsSimulate(s); 0039 % data=dsCalcFR(data,'variable','*_v'); 0040 % data % contains firing rates for E and I pops in .E_v_FR and .I_v_FR. 0041 % 0042 % See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables 0043 % 0044 % Author: Jason Sherfey, PhD <jssherfey@gmail.com> 0045 % Copyright (C) 2016 Jason Sherfey, Boston University, USA 0046 0047 %% 1.0 Check inputs 0048 options=dsCheckOptions(varargin,{... 0049 'variable',[],[],... 0050 'time_limits',[-inf inf],[],... 0051 'threshold',1e-5,[],... % slightly above zero in case variable is point process *_spikes {0,1} 0052 'bin_size',.05,[],... 0053 'bin_shift',.01,[],... 0054 'exclude_data_flag',0,{0,1},... 0055 'output_suffix','',[],... 0056 'auto_gen_test_data_flag',0,{0,1},... 0057 },false); 0058 0059 %% auto_gen_test_data_flag argin 0060 if options.auto_gen_test_data_flag 0061 varargs = varargin; 0062 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0063 varargs(end+1:end+2) = {'unit_test_flag',1}; 0064 argin = [{data}, varargs]; % specific to this function 0065 end 0066 0067 data = dsCheckData(data, varargin{:}); 0068 % note: calling dsCheckData() at beginning enables analysis function to 0069 % accept data matrix [time x cells] in addition to DynaSim data structure. 0070 0071 if numel(data)>1 0072 % use dsAnalyzeStudy to recursively call dsCalcFR on each data set 0073 data=dsAnalyzeStudy(data,@dsCalcFR,varargin{:}); 0074 return; 0075 end 0076 0077 % time info 0078 time = data.time; 0079 dt = time(2)-time(1); 0080 t1=nearest(time,options.time_limits(1)); % index to first sample 0081 t2=nearest(time,options.time_limits(2)); % index to last sample 0082 ntime=t2-t1+1; 0083 0084 % set defaults 0085 % default variable to process 0086 if isempty(options.variable) 0087 if any(~cellfun(@isempty,regexp(data.labels,'_spikes$'))) 0088 % use results from DynaSim spike monitor 0089 options.variable=data.labels(~cellfun(@isempty,regexp(data.labels,'_spikes$'))); 0090 if length(options.variable)==1 % store in string 0091 options.variable=options.variable{1}; 0092 end 0093 else 0094 % use first state variable in model 0095 %options.variable=data.labels{1}; 0096 end 0097 end 0098 0099 % check bin_size 0100 if options.bin_size>1 0101 % convert from ms to time points 0102 options.bin_size=ceil(options.bin_size/dt); 0103 else 0104 % convert from fraction to time points 0105 options.bin_size=ceil(options.bin_size*ntime); 0106 end 0107 0108 % constrain bin_size to entire data set 0109 if options.bin_size>ntime 0110 options.bin_size=ntime; 0111 end 0112 0113 % check bin_shift 0114 if options.bin_shift>1 0115 % convert from ms to time points 0116 options.bin_shift=ceil(options.bin_shift/dt); 0117 else 0118 % convert from fraction to time points 0119 options.bin_shift=ceil(options.bin_shift*ntime); 0120 end 0121 0122 %% 2.0 set list of variables to process as cell array of strings 0123 options.variable=dsSelectVariables(data(1).labels,options.variable, varargin{:}); 0124 0125 %% 3.0 calculate firing rates for each variable 0126 if ~isfield(data,'results') 0127 data.results={}; 0128 end 0129 0130 % 3.1 calc bin info 0131 % samples at which bins begin 0132 bin_index_begs=t1:options.bin_shift:t2; 0133 % samples at which bins end 0134 bin_index_ends=bin_index_begs+options.bin_size; 0135 0136 if bin_index_ends(end)>t2 0137 if length(bin_index_ends) > 1 %multiple bins 0138 % remove final bin if extends beyond data 0139 bin_index_begs=bin_index_begs(bin_index_ends<=t2); 0140 bin_index_ends=bin_index_ends(bin_index_ends<=t2); 0141 else %1 bin 0142 bin_index_ends = t2; 0143 end 0144 end 0145 0146 % times at which bins begin 0147 bin_times=time(bin_index_begs); 0148 % number of bins 0149 nbins=length(bin_index_begs); 0150 % time width of a single bin in seconds 0151 bin_width=(dt/1000)*options.bin_size; 0152 0153 % 3.2 loop over variables to process 0154 for v=1:length(options.variable) 0155 % extract this data set 0156 var=options.variable{v}; 0157 dat=data.(var); 0158 % determine how many cells are in this data set 0159 ncells=size(dat,2); 0160 % loop over cells 0161 FR=zeros(nbins,ncells); 0162 spike_times=cell(1,ncells); 0163 for i=1:ncells 0164 % get spikes in this cell 0165 spike_inds=1+find((dat(2:end,i)>=options.threshold & dat(1:end-1,i)<options.threshold)); 0166 spikes=zeros(ntime,1); 0167 spike_times{i}=time(spike_inds); 0168 if any(spike_inds) 0169 spikes(spike_inds)=1; 0170 % calculate firing rates 0171 for bin=1:nbins 0172 % (# spikes in bin) / (duration of bin in seconds) 0173 FR(bin,i)=sum(spikes(bin_index_begs(bin):bin_index_ends(bin)))/bin_width; 0174 end 0175 end 0176 end 0177 % add firing rates to data structure 0178 data.([var '_FR' options.output_suffix])=FR; 0179 data.([var '_spike_times' options.output_suffix])=spike_times; 0180 if ~ismember([var '_FR'],data.results) 0181 data.results{end+1}=[var '_FR' options.output_suffix]; 0182 data.results{end+1}=[var '_spike_times' options.output_suffix]; 0183 end 0184 end 0185 % add bin times to data 0186 data.(['time_FR' options.output_suffix])=bin_times; 0187 if ~ismember(['time_FR' options.output_suffix],data.results) 0188 data.results{end+1}=['time_FR' options.output_suffix]; 0189 end 0190 if options.exclude_data_flag 0191 for l=1:length(data.labels) 0192 data=rmfield(data,data.labels{l}); 0193 end 0194 end 0195 0196 %% auto_gen_test_data_flag argout 0197 if options.auto_gen_test_data_flag 0198 argout = {data}; % specific to this function 0199 0200 dsUnitSaveAutoGenTestData(argin, argout); 0201 end