CALCFRMULTI - extends dsCalcFR to get SUA and MUA firing rates
Usage:
data = dsCalcFRmulti(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)
'time_limits' : [beg,end] (units of data.time)
'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)
'output_suffix' : suffix to attach to output variable names (default: '')
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_SUA/MUA and .I_v_FR_SUA/MUA.
See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables0001 function data = dsCalcFRmulti(data, varargin) 0002 %CALCFRMULTI - extends dsCalcFR to get SUA and MUA firing rates 0003 % 0004 % Usage: 0005 % data = dsCalcFRmulti(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 % 'time_limits' : [beg,end] (units of data.time) 0013 % 'threshold' : scalar threshold value for detecting events (default: 0) 0014 % 'bin_size' : size of temporal window over which to calculate rate 0015 % [ms or fraction of data set] (default: 5% of the data set) 0016 % 'bin_shift' : how much to shift the bin before calculating rate again [ms 0017 % or fraction of data set] (default: 1% of the data set) 0018 % 'exclude_data_flag': whether to remove simulated data from result 0019 % structure (default: 0) 0020 % 'output_suffix' : suffix to attach to output variable names (default: '') 0021 % 0022 % Outputs: 0023 % - data: data structure with firing rates [Hz] in .variable_FR 0024 % 0025 % Notes: 0026 % - "variable" can be specified as the name of a variable listed in 0027 % data.labels, a cell array of string listing variable names, or as a 0028 % regular expression pattern for identifying variables to process. 0029 % See dsSelectVariables for more info on supported specifications. 0030 % - DynaSim spike monitor returns spike data in variables *_spikes. 0031 % - e.g., `data=dsSimulate('dv/dt=@current+10; {iNa,iK}; monitor v.spikes');` 0032 % returns spikes in data.pop1_v_spikes (where 'pop1' is the default 0033 % population name if not specified by the user). 0034 % 0035 % Examples: 0036 % s.populations(1).name='E'; 0037 % s.populations(1).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65'; 0038 % s.populations(2).name='I'; 0039 % s.populations(2).equations='dv/dt=@current+10; {iNa,iK}; v(0)=-65'; 0040 % data=dsSimulate(s); 0041 % data=dsCalcFR(data,'variable','*_v'); 0042 % data % contains firing rates for E and I pops in .E_v_FR_SUA/MUA and .I_v_FR_SUA/MUA. 0043 % 0044 % See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables 0045 0046 %% 1.0 Check inputs 0047 options=dsCheckOptions(varargin,{... 0048 'variable',[],[],... 0049 'time_limits',[-inf inf],[],... 0050 'threshold',1e-5,[],... % slightly above zero in case variable is point process *_spikes {0,1} 0051 'bin_size',.05,[],... % 30 0052 'bin_shift',.01,[],... % 10 0053 'exclude_data_flag',0,{0,1},... 0054 'output_suffix','',[],... 0055 'auto_gen_test_data_flag',0,{0,1},... 0056 },false); 0057 0058 %% auto_gen_test_data_flag argin 0059 if options.auto_gen_test_data_flag 0060 varargs = varargin; 0061 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0062 varargs(end+1:end+2) = {'unit_test_flag',1}; 0063 argin = [{data}, varargs]; % specific to this function 0064 end 0065 0066 data = dsCheckData(data, varargin{:}); 0067 % note: calling dsCheckData() at beginning enables analysis function to 0068 % accept data matrix [time x cells] in addition to DynaSim data structure. 0069 0070 if numel(data)>1 0071 % use dsAnalyzeStudy to recursively call dsCalcFRmulti on each data set 0072 data=dsAnalyzeStudy(data,@dsCalcFRmulti,varargin{:}); 0073 return; 0074 end 0075 0076 % time info 0077 time = data.time; 0078 dt = time(2)-time(1); 0079 ntime=length(time); 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_SUA=zeros(nbins,ncells); 0162 FR_MUA=zeros(nbins,1); 0163 % spike_times=cell(1,ncells); 0164 % spike_inds = zeros(ntime, ncells); 0165 0166 % for i=1:ncells 0167 % % get spikes in this cell 0168 % spike_inds=1+find((dat(2:end,i)>=options.threshold & dat(1:end-1,i)<options.threshold)); 0169 % spikes=zeros(ntime,1); 0170 % spike_times{i}=time(spike_inds); 0171 % end 0172 0173 spikes = [zeros(1,ncells); (dat(2:end,:)>=options.threshold & dat(1:end-1,:)<options.threshold)]; 0174 0175 if any(spikes(:)) 0176 % calculate firing rates 0177 for bin=1:nbins 0178 % (# spikes in bin) / (duration of bin in seconds) 0179 FR_SUA(bin,:)=sum(spikes(bin_index_begs(bin):bin_index_ends(bin), :))/bin_width; 0180 FR_MUA(bin)=sum(sum(spikes(bin_index_begs(bin):bin_index_ends(bin), :)))/bin_width; 0181 end 0182 end 0183 0184 % add firing rates to data structure 0185 data.([var '_FR_SUA' options.output_suffix])=FR_SUA; 0186 data.([var '_FR_MUA' options.output_suffix])=FR_MUA; 0187 % data.([var '_spike_times'])=spike_times; 0188 if ~ismember([var '_FR_SUA' options.output_suffix], data.results) 0189 data.results{end+1}=[var '_FR_SUA' options.output_suffix]; 0190 % data.results{end+1}=[var '_spike_times']; 0191 end 0192 if ~ismember([var '_FR_MUA' options.output_suffix], data.results) 0193 data.results{end+1}=[var '_FR_MUA' options.output_suffix]; 0194 end 0195 end 0196 % add bin times to data 0197 data.(['time_FR' options.output_suffix])=bin_times; 0198 if ~ismember(['time_FR' options.output_suffix], data.results) 0199 data.results{end+1}=['time_FR' options.output_suffix]; 0200 end 0201 if options.exclude_data_flag 0202 for l=1:length(data.labels) 0203 data=rmfield(data,data.labels{l}); 0204 end 0205 end 0206 0207 %% auto_gen_test_data_flag argout 0208 if options.auto_gen_test_data_flag 0209 argout = {data}; % specific to this function 0210 0211 dsUnitSaveAutoGenTestData(argin, argout); 0212 end