CALCACF - Calculate the autocorrelation function.
Usage:
data = dsCalcACF(data,'option',value)
Inputs:
- data: DynaSim data structure (see dsCheckData)
- options:
'variable' : name of field containing data on which to calculate
ACFs (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)
'muaSmearTimeWidth': gaussian convolution smearing in time of MUA spikes
for ACF (default: 10)
'output_suffix' : suffix to attach to output variable names (default: '')
Outputs:
- data: data structure with ACFs [ms] in .variable_ACF
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=dsCalcACF(data,'variable','*_v');
data % contains ACFs for E and I pops in .E_v_ACF and .I_v_ACF.
See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables0001 function data = dsCalcACF(data, varargin) 0002 %CALCACF - Calculate the autocorrelation function. 0003 % 0004 % Usage: 0005 % data = dsCalcACF(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 % ACFs (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 % 'muaSmearTimeWidth': gaussian convolution smearing in time of MUA spikes 0016 % for ACF (default: 10) 0017 % 'output_suffix' : suffix to attach to output variable names (default: '') 0018 % 0019 % Outputs: 0020 % - data: data structure with ACFs [ms] in .variable_ACF 0021 % 0022 % Notes: 0023 % - "variable" can be specified as the name of a variable listed in 0024 % data.labels, a cell array of string listing variable names, or as a 0025 % regular expression pattern for identifying variables to process. 0026 % See dsSelectVariables for more info on supported specifications. 0027 % 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=dsCalcACF(data,'variable','*_v'); 0040 % data % contains ACFs for E and I pops in .E_v_ACF and .I_v_ACF. 0041 % 0042 % See also: dsPlotFR, dsAnalyzeStudy, dsSimulate, dsCheckData, dsSelectVariables 0043 0044 %% 1.0 Check inputs 0045 options=dsCheckOptions(varargin,{... 0046 'variable',[],[],... 0047 'threshold',1e-5,[],... % slightly above zero in case variable is point process *_spikes {0,1} 0048 'exclude_data_flag',0,{0,1},... 0049 'numLags',1000,[],... 0050 'muaSmearTimeWidth',10,[],... 0051 'output_suffix','',[],... 0052 'auto_gen_test_data_flag',0,{0,1},... 0053 },false); 0054 0055 0056 %% auto_gen_test_data_flag argin 0057 if options.auto_gen_test_data_flag 0058 varargs = varargin; 0059 varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0; 0060 varargs(end+1:end+2) = {'unit_test_flag',1}; 0061 argin = [{data}, varargs]; % specific to this function 0062 end 0063 0064 data = dsCheckData(data, varargin{:}); 0065 % note: calling dsCheckData() at beginning enables analysis function to 0066 % accept data matrix [time x cells] in addition to DynaSim data structure. 0067 0068 if numel(data)>1 0069 % use dsAnalyzeStudy to recursively call dsCalcACF on each data set 0070 data=dsAnalyzeStudy(data,@dsCalcACF,varargin{:}); 0071 return; 0072 end 0073 0074 % time info 0075 time = data.time; 0076 dt = time(2)-time(1); 0077 ntime=length(time); 0078 0079 % set defaults 0080 % default variable to process 0081 if isempty(options.variable) 0082 if any(~cellfun(@isempty,regexp(data.labels,'_spikes$'))) 0083 % use results from DynaSim spike monitor 0084 options.variable=data.labels(~cellfun(@isempty,regexp(data.labels,'_spikes$'))); 0085 if length(options.variable)==1 % store in string 0086 options.variable=options.variable{1}; 0087 end 0088 else 0089 % use first state variable in model 0090 % options.variable=data.labels{1}; 0091 end 0092 end 0093 0094 numLags = min(options.numLags, ntime-1); 0095 0096 %% 2.0 set list of variables to process as cell array of strings 0097 options.variable=dsSelectVariables(data(1).labels,options.variable, varargin{:}); 0098 0099 %% 3.0 calculate ACFs for each variable 0100 if ~isfield(data,'results') 0101 data.results={}; 0102 end 0103 0104 % 3.2 loop over variables to process 0105 for v=1:length(options.variable) 0106 % extract this data set 0107 var=options.variable{v}; 0108 dat=data.(var); 0109 % determine how many cells are in this data set 0110 ncells=size(dat,2); 0111 % loop over cells 0112 ACF_SUA = zeros(numLags+1, size(dat,2)); 0113 spikes_MUA=zeros(ntime,1); 0114 spike_times = cell(1,ncells); 0115 for i=1:ncells 0116 % get spikes in this cell 0117 spike_inds=1+find((dat(2:end,i)>=options.threshold & dat(1:end-1,i)<options.threshold)); 0118 spikes=zeros(ntime,1); 0119 spike_times{i}=time(spike_inds); 0120 if any(spike_inds) 0121 spikes(spike_inds)=1; 0122 % calculate ACFs 0123 ACF_SUA(:,i)= autocorr(spikes, numLags); 0124 end 0125 spikes_MUA = spikes_MUA + spikes; 0126 end 0127 0128 ACF_SUA(1,:) = []; %remove ACF=1 at 0 lag 0129 0130 if options.muaSmearTimeWidth ~= 0 0131 gausWin = gausswin(options.muaSmearTimeWidth/dt); 0132 gausWin = gausWin/trapz(gausWin); %0 area 0133 end 0134 0135 spikes_MUA = conv(spikes_MUA, gausWin,'same'); 0136 ACF_MUA = autocorr(spikes_MUA, numLags); 0137 ACF_MUA(1,:) = []; %remove ACF=1 at 0 lag 0138 0139 % ACF_MUA = mean(ACF_SUA,2); 0140 0141 % add firing rates to data structure 0142 data.([var '_ACF_SUA' options.output_suffix])=ACF_SUA; 0143 data.([var '_ACF_MUA' options.output_suffix])=ACF_MUA; 0144 data.([var '_spike_times' options.output_suffix])=spike_times; 0145 if ~ismember([var '_ACF_SUA' options.output_suffix],data.results) 0146 data.results{end+1}=[var '_ACF_SUA' options.output_suffix]; 0147 end 0148 if ~ismember([var '_ACF_MUA' options.output_suffix], data.results) 0149 data.results{end+1}=[var '_ACF_MUA' options.output_suffix]; 0150 end 0151 if ~ismember([var '_spike_times' options.output_suffix], data.results) 0152 data.results{end+1}=[var '_spike_times' options.output_suffix]; 0153 end 0154 end 0155 0156 if options.exclude_data_flag 0157 for l=1:length(data.labels) 0158 data=rmfield(data,data.labels{l}); 0159 end 0160 end 0161 0162 %% auto_gen_test_data_flag argout 0163 if options.auto_gen_test_data_flag 0164 argout = {data, modifications}; % specific to this function 0165 0166 dsUnitSaveAutoGenTestData(argin, argout); 0167 end