merge into master

.gitignore

@@ -0,0 +1,2 @@
+
+SM_SpikeAnalysis.m

AR_Process.m

@@ -0,0 +1,135 @@
+function [B,Pvals,BurstRate,threshVec] = AR_Process(P_Av_Data,Fs,maxThreshold)
+%AR_Process.m
+% Fit a generalized-linear model with a log link function to the burst 
+%  data (which is preprocessed to consider burst activity as a point process,
+%  i.e. as a series of ones and zeros.  The firing rate, mu, is modelled as 
+%  a function of the activity during 20 seconds of lags up to the present 
+%  like ln(mu) = X*b, where X is a matrix of the previous 20 seconds of 
+%  activity.  Thus, the log odds of firing a burst at the present moment 
+%  depends on a linear combination of the previous 20 seconds of bursting 
+%  activity. The output vector b(lag) represents the odds of firing in the 
+%  present given a burst lag-1 bins ago.  So, the value of b(2) gives the 
+%  log odds of firing a burst given a burst 1 time step ago and no bursts 
+%  at any other time
+% See Pillow et al. Nature 2008
+% See also Gerhard et al. PLoS Comp Bio 2013
+%INPUT: P_Av_Data - output of the Ca_im_DataStream code, i.e.
+%        Processed_Data.Av_Data
+%        This code assumes that, for a given night, the same number of 
+%        frames were taken for each video, i.e. each video was recorded for
+%        the same amount of time
+%       Fs - sampling frequency, Hertz
+%       maxThreshold - fraction from 0 to 1, the code will progressively 
+%        threshold the burst data taking only those bursts with greater than
+%        threshold fraction of the total number of active ROIs, and 
+%        maxThreshold tells this threshold procedure when to stop
+%OUTPUT: B - list of coefficients for the link function ln(mu) = X*b
+%         Each row of B represents a burst threshold (see above), while
+%         each column represents the number of lags
+%          B(1,1) represents the log odds of firing a burst given no activity
+%          in the previous 20 seconds and a burst threshold of only a single
+%          ROI, B(1,2) is a lag of one frame, B(1,3) is a lag of two frames
+%       Pvals - statistics on the fit of the model
+%         equal to stats.p, for the p-value of the t test, which
+%         evaluates the null hypothesis that the coefficient b(lag) is equal
+%         to zero
+%       BurstRate - the baseline rate of bursting at each threshold (Hertz) 
+%       threshVec - vector of thresholds up to maxThreshold
+%
+%Created: 2016/02/19 at 24 Cummington, Boston
+%   Byron Price
+%Updated: 2016/04/06
+%  By: Byron Price
+
+threshVec = 0:0.1:maxThreshold;
+%timeLag = 2;
+%numLags = round(timeLag*Fs);
+numLags = 600;
+timeLag = numLags/Fs;
+
+a =15;c=65;
+phiVec = 20*pi:pi/2:31.5*pi;
+numBases = length(phiVec);
+BASIS = zeros(numLags,numBases);
+count = 1;
+t = 1:numLags;
+for phi = phiVec
+    for tt=1:numLags
+        if a*log(t(tt)+c) > (phi-pi) && a*log(t(tt)+c) < (phi+pi)
+            BASIS(tt,count) = 0.5.*cos(a*log(t(tt)+c)-phi)+0.5;
+        else
+            BASIS(tt,count) = 0;
+        end
+    end
+%     plot(t,b)
+%     hold on;
+%     title(sprintf('phi=%d',phi));
+%     pause(0.2)
+    count = count+1;
+end
+
+numNights = size(P_Av_Data,2);
+numVideos = zeros(numNights,1);
+numFrames = zeros(numNights,1);
+totalLength = 0;
+for ii=1:numNights
+    numVideos(ii) = size(P_Av_Data{ii},1);
+    numFrames(ii) = size(P_Av_Data{ii},2);
+    if numFrames(ii) > numLags
+        totalLength = totalLength + numVideos(ii)*(numFrames(ii)-numLags);
+    else
+        display('Must decrease the variable timeLag')
+        return;
+    end
+end
+clear numVideos numFrames;
+
+Time = linspace(1/Fs,timeLag,numLags);
+B = zeros(length(threshVec),numBases+1);
+Pvals = zeros(length(threshVec),numBases+1);
+threshcount = 1;
+for threshold=threshVec
+    Y = zeros(totalLength,1);
+    HIST = zeros(totalLength,numLags);
+    count = 1;
+    for zz=1:numNights
+        Av_Data = P_Av_Data{zz};
+        numVideos = size(Av_Data,1);
+        numFrames = size(Av_Data,2);
+        for jj=1:numVideos
+            newData = Av_Data(jj,:)./max(Av_Data(jj,:));
+            newData(newData<threshold) = 0; % clear all activity below the threshold
+            newData(newData>0) = 1;
+            if count == 1
+                top = 1;
+            else
+                top = bottom+1;
+            end
+            bottom = top+length(newData((numLags+1):end))-1;
+            Y(top:bottom,1) = newData((numLags+1):end);
+            for kk=1:numLags
+                HIST(top:bottom,kk) = newData((numLags-(kk-1)):(numFrames-kk));
+            end
+            count = count+1;
+            clear newData;
+        end
+    end
+    X = HIST*BASIS;
+    [b,~,stats] = glmfit(X,Y,'poisson');
+    B(threshcount,:) = b';
+    Pvals(threshcount,:) = stats.p;
+    threshcount = threshcount+1;
+end
+figure();plot(Time,BASIS*(B(1,2:end)'));xlabel('Time Lag (seconds)');ylabel('Probability of Bursting');
+title('Probability a burst will occur time-lag seconds after a previous burst (Threshold for Bursting of 1 ROI');
+figure();plot(1:numBases,Pvals(1,2:end));xlabel('Parameter');ylabel('P-value');
+title('P-values for each of the model parameters');
+figure();imagesc(Time,threshVec,(BASIS*(B(:,2:end)'))');xlabel('Time Lag (seconds)');ylabel('Threshold');
+title('Probability of bursting (all thresholds)');
+%figure();plot(Time,B(1,2:end));figure();plot(Time,B(end,2:end));
+BurstRate = exp(B(:,1)).*Fs;
+figure();plot(threshVec,BurstRate);ylabel('Baseline Burst Rate (Hz)');
+xlabel('Threshold (fraction of total ROIs active)');
+title('Burst rate as a function of burst threshold');
+end
+

Byrons_Avalanche_Analysis.m

@@ -0,0 +1,182 @@
+function [Av_Size,mean_Size,Av_IEI,mean_IEI,StDev_IEI,Time] = Byrons_Avalanche_Analysis(Av_Data,window,maxThreshold)
+%Avalanche_Analysis.m
+%   Take "Av_Data" matrix from "From_Raw_To_Traces.m" output and count 
+%    avalanches using different definitions and thresholds.
+% Created: 2015/12/15 at 24 Cummington, Boston
+%   Byron Price
+% Updated: 2015/12/15
+%  by: Byron Price
+% 
+% INPUT: Av_Data - a single matrix size videos by frames
+%
+%        window - for avalanche definition
+%        maxThreshold - also for avalanche definition, the code will
+%        iterate from a threshold of 1 up to maxThreshold, where each
+%        threshold signifies the minimum number of ROIs that must be active
+%        in a single frame for an avalanche to be counted 
+%
+% OUTPUT: Av_Size - a matrix with count data for each of the possible
+%         thresholds and avalanche sizes ... if Av_Size(1,1) = 50, then
+%         there were 50 instances of an avalanche of size 1
+%           * matrix size threshold by max number of ROIs
+%
+%         mean_Av_Size - the mean size of an avalanche at each threshold
+%         ... this isn't so interesting as we force the minimum up with
+%         each increase in the treshold
+%
+%         Av_IEI - a matrix with count data for each of the possible
+%         avalanche thresholds and inter-event intervals 
+%         = There are # avalanches with an IEI of 1 frame; 2 frames; etc.
+%
+%         mean_Av_IEI - the mean inter-event interval between avalanches of
+%         size at least "threshold"
+%
+%         Time - vector of times corresponding to the bins in Av_IEI,
+%         try plotting  plot(Time, Av_IEI(1,:))  for a plot of the number
+%         of instances of an inter-event interval at each of the times in
+%         Time assuming only those bins with greater than or equal to 1
+%         active ROI is counted as an avalanche ... Av_IEI(2,:) would give
+%         you the same plot but with IEI's between those avalanches greater
+%         than or equal to a size of two ROIs
+%
+% AVALANCHE DEFINITION
+%   If window = 0, then the definition is ... the summed activity across
+%   all ROIs in a single time bin, if that sum is greater than the value of 
+%   threshold.  So, if 10 ROIs are active in time bin 20,
+%   then that's an avalanche of size 10. Each bin is an island and
+%   contiguous bins are never counted as part of the same avalanche.  As
+%   threshold increases, the minimum number of active ROIs per bin
+%   increases.
+%      Example: For a threshold of 1, a sequence from Av_Data like 
+%      [0,0,1,3,4,0,5,6,9,0,0,0] would have 4 inter-event intervals of 
+%      0 seconds (contiguous bins active), 1 IEI of 0.0333 seconds (1 bin 
+%      separation), 1 IEI of 0.0666 seconds and at the end 1 of 0.0999 seconds.
+%      That would be for threshold 1 ... Av_IEI(1,:). For a threshold of
+%      2, stored in Av_IEI(2,:), you would have only 4 avalanches and the
+%      IEIs would change.
+%
+%   If window = 1, then the definition is ... the summed activity across
+%   all ROIs and time bins for which at least "threshold" ROIs are active.
+%   So, for threshold = 1, contiguous bins with at least one ROI active are
+%   counted as a single avalanche.
+%      Example: Given the sequence above [0,0,1,3,4,0,5,6,9,0,0,0] and a
+%      threshold of 1, the inter-event inetervals would be the same, except
+%      there would be no IEIs of 0 seconds because IEIs of 0 seconds mean
+%      two consecutive bins were active and both bins are counted in the
+%      same avalanche.  The sizes, however, would change, with 1 avalanche
+%      of size 1+3+4=8 ROIs and another avalanche of size 5+6+9=20 ROIs.
+%
+%   As 'window' increases, the number of bins with activity below threshold
+%   between bins with activity above or equal to threshold grows.  Thus,
+%   for window = 2 and the sequence above [0,0,1,3,4,0,5,6,9,0,0,0], there
+%   would now be only one avalanche of size 28.  The IEIs would be
+%   basically the same but with 1 IEI of 2 bins (0.0666 seconds) and one of
+%   3 bins (0.0999 seconds).  
+%
+Fs = 30;
+numVideos = size(Av_Data,1);
+T = size(Av_Data,2);
+Time = 0:1/Fs:T/Fs;
+maxROIs = 500;   
+
+Av_Size = zeros(maxThreshold,maxROIs); % avalanche size 
+%Av_Size(:,:,2) = ones(maxThreshold,1)*(1:maxROIs);
+mean_Size = zeros(maxThreshold,1);
+
+Av_IEI = zeros(maxThreshold,length(Time)); % inter-event interval (seconds), or time 
+                % between avalanches
+mean_IEI = zeros(maxThreshold,1);
+
+
+for threshold = 1:maxThreshold
+    Sizes = [];
+    IEIs = [];
+    for i=1:numVideos
+            icount = 0;
+            scount = 0;
+            if window == 0
+                for t=1:T
+                    if Av_Data(i,t) >= threshold || t == 1
+                        if Av_Data(i,t) >= threshold
+                            Av_Size(threshold,Av_Data(i,t),1) = Av_Size(threshold,Av_Data(i,t),1)+1;
+                            Sizes = [Sizes,Av_Data(i,t)];
+                        end
+                        count = 0;
+                        k = t+1;
+                        while Av_Data(i,k) < threshold && k < T
+                            k = k+1;
+                            count = count+1;
+                        end
+                        if count < T-2
+                            Av_IEI(threshold,count+1,1) = Av_IEI(threshold,count+1,1)+1;
+                            IEIs = [IEIs,Time(count+1)];
+                        end
+                    end
+                end
+
+            else
+                for t=1:T
+                    if t >= (T-(window-1))
+                        if t == T
+                            if scount > 0
+                                scount = sum(squeeze(Av_Data(i,t-scount:end)));
+                                Av_Size(threshold,scount,1) = Av_Size(threshold,scount,1)+1;
+                                Sizes = [Sizes,scount];
+                            end
+                            if icount < T-1
+                                icount = icount+(window-1);
+                                Av_IEI(threshold,icount+1,1) = Av_IEI(threshold,icount+1,1)+1;
+                                IEIs = [IEIs,Time(icount+1)];
+                            end
+                        else
+                            len = length(Av_Data(i,t:end));
+                            if (squeeze(Av_Data(i,t:end)) < threshold) == ones(1,len)
+                                icount = icount+1;
+                                if scount > 0
+                                    scount = sum(squeeze(Av_Data(i,t-scount:t-1)));
+                                    Av_Size(threshold,scount,1) = Av_Size(threshold,scount,1)+1;
+                                    Sizes = [Sizes,scount];
+                                    scount = 0;
+                                end
+                            else
+                                scount = scount+1;
+                                if icount > 0
+                                    icount = icount+(window-1);
+                                    Av_IEI(threshold,icount+1,1) = Av_IEI(threshold,icount+1,1)+1;
+                                    IEIs = [IEIs,Time(icount+1)];
+                                    icount = 0;
+                                end
+                            end
+                        end
+                    else
+                        if (squeeze(Av_Data(i,t:(t+(window-1)))) < threshold) == ones(1,window)
+                            icount = icount+1;
+                            if scount > 0
+                                scount = sum(squeeze(Av_Data(i,t-scount:t-1)));
+                                Av_Size(threshold,scount,1) = Av_Size(threshold,scount,1)+1;
+                                Sizes = [Sizes,scount];
+                                scount = 0;
+                            end
+                        else
+                            scount = scount+1;
+                            if icount > 0
+                                icount = icount+(window-1);
+                                Av_IEI(threshold,icount+1,1) = Av_IEI(threshold,icount+1,1)+1;
+                                IEIs = [IEIs,Time(icount+1)];
+                                icount = 0;
+                            end
+                        end
+                    end
+                end
+            end
+    end
+    mean_Size(threshold,1) = mean(Sizes);
+    mean_IEI(threshold,1) = mean(IEIs);
+    StDev_IEI(threshold,1) = std(IEIs);
+
+
+end
+save('Spontaneous_Data_Averages2_2.mat','Av_Data','Av_IEI','Av_Size','mean_Size','mean_IEI','StDev_IEI');
+end
+
+