Description of dsCalcSpikeSync

0001 function stats = dsCalcSpikeSync(data, varargin)
0002 %CALCSPIKESYNC - Compute spike synchronization between spiketrains
0003 %
0004 % Usage:
0005 %   stats = dsCalcSpikeSync(data,'option',value)
0006 %
0007 % Inputs:
0008 %   - data: DynaSim data structure (see dsCheckData)
0009 %   - options:
0010 %     'ROI_pairs'   : {'var1',roi1,'var2',roi2; ...}
0011 %     'kernel_width': ms, width of gaussian for kernel regression (default: 1)
0012 %     'Ts'          : ms, set to this effective time step for rate process
0013 %                     before regression (default: 1)
0014 %     'maxlag_time' : ms, max lag time for cross correlation (default: 10)
0015 %     'spike_threshold',0,[],...:  threshold for spike detection
0016 %       - Note: Fractional roi=[a b] selects the interval (a,b]
0017 %         - Example: N=10: [0 .5] -> [1:5], [.5 1]->[6:10]
0018 %
0019 % Examples:
0020 %   ROI_pairs={'E_v',[0 1],'E_v',[0 1]; 'E_v',[0 1],'I_v',[0 1]};
0021 %   ROI_pairs={'E_v',[0 .49],'E_v',[.5 1]};
0022 %   ROI_pairs={'E_v',1:4,'E_v',4:8};
0023 %
0024 %   spike_threshold=0; % same for all ROIs
0025 %   spike_threshold=[0 .5]; % use 0 for all ROI1s and .5 for all ROI2s
0026 %   spike_threshold=[0 .5; 0 .25];
0027 %
0028 % Author: Jason Sherfey, PhD <jssherfey@gmail.com>
0029 % Copyright (C) 2016 Jason Sherfey, Boston University, USA
0030 
0031 %% 1.0 Check inputs
0032 options=dsCheckOptions(varargin,{...
0033   'ROI_pairs',[],[],...
0034   'spike_threshold',0,[],... % threshold for spike detection
0035   'kernel_width',1,[],... % ms, width of gaussian for kernel regression
0036   'Ts',1,[],...           % ms, set to this effective time step for rate process before regression
0037   'maxlag_time',10,[],... % ms, max lag time for cross correlation
0038   'time_limits',[100 inf],[],... % time limits for spectral analysis
0039   'auto_gen_test_data_flag',0,{0,1},...
0040   },false);
0041 
0042 %% auto_gen_test_data_flag argin
0043 if options.auto_gen_test_data_flag
0044   varargs = varargin;
0045   varargs{find(strcmp(varargs, 'auto_gen_test_data_flag'))+1} = 0;
0046   varargs(end+1:end+2) = {'unit_test_flag',1};
0047   argin = [{data}, varargs]; % specific to this function
0048 end
0049 
0050 data=dsCheckData(data, varargin{:});
0051 
0052 if numel(data)>1
0053   error('dsCalcSpikeSync currently only supports one data set at a time');
0054 end
0055 if isempty(options.ROI_pairs)
0056   % set default ROIs
0057   var=data.labels{1};
0058   roi=[0 1];
0059   options.ROI_pairs={var,roi,var,roi};
0060 end
0061 if size(options.ROI_pairs,2)<4
0062   error('ROIs must include four elements: {var1,roi1,var2,roi2}');
0063 end
0064 
0065 % determine ROI pairs for which to calculate spike synchrony
0066 pop_names={data.model.specification.populations.name};
0067 pop_sizes=[data.model.specification.populations.size];
0068 npairs=size(options.ROI_pairs,1);
0069 VAR1=cell(npairs,1); ROI1=VAR1;
0070 VAR2=cell(npairs,1); ROI2=VAR2;
0071 exclude=[];
0072 for i=1:npairs
0073   VAR1{i}=options.ROI_pairs{i,1};
0074   VAR2{i}=options.ROI_pairs{i,3};
0075   roi1=options.ROI_pairs{i,2};
0076   roi2=options.ROI_pairs{i,4};
0077   % set ROI in cell indices for var1
0078   name=regexp(VAR1{i},'^([a-zA-Z0-9]+)_','tokens','once');
0079   if ~ismember(name{1},pop_names)
0080     exclude=[exclude i];
0081     continue;
0082   end
0083   N=pop_sizes(strcmp(pop_names,name{1}));
0084   if numel(roi1)==2 && any(roi1<1)
0085     x=floor(N*roi1); a=x(1)+1; b=x(2);
0086     %a=round(roi1(1)*N);%round(1+roi1(1)*(N-1));
0087     %b=round(roi1(2)*N);%round(1+roi1(2)*(N-1));
0088   else
0089     a=roi1(1);
0090     b=roi1(2);
0091   end
0092   ROI1{i}=a:b;
0093   % set ROI in cell indices for var2
0094   name=regexp(VAR2{i},'^([a-zA-Z0-9]+)_','tokens','once');
0095   if ~ismember(name{1},pop_names)
0096     exclude=[exclude i];
0097     continue;
0098   end
0099   N=pop_sizes(strcmp(pop_names,name{1}));
0100   if numel(roi2)==2 && any(roi2<1)
0101     x=floor(N*roi2); a=x(1)+1; b=x(2);
0102     %a=round(roi2(1)*N);%round(1+roi2(1)*(N-1));
0103     %b=round(roi2(2)*N);%round(1+roi2(2)*(N-1));
0104   else
0105     a=roi2(1);
0106     b=roi2(2);
0107   end
0108   ROI2{i}=a:b;
0109 end
0110 if ~isempty(exclude)
0111   options.ROI_pairs(exclude,:)=[];
0112   VAR1(exclude)=[];
0113   VAR2(exclude)=[];
0114   ROI1(exclude)=[];
0115   ROI2(exclude)=[];
0116   npairs=length(VAR1);
0117 end
0118 
0119 % determine thresholds for each ROI
0120 thresholds=zeros(npairs,2);
0121 if numel(options.spike_threshold)==1
0122   thresholds(:)=options.spike_threshold;
0123 elseif numel(options.spike_threshold)==2
0124   thresholds(:,1)=options.spike_threshold(1);
0125   thresholds(:,2)=options.spike_threshold(2);
0126 end
0127 
0128 kwidth=options.kernel_width;
0129 Ts=options.Ts;
0130 maxlag_time=options.maxlag_time;
0131 
0132 stats=[]; cnt=0;
0133 for pair=1:npairs
0134 
0135   var1=VAR1{pair};
0136   var2=VAR2{pair};
0137   roi1=ROI1{pair};
0138   roi2=ROI2{pair};
0139   fldprefix=[var1 '_' var2];
0140   equal_rois=(isequal(var1,var2) && isequal(roi1,roi2));
0141 
0142   % do nothing if the variables are not present in data
0143   if ~ismember(var1,data.labels) || ~ismember(var2,data.labels)
0144     continue;
0145   end
0146   cnt=cnt+1;
0147 
0148   stats.pairs(cnt).var1=var1;
0149   stats.pairs(cnt).indices1=roi1;
0150   stats.pairs(cnt).var2=var2;
0151   stats.pairs(cnt).indices2=roi2;
0152   stats.pairs(cnt).roi1=options.ROI_pairs{pair,2};
0153   stats.pairs(cnt).roi2=options.ROI_pairs{pair,4};
0154 
0155   t=data.time;
0156   V1=data.(var1)(:,roi1);
0157   V2=data.(var2)(:,roi2);
0158 
0159   n1=length(roi1);
0160   n2=length(roi2);
0161   nt=length(min(t):Ts:max(t));
0162 
0163   % get spike rasters
0164   raster1=dsComputeRaster(t,V1,thresholds(pair,1));
0165   if equal_rois
0166     raster2=raster1;
0167   else
0168     raster2=dsComputeRaster(t,V2,thresholds(pair,2));
0169   end
0170   % raster(:,1) -> spike times
0171   % raster(:,2) -> cell index for each spike
0172 
0173   % calculate fraction of 10ms bins with spikes in both populations
0174   edges=0:maxlag_time:max(t);
0175   spiked1=zeros(1,length(edges));
0176   spiked2=zeros(1,length(edges));
0177   for i=1:length(edges)-1
0178     if size(raster1,1)>0
0179       spiked1(i)=any(raster1(:,1)>edges(i)&raster1(:,1)<=edges(i+1));
0180     end
0181     if size(raster2,1)>0
0182       spiked2(i)=any(raster2(:,1)>edges(i)&raster2(:,1)<=edges(i+1));
0183     end
0184   end
0185   coactive=length(find(spiked1&spiked2))/length(find(spiked1|spiked2));
0186 
0187   % calculate fraction of 10ms bins with spikes in both populations
0188   edges=0:maxlag_time:max(t);
0189   nspiked1=zeros(1,length(edges));
0190   nspiked2=zeros(1,length(edges));
0191   for i=1:length(edges)-1
0192     if size(raster1,1)>0
0193       nspiked1(i)=length(find(raster1(:,1)>edges(i)&raster1(:,1)<=edges(i+1)));
0194     end
0195     if size(raster2,1)>0
0196       nspiked2(i)=length(find(raster2(:,1)>edges(i)&raster2(:,1)<=edges(i+1)));
0197     end
0198   end
0199   th=99;
0200   tmp=nspiked1.*nspiked2;
0201   rm=tmp>prctile(tmp,th);
0202   tmp(rm)=0;
0203   ncoactive=sum(tmp)/(sum(nspiked1(~rm))*sum(nspiked2(~rm)));
0204 
0205   % Calculate instantaneous firing rates for each cell
0206   r1=zeros(nt,n1);
0207   if ~isempty(raster1)
0208     for i=1:n1
0209       [ri,time]=dsNwGaussKernelRegr(t,raster1,i,kwidth,Ts);
0210       r1(:,i)=ri;
0211     end
0212   else
0213     time=0:Ts:max(t)+Ts;
0214     time=time(nearest(time,min(t))):Ts:time(nearest(time,max(t)));
0215   end
0216   if equal_rois
0217     r2=r1;
0218   else
0219     r2=zeros(nt,n2);
0220     if ~isempty(raster2)
0221       for i=1:n2
0222         [ri,time]=dsNwGaussKernelRegr(t,raster2,i,kwidth,Ts);
0223         r2(:,i)=ri;
0224       end
0225     end
0226   end
0227 
0228   % Calculate pairwise cross correlations
0229   maxlags=maxlag_time/(time(2)-time(1));
0230   allxc_possums=nan(n1,n2);
0231   allxc_avgs=nan(n1,n2);
0232   allxc_maxs=nan(n1,n2);
0233   for i=1:n1
0234     xi=r1(:,i);
0235     for j=1:n2
0236       if equal_rois && j>=i, continue; end % exclude symmetric upper triangular matrix and self-correlation
0237       xj=r2(:,j);
0238       [xc,lags]=xcov(xi,xj,maxlags,'coeff');
0239       % take sum over the positive part of the curve
0240       allxc_possums(i,j)=sum(xc(xc>0));
0241       allxc_avgs(i,j)=mean(xc);
0242       allxc_maxs(i,j)=max(xc);
0243     end
0244   end
0245 
0246   % Calculate competition measures
0247   num_spikes1=size(raster1,1);
0248   num_spikes2=size(raster2,1);
0249   dN=num_spikes1-num_spikes2;
0250   dNsumN=dN/(num_spikes1+num_spikes2);
0251   minN=min(num_spikes1,num_spikes2);
0252   maxN=max(num_spikes1,num_spikes2);
0253 
0254   % spectral analysis
0255   dat=dsSelect(data,'roi',{var1,roi1});
0256   dat=dsCalcPower(dat,'time_limits',options.time_limits);
0257   Power_MUA1=dat.([var1 '_Power_MUA']);
0258   Power_SUA1=dat.([var1 '_Power_SUA']);
0259   if ~strcmp(reportUI,'matlab') && exist('nanmean') ~= 2 % 'nanmean is not in Octave's path
0260     try
0261       pkg load statistics; % trying to load octave forge 'statistics' package before using nanmean function
0262     catch
0263       error('nanmean function is needed, please install the statistics package from Octave Forge');
0264     end
0265   end
0266   Power_SUA1.Pxx_mu=nanmean(Power_SUA1.Pxx,2);
0267   Power_SUA1.Pxx_sd=nanstd(Power_SUA1.Pxx,[],2);
0268 
0269   dat=dsSelect(data,'roi',{var2,roi2});
0270   dat=dsCalcPower(dat,'time_limits',options.time_limits);
0271   Power_MUA2=dat.([var2 '_Power_MUA']);
0272   Power_SUA2=dat.([var2 '_Power_SUA']);
0273   Power_SUA2.Pxx_mu=nanmean(Power_SUA2.Pxx,2);
0274   Power_SUA2.Pxx_sd=nanstd(Power_SUA2.Pxx,[],2);
0275 
0276   % store measures
0277   stats.pairs(cnt).raster1=raster1;
0278   stats.pairs(cnt).raster2=raster2;
0279   stats.pairs(cnt).coactivity=coactive;
0280   stats.pairs(cnt).ncoactivity=ncoactive;
0281   stats.pairs(cnt).binned_spikes1=nspiked1;
0282   stats.pairs(cnt).binned_spikes2=nspiked2;
0283   stats.pairs(cnt).bin_edges=edges;
0284   stats.pairs(cnt).r1=r1;
0285   stats.pairs(cnt).r2=r2;
0286   stats.pairs(cnt).Nspikes1=num_spikes1;
0287   stats.pairs(cnt).Nspikes2=num_spikes2;
0288   stats.pairs(cnt).dN=dN;
0289   stats.pairs(cnt).dNsumN=dNsumN;
0290   stats.pairs(cnt).minN=minN;
0291   stats.pairs(cnt).maxN=maxN;
0292   stats.pairs(cnt).xcsum_cells=allxc_possums;
0293   stats.pairs(cnt).xcavg_cells=allxc_avgs;
0294   stats.pairs(cnt).xcmax_cells=allxc_maxs;
0295   stats.pairs(cnt).xcsum_pops=nanmean(allxc_possums(:));
0296   stats.pairs(cnt).xcavg_pops=nanmean(allxc_avgs(:));
0297   stats.pairs(cnt).xcmax_pops=nanmean(allxc_maxs(:));
0298   stats.pairs(cnt).Power_MUA1=Power_MUA1;
0299   stats.pairs(cnt).Power_SUA1=Power_SUA1;
0300   stats.pairs(cnt).Power_MUA2=Power_MUA2;
0301   stats.pairs(cnt).Power_SUA2=Power_SUA2;
0302 
0303 %   stats.([var1 '_roi'])=roi1;
0304 %   stats.([var2 '_roi'])=roi2;
0305 %   stats.([var1 '_raster'])=raster1;
0306 %   stats.([var2 '_raster'])=raster2;
0307 %   stats.([var1 '_inst_rate_cells'])=r1;
0308 %   stats.([var2 '_inst_rate_cells'])=r2;
0309 %   stats.([var1 '_num_spikes'])=num_spikes1;
0310 %   stats.([var2 '_num_spikes'])=num_spikes2;
0311 %   stats.([fldprefix '_dN'])=dN;
0312 %   stats.([fldprefix '_dNsumN'])=dNsumN;
0313 %   stats.([fldprefix '_minN'])=minN;
0314 %   stats.([fldprefix '_xcsum_cells'])=allxc_possums;
0315 %   stats.([fldprefix '_xcavg_cells'])=allxc_avgs;
0316 %   stats.([fldprefix '_xcmax_cells'])=allxc_maxs;
0317 %   stats.([fldprefix '_xcsum_pops'])=nanmean(allxc_possums(:));
0318 %   stats.([fldprefix '_xcavg_pops'])=nanmean(allxc_avgs(:));
0319 %   stats.([fldprefix '_xcmax_pops'])=nanmean(allxc_maxs(:));
0320 end
0321 stats.options=options;
0322 
0323 % calculate instantaneous population firing rate
0324 % pop=unique(raster(:,2));
0325 % [inst_pop_rate,time]=dsNwGaussKernelRegr(t,raster,pop,kwidth,Ts);
0326 
0327 %% auto_gen_test_data_flag argout
0328 if options.auto_gen_test_data_flag
0329   argout = {stats};
0330 
0331   dsUnitSaveAutoGenTestData(argin, argout);
0332 end
dsCalcSpikeSync

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
