gin commit from DESKTOP-593CF0Q

New files: 26

codes/Behavioral analysis/Fig1NormLearningScoreAllData.m

@@ -0,0 +1,59 @@
+load('CtrlDreaddCumValsNew.mat', 'untrained','pre_norm_last', 'post_norm_last','Ctrl','Exp')
+
+
+%%
+%%                                                 % Mean
+%dci  = std(data)*tinv(0.975,size(data,1)-1);                        % Confidence Intervals
+xt = [1:5];                                                         % X-Ticks
+
+lb = [xt'-ones(5,1)*0.2,  xt'+ones(5,1)*0.2]; % Long Bar X
+figure
+plot(xt(1), untrained, 'o', 'color', [0.5 0.5 0.5],'MarkerSize', 10)
+hold on
+plot(xt(2),Ctrl , 'ko','MarkerSize', 10)
+plot(xt(3), Exp, 'ro','MarkerSize', 10)
+plot(xt(4), pre_norm_last, 'ro','MarkerSize', 10)
+plot(xt(5), post_norm_last, 'ro','MarkerSize', 10)
+
+
+dmean = [mean(untrained), mean(Ctrl), mean(Exp),  mean(pre_norm_last), mean(post_norm_last)];
+%dmedian = [median(untrained), median(Ctrl), median(Exp),  median(expert)];
+
+for k1 = 1:5
+    %plot(lb(k1,:), [1 1]*dmedian(k1), '-k')
+    plot(lb(k1,:), [1 1]*dmean(k1), '-b')
+end
+hold off
+set(gca, 'XTick', xt, 'XTickLabel', {'Untrained','Ctrl','hM4Di/CNO', 'Expert preCNO','Expert postCNO'})
+xlabel('Group')
+ylabel('Normalized learning score')
+
+xlim([0.5 5.5])
+ylim([-1 1])
+set(gca,'YTick', [-1:0.5:1])
+
+%% Statistics
+[h, p] = swtest(Ctrl)
+[h, p] = swtest(Exp)
+
+
+% Parametric
+group = [ones(20,1); ones(12,1)+1; ones(5,1)+2];
+[p,tbl,stats] = anova1([Ctrl; Exp ; untrained],group,'off')
+
+c = multcompare(stats,'CType', 'bonferroni')
+
+
+[h p] =ttest2(Ctrl, Exp)
+
+%% non-parametric
+[p,tbl,stats] = kruskalwallis([Ctrl; Exp ; untrained],group,'off')
+
+c = multcompare(stats,'CType', 'dunn-sidak')
+
+
+[p h] =ranksum(Ctrl, Exp)
+
+[p h] = ranksum(Ctrl, untrained)
+[p h] = ranksum(untrained, Exp)  
+[p h] = ranksum(Ctrl, post_norm_last)

codes/Behavioral analysis/Fig1PlotTotalCumSumVals.m

@@ -0,0 +1,38 @@
+load('CtrlDreaddCumValsNew.mat', 'Control_cumvals', 'Exp_cumvals','Naive_cumvals' )
+
+%%
+figure
+errorbar(nanmean(Control_cumvals,2),nanstd(Control_cumvals,[],2)/sqrt(20),'k')
+hold on
+errorbar(nanmean(Exp_cumvals(1:100,:),2),nanstd(Exp_cumvals(1:100,:),[],2)/sqrt(12),'r')
+hold on
+errorbar(nanmean(Naive_cumvals(1:100,:),2),nanstd(Naive_cumvals(1:100,:),[],2)/sqrt(5),'b')
+
+ylabel('Cummulative Sum')
+xlabel('trial #')
+
+xlim([0 100])
+
+
+%%
+for i =1:size(Control_cumvals,2)
+    plot(Control_cumvals(:,i),'k')
+    hold on
+end
+
+
+for i =1:size(Exp_cumvals,2)
+    plot(Exp_cumvals(:,i),'r')
+    hold on
+end
+
+
+
+for i =1:size(Naive_cumvals,2)
+    plot(Naive_cumvals(:,i),'b')
+    hold on
+end
+
+legend({'Control (n=20)', 'hM4Di (n=12)',  'Naive (n=5)'})
+ylim([-100 100])
+set(gca,'ytick', -100:50:100)

codes/Behavioral analysis/Fig3CtrlBacSST.m

@@ -0,0 +1,77 @@
+load('Figure3Data.mat')
+
+%%
+SST_cumvals = nan(150,size(SSTCumVals,1));
+for i = 1:size(SSTCumVals,1)
+    cumvals = SSTCumVals{i};
+    if size(cumvals,1) < 150
+        SST_cumvals(1:size(cumvals,1),i) = cumvals;
+    else
+        SST_cumvals(:,i) = cumvals(1:150);
+    end
+end
+
+%%
+BacL1_cumvals = BacL1_cumvals(1:100,:);
+SST_cumvals = SST_cumvals(1:100,:);
+
+%%
+figure
+errorbar(nanmean(Control_cumvals,2),nanstd(Control_cumvals,[],2)/sqrt(11),'k')
+
+hold on
+errorbar(nanmean(BacL1_cumvals,2),nanstd(BacL1_cumvals,[],2)/sqrt(6),'r')
+
+errorbar(nanmean(SST_cumvals,2),nanstd(SST_cumvals,[],2)/sqrt(6),'b')
+
+xlim([0 100])
+ylabel('Cummulative Sum')
+xlabel('trial #')
+
+
+
+%%
+
+
+for i =1:size(BacL1_cumvals,2)
+    plot(BacL1_cumvals(:,i),'r')
+    hold on
+end
+
+
+for i =1:size(SST_cumvals,2)
+    plot(SST_cumvals(:,i),'b')
+    hold on
+end
+legend({'Control (n=20)', 'BacL1 (n=6)',  'SST (n=6)'})
+ylim([-100 100])
+set(gca,'YTick', -100:50:100)
+%%
+
+figure
+PlotAllDataPoints3(Ctrl, BacL1, SSTNormCumSumLastVal)
+
+
+set(gca, 'XTick',1:3, 'XTickLabel', {'Ctrl','Baclofen', 'SST'})
+
+ylabel('Normalized learning score')
+xlim([0.5 3.5])
+
+%%
+group = [ones(20,1); ones(6,1)+1; ones(6,1)+2];
+
+[p,tbl,stats] = kruskalwallis([Ctrl; BacL1; SSTNormCumSumLastVal'],group,'off')
+
+c = multcompare(stats,'CType', 'dunn-sidak')
+
+
+%%
+[p1 h] = ranksum(Ctrl, BacL1)
+
+[p2 h] = ranksum(Ctrl, SSTNormCumSumLastVal)
+
+text(1.8, 0.7, ['p = ', num2str(p1)])
+text(2.8, 0.7, ['p = ', num2str(p2)])
+
+ylim([-1 1])
+set(gca,'YTick', [-1:0.5:1])

codes/Behavioral analysis/Fig3NormLearningScoreCtrlSST.m

@@ -0,0 +1,10 @@
+%Ctrl = [0.317596566523605; 0.281105990783410; 0.574712643678161; 0.487603305785124 ;0.461538461538462; 0.707006369426752];
+%SST = [-0.651785714285714; -0.113772455089820; -0.467032967032967; 0.166666666666667; 0.0941828254847646; 0];
+
+                                                 % Mean
+figure                                                 
+PlotAllDataPoints(Ctrl,SST)
+set(gca, 'XTickLabel', {'Ctrl','SST'})
+ylabel('Normalized learning score')
+ylim([-1 1])
+set(gca,'YTick', [-1 -0.5 0 0.5 1])

codes/Behavioral analysis/GetNormCumSumLastVal.m

@@ -0,0 +1,17 @@
+function [CumSumLastVal, NormCumSumLastVal] = GetNormCumSumLastVal (Group)
+ %input: group is  a cell containing cumulative sum values of the group 
+
+AnimalNum = size(Group,1);
+TrialNum = cellfun(@length, Group);
+
+CumSumLastVal = zeros(AnimalNum,1);
+for i = 1:AnimalNum
+    CumSumLastVal(i) = Group{i}(end);
+end
+
+NormCumSumLastVal = CumSumLastVal./TrialNum;
+
+
+
+
+

codes/Behavioral analysis/HitMissMicStim.m

@@ -0,0 +1,16 @@
+function MicStim = HitMissMicStim(MicStim)
+    reactiontime = MicStim(:,2);
+    
+    ReactionTimeLow = 0.1;
+    ReactionTimeHigh = 1.2;
+    
+    Aborted = find(reactiontime >0 & reactiontime <= ReactionTimeLow);
+    MicStim(Aborted,3) = -1; % Aborted trials = -1
+    
+    Hits = find(reactiontime>ReactionTimeLow & reactiontime<=ReactionTimeHigh);
+    MicStim(Hits,3) = 1; % Hit trials = 1;
+    
+    Misses = find(reactiontime<0 | reactiontime>ReactionTimeHigh);
+    MicStim(Misses,3) = 0; %Miss trials = 0
+    
+end

codes/Behavioral analysis/JScumvals.m

@@ -0,0 +1,26 @@
+function cumvals = JScumvals(MicStim)
+%% Cummulative response
+
+
+trialNum = length(MicStim);
+
+
+cum = 0;
+cumvals=[];
+
+for i = 1:trialNum
+    
+    if MicStim(i,3) ==1  %Hit
+        cum = cum+1;
+    elseif MicStim(i,3) == 0 %Miss
+        cum = cum-1;
+    else
+        cum = cum;
+    end
+    cumvals = [cumvals; cum];
+    
+end
+
+figure
+plot(cumvals)
+end

codes/Behavioral analysis/Supp2NormLearningScoreCtrlHPC.m

@@ -0,0 +1,57 @@
+load('SuppLidoHCCumVals.mat')
+%%
+xt = [1:4];                                                         % X-Ticks
+
+lb = [xt'-ones(4,1)*0.2,  xt'+ones(4,1)*0.2]; % Long Bar X
+
+plot(xt(1), RatCtrlNormCumSumLastVal, 'ko','MarkerSize', 10)
+hold on
+plot(xt(2), RatHPCExpNormCumSumLastVal , 'ro','MarkerSize', 10)
+plot(xt(3), MouseCtrl , 'ko','MarkerSize', 10)
+plot(xt(4), MouseHPC , 'ro','MarkerSize', 10)
+
+
+dmean = [mean(RatCtrlNormCumSumLastVal), mean(RatHPCExpNormCumSumLastVal), mean(MouseCtrl), mean(MouseHPC)];
+
+
+for k1 = 1:4
+    
+    plot(lb(k1,:), [1 1]*dmean(k1), '-b')
+end
+
+xlim([0.5 4.5])
+ylim([-1 1])
+set(gca, 'XTick', xt, 'XTickLabel', {'Rat Ctrl','Rat HPC', 'Mouse Ctrl', 'MouseHPC'})
+
+set(gca, 'YTick', -1:0.5:1)
+
+%%
+[p1 h] = ranksum(RatCtrlNormCumSumLastVal, RatHPCExpNormCumSumLastVal)
+
+[p2 h] = ranksum(MouseCtrl, MouseHPC)
+
+[p h] = ranksum(RatCtrlNormCumSumLastVal, MouseCtrl)
+
+text(0.8, 0.8, ['p=', num2str(p1)])
+text(2.8, 0.8, ['p=', num2str(p2)])
+
+%%
+%%
+figure
+errorbar(nanmean(MouseCtrl_cumvals,2),nanstd(MouseCtrl_cumvals,[],2)/sqrt(20),'k')
+hold on
+
+errorbar(nanmean(MouseHPC_cumvals(1:100,:),2),nanstd(MouseHPC_cumvals(1:100,:),[],2)/sqrt(3),'b')
+
+
+errorbar(nanmean(RatCtrl_cumvals(1:100,:),2),nanstd(RatCtrl_cumvals(1:100,:),[],2)/sqrt(24),'b')
+
+errorbar(RatHPCExpCumSumValsAvg(1:100), RatHPCExpCumSumValsSEM(1:100),'m')
+
+xlim([0 100])
+ylabel('Cummulative Sum')
+xlabel('trial #')
+
+ylim([-50 100])
+
+legend({'Mouse Ctrl (n=20)',   'Mouse HPC (n=3)', 'Rat Ctrl (n=24)', 'RatHPC (n=7)'})

codes/Behavioral analysis/Supp2RatCtrlPRh.m

@@ -0,0 +1,45 @@
+load('SuppLidoPRhCumVals.mat')
+%%
+PlotAllDataPoints3(RatCtrlNormCumSumLastVal, PRhExp, PRhPost)
+
+
+ylim([-1 1])
+set(gca,'xtick',1:3, 'XTickLabel', {'Rat Ctrl','PRh', 'PRhtrained'})
+
+set(gca, 'YTick', -1:0.5:1)
+
+%%
+[p1 h] = ranksum(RatCtrlNormCumSumLastVal, PRhExp)
+
+[p2 h] = ranksum(RatCtrlNormCumSumLastVal, PRhPost)
+
+
+%text(0.8, 0.8, ['p=', num2str(p1)])
+%text(2.8, 0.8, ['p=', num2str(p2)])
+
+%% non-parametric
+group = [ones(24,1); ones(4,1)+1; ones(6,1)+2];
+
+[p,tbl,stats] = kruskalwallis([RatCtrlNormCumSumLastVal; PRhExp ; PRhPost],group,'off')
+
+c = multcompare(stats,'CType', 'dunn-sidak')
+
+
+%%
+figure
+errorbar(nanmean(RatCtrl_cumvals,2),nanstd(RatCtrl_cumvals,[],2)/sqrt(24),'k')
+hold on
+
+errorbar(nanmean(PRhExp_cumvals,2),nanstd(PRhExp_cumvals,[],2)/sqrt(4),'b')
+
+
+errorbar(nanmean(PRhPost_cumvals,2),nanstd(PRhPost_cumvals,[],2)/sqrt(6),'g')
+
+
+xlim([0 100])
+ylabel('Cummulative Sum')
+xlabel('trial #')
+
+ylim([-50 100])
+
+legend({'Rat Ctrl (n=24)',   'PRh (n=4)', 'PRh trained (n=6)'})

codes/Behavioral analysis/Supp6CtrlPOm.m

@@ -0,0 +1,74 @@
+load('SuppPOmData.mat')
+
+%%
+POm_cumvals = nan(150,size(POmCumVals,1));
+for i = 1:size(POmCumVals,1)
+    cumvals = POmCumVals{i}';
+    if size(cumvals,1) < 150
+        POm_cumvals(1:size(cumvals,1),i) = cumvals;
+    else
+        POm_cumvals(:,i) = cumvals(1:150);
+    end
+end
+
+%%
+POm_cumvals = POm_cumvals(1:100,:);
+
+%%
+figure
+errorbar(nanmean(Control_cumvals,2),nanstd(Control_cumvals,[],2)/sqrt(20),'k')
+
+hold on
+
+errorbar(nanmean(POm_cumvals,2),nanstd(POm_cumvals,[],2)/sqrt(7),'b')
+
+xlim([0 100])
+ylabel('Cummulative Sum')
+xlabel('trial #')
+
+
+%%
+% figure
+% boundedline(1:150, nanmean(Control_cumvals,2), nanstd(Control_cumvals,[],2)/sqrt(11),'k')
+% hold on
+% boundedline(1:150, nanmean(BacL1_cumvals,2),nanstd(BacL1_cumvals,[],2)/sqrt(12),'r')
+% hold on
+% boundedline(1:150, nanmean(POm_cumvals,2),nanstd(POm_cumvals,[],2)/sqrt(5),'b')
+% 
+% ylabel('Cummulative Sum')
+% xlabel('trial #')
+
+%%
+for i =1:size(Control_cumvals,2)
+    plot(Control_cumvals(:,i),'k')
+    hold on
+end
+
+
+
+for i =1:size(POm_cumvals,2)
+    plot(POm_cumvals(:,i),'b')
+    hold on
+end
+legend({'Control (n=20)',   'POm (n=7)'})
+
+ylim([-50 100])
+%%
+
+figure
+PlotAllDataPoints(Ctrl, POm)
+
+
+set(gca, 'XTick',1:2, 'XTickLabel', {'Ctrl', 'POm'})
+
+ylabel('Normalized learning score')
+xlim([0.5 2.5])
+
+
+[p h] = ranksum(Ctrl, POm)
+
+text(1.8, 0.7, ['p = ', num2str(p)])
+
+
+ylim([-1 1])
+set(gca,'YTick', [-1:0.5:1])

codes/Behavioral analysis/Supp6WTCNO.m

@@ -0,0 +1,34 @@
+load('CtrlDreaddCumValsNew.mat', 'Ctrl', 'CNO')
+%%
+Expert = Ctrl(1:3);
+WT = setdiff(Ctrl, [CNO;Expert]);
+%%
+
+figure
+PlotAllDataPoints3(WT, CNO, Expert)
+
+
+set(gca, 'XTick',1:3, 'XTickLabel', {'WT', 'CNO','Expert'})
+
+ylabel('Normalized learning score')
+xlim([0.5 3.5])
+
+%%
+x = [WT; CNO; Expert];
+group  = [ones(length(WT),1); ones(length(CNO),1)+1; ones(length(Expert),1)+2]
+p = kruskalwallis(x,group)
+
+
+
+
+%%
+[p1 h] = ranksum(WT, CNO)
+
+text(1.2, 0.7, ['p = ', num2str(p1)])
+
+[p2 h] = ranksum(WT, Expert)
+text(2.2, 0.7, ['p = ', num2str(p2)])
+
+
+ylim([-0.5 1])
+set(gca,'YTick', [-0.5:0.5:1])

codes/Electrophysiology analysis/AllCellsBurstPSTH.m

@@ -0,0 +1,52 @@
+function Group = AllCellsBurstPSTH(Group,preStim,postStim, binsize, mkplt)
+% ex. Controls = AllCellsBurstPSTH(Controls)
+% Compute burst psth and add a field of the input structure
+%       INPUT: 
+%               Group - (required) Experimental group structure containing
+%               all cell information
+%               mkplt - 1 to plot 0 to not plot
+%       OUTPUT:
+%               Group - original input structure plus a new field bpsth
+            
+
+for i = 1:length(Group)
+    BurstTimes = Group(i).BurstTimes;
+    cellname = Group(i).CellName;
+    
+    if length(Group(i).Mixes) ~=0
+        Mixes = Group(i).Mixes;
+        Mixes(:,3) = 10;
+        Microstims = [Group(i).Microstims; Mixes];
+    else 
+        Microstims = [Group(i).Microstims];
+    end
+    
+    %HighIntensity = find(Microstims(:,3) >= 40);
+    %LowIntensity = find(Microstims(:,3) < 40);
+    HitTimes = find(Microstims(:,2) <= 1.2 & Microstims (:,2) > 0.1);
+    
+    Hits = Microstims(HitTimes,1);
+    %HighStimTimes = Microstims(intersect(HighIntensity, HitTimes),1);
+    %LowStimTimes = Microstims(intersect(LowIntensity, HitTimes),1);
+    
+ 
+    
+    
+    %[bpsthNormHigh,numTrialsHigh edges] = RasterPSTH (HighStimTimes, BurstTimes(1:2:end), preStim, postStim, binsize,'mkplt', false );
+    %saveas(gcf, [cd, '\RasterPSTH\', cellname, 'RasterPSTHhigh.fig'])
+
+    
+    %[bpsthNormLow,numTrialsLow edges] = RasterPSTH (LowStimTimes, BurstTimes(1:2:end), preStim, postStim, binsize,'mkplt', false );
+    %saveas(gcf, [cd, '\RasterPSTH\', cellname, 'RasterPSTHlow.fig'])
+    
+    % Total microstim
+    [bpsthNorm,numHits, edges] = RasterPSTH (Hits, BurstTimes(1:2:end), preStim, postStim, binsize,'mkplt', mkplt);
+    
+    %Group(i).bpsthNormHigh= bpsthNormHigh;
+    %Group(i).numTrialsHigh = numTrialsHigh;
+    %Group(i).bpsthNormLow = bpsthNormLow;
+    %Group(i).numTrialsLow = numTrialsLow;
+    Group(i).bpsthNorm= bpsthNorm;
+    Group(i).numHits= numHits;
+    %close all
+end

codes/Electrophysiology analysis/AllDoron_ForScience2020.py

@@ -0,0 +1,25 @@
+import matplotlib as mpl
+mpl.use('Agg')
+#import matplotlib.pyplot as plt
+
+from pylab import *
+from numpy import *
+from scipy import io
+import os
+import seaborn as sns
+sns.set(style='ticks')
+
+# 3 sigma was chosen to fit opto protocol
+nsigmaON = 3.0
+nsigmaOFF = 3.0
+task = 'Figures_ONOFF_later_3typesFRHits_SimplerCode/' #output path
+datapath = 'ExpertS1Cells/' # input path
+execfile('CellCategories_3typesFRHits_SimpleCode.py')
+
+task = 'Figures3types_Naive_FRHits_SimplerCode/' #output
+datapath = 'NaiveS1Cells/' # input
+execfile('CellCategories_3typesFRHits_SimpleCode.py')
+
+
+
+

codes/Electrophysiology analysis/CellCategories_3typesFRHits_SimpleCode.py

@@ -0,0 +1,216 @@
+
+def create_2tracefig(ylims):
+    from pylab import figure
+    fig1 = figure(num=1,figsize = (5,1.8), facecolor = 'w', dpi = 150, edgecolor = 'w')
+    fig1.clf()
+    ax1=axes([0.12,0.3,0.75,0.65])
+    sns.despine()
+    ax2 = ax1.twinx()
+    ax2.spines["right"].set_edgecolor('r')
+    ax1.set_ylabel('Rate [Hz]')
+    ax1.set_xlabel('Time [s]')
+    ax2.set_ylabel('Burst. Prob.')
+    ax1.set_xticks([-1.5,-1,-.5,0,.5,1,1.5])
+    ax2.tick_params(axis='y', colors='r')
+    ax1.spines['top'].set_visible(False)
+    ax2.spines['top'].set_visible(False)
+    ax2.set_ylim((-0.05,.9))
+    ax1.set_xlim((-1.4,1.4))
+    ax1.set_ylim(ylims)
+    return fig1,ax1,ax2
+def train2binned(train,dt,smboxsize):
+    if train.shape[0]==0:
+        smrate = zeros(train.shape[1])
+    else:
+        rate = mean(train,axis=0)/dt
+        nbins = floor(float(int(len(rate)))/smboxsize)
+        smrate = rate[:int(nbins*smboxsize)].reshape((int(nbins),smboxsize))
+        smrate = array([mean(smrate,axis=1)])
+        smrate = smrate.repeat(smboxsize,axis=0)
+        smrate = transpose(smrate).reshape(int(nbins*smboxsize),1)
+        smrate = append(smrate,ones(len(rate)-int(nbins*smboxsize))*mean(rate[int(nbins*smboxsize):]))
+    return smrate
+def create_1tracefig(ylims):
+    fig1 = figure(num=1,figsize = (5,1.8), facecolor = 'w', dpi = 150, edgecolor = 'w')
+    fig1.clf()
+    ax1=axes([0.12,0.3,0.75,0.65])
+    sns.despine()
+    ax1.set_ylabel('Rate [Hz]')
+    ax1.set_xlabel('Time [s]')
+    ax1.spines['top'].set_visible(False)
+    ax1.set_ylim((-1,14))
+    ax1.set_xticks([-1,-.5,0,.5,1])
+    ax1.set_ylim(ylims)
+    return fig1,ax1
+def S2Btr(spktre,t,dt,npyr,BstTh):
+    """
+    transform set of spike trains into set of burst trains
+    """
+    if spktre.shape[0]==0:
+        bstr = zeros((1,spktre.shape[1]))
+        sstr = zeros((1,spktre.shape[1]))
+    else:
+        bstr = spktre*1.0
+        sstr = spktre*1.0
+        for ind in range(spktre.shape[0]):
+            isise = array([1e3])
+            spkse = where(spktre[ind,:])[0]
+            isise = append(isise,diff(spkse)*dt)
+            isise = append(isise,1e3)
+            if len(spkse)>1:
+                WhichAPinBst = where(logical_or(isise[1:]<BstTh,isise[:-1]<BstTh))[0] # spk is in bst if either isi after or isi before is smalle than bst isi-threshold
+                Is1stinBst = logical_and(isise[1:]<BstTh,isise[:-1]>BstTh)
+                sstr[ind,spkse[WhichAPinBst]] = sstr[ind,spkse[WhichAPinBst]]*0
+                WhichAPnot1stBsty = where(Is1stinBst!=1)[0] # not first spike in burst isi after bursty but not before
+                bstr[ind,spkse[WhichAPnot1stBsty]] = bstr[ind,spkse[WhichAPnot1stBsty]]*0 # 
+            else:
+                bstr[ind,spkse] = bstr[ind,spkse]*0
+
+    return bstr,sstr
+
+
+legroom1 = -3000 # half window around stim for deciding ON or OFF, T/2-legroom1:T/2+legroom2. So -3000:4000 means +0.3:0.4s
+legroom2 = 4000
+BstTh = 15e-3
+smboxsize = 1001 # 100 ms binning
+lim1=[0.12,0.3,0.75,0.65]
+dt = 1e-4 # discretize in 0.1 ms time bins, dt is in seconds now
+T = 4.0
+t = arange(0,T,dt)
+allhitreps=[]
+L=int(T/dt)
+
+# Loading data
+cells = os.listdir(datapath)
+ncells = len(cells)
+firingcells = array([])
+for i in range(ncells):
+    if cells[i][-4:]=='.mat':
+        WS = io.loadmat(datapath+cells[i])
+        StimTimes = WS.get('MicrostimTimes')[:,0]
+        if logical_and(len(StimTimes)>0,cells[i][-4:]=='.mat'): # only cells containing at least one stim
+            firingcells=append(firingcells,i)
+ncells = len(firingcells)
+
+# prepare data arrays and counters
+FratesNR_Hits = zeros((ncells,L))
+FratesON_Hits = zeros((ncells,L))
+FratesOFF_Hits = zeros((ncells,L))
+BratesNR_Hits = zeros((ncells,L))
+BratesON_Hits = zeros((ncells,L))
+BratesOFF_Hits = zeros((ncells,L))
+nHits = zeros((ncells,1))
+nMissed = zeros((ncells,1))
+nrep = zeros((ncells,1))
+counton = 0
+countoff = 0
+countnr = 0
+nON_Hits=0
+nON_Missed=0
+nNR_Hits=0
+nNR_Missed=0
+nOFF_Hits =0
+nOFF_Missed=0
+oninds =array([])
+offinds = array([])
+
+# process each cell type
+for i in firingcells.astype('int'):
+    WS = io.loadmat(datapath+cells[i])
+    spks = WS.get('SpikeTimes')[:,0]
+    StimTimes = WS.get('MicrostimTimes')[:,0]
+    Hit = WS.get('Hit')[:,0]
+    ind = where(firingcells==i)
+    hitreps = where(Hit==1)[0]
+    nHits[ind] = sum(Hit==1)
+    nMissed[ind] = sum(Hit!=1)
+    nrep[ind] = nHits[ind]+nMissed[ind]
+    allhitreps = [allhitreps,hitreps]
+    spktr = zeros((int(len(Hit)),int(T/dt)),dtype=dtype('bool') )
+    for j in range(len(spks)):
+    	dvec = abs(spks[j]-StimTimes)
+    	mdvec = min(dvec)
+    	rep = where(dvec==mdvec)[0][0]
+    	if logical_and(mdvec<T/2,Hit[rep]==1):	########### Hit or Miss selection	
+    		spktr[rep,int((spks[j]-StimTimes[rep]+T/2)/dt)]=True		
+    	elif mdvec<T/2:
+    		spktr[rep,int((spks[j]-StimTimes[rep]+T/2)/dt)]=True
+
+    btr,sstr = S2Btr(spktr,t,dt,nrep[ind],BstTh)
+    frate = train2binned(spktr,dt,smboxsize)
+    brate = train2binned(btr,dt,smboxsize)
+    spktrHits = spktr[hitreps,:]
+    btrHits,sstrHits = S2Btr(spktrHits,t,dt,nHits[ind],BstTh) # separates bursts from single spikes based in ISI interval
+    frateHits = train2binned(spktrHits,dt,smboxsize)
+    brateHits = train2binned(btrHits,dt,smboxsize)
+
+    e0 = mean(frateHits[:int((T/2-.02)/dt)])
+    sde = std(frateHits[:int((T/2-.02)/dt)])
+    close('all')
+    if any(frateHits[int(T/2/dt)-legroom1:int(T/2/dt)+legroom2]>e0+nsigmaON*sde):	
+        BratesON_Hits[counton,:]=brateHits*1.0
+        FratesON_Hits[counton,:]=frateHits*1.0
+        nON_Hits += nHits[ind][0][0]
+        nON_Missed += nHits[ind][0][0]
+
+        cellno=counton
+        oninds = append(oninds,i)
+        counton +=1
+    elif any(frateHits[int(T/2/dt)-legroom1:int(T/2/dt)+legroom2]<e0-nsigmaOFF*sde): 
+        BratesOFF_Hits[countoff,:]=brateHits*1.0
+        FratesOFF_Hits[countoff,:]=frateHits*1.0
+        nOFF_Hits += nHits[ind][0][0]
+        nOFF_Missed += nHits[ind][0][0]
+        cellno=countoff
+        offinds= append(offinds,i)
+        countoff +=1
+    else:
+        BratesNR_Hits[countnr,:]=brateHits*1.0
+        FratesNR_Hits[countnr,:]=frateHits*1.0
+        nNR_Hits += nHits[ind][0][0]
+        nNR_Missed += nHits[ind][0][0]
+        countnr +=1
+        
+
+# Plotting
+close()
+fig1,ax1=create_1tracefig((-.1,6))
+mu = mean(BratesON_Hits[:counton,:],axis=0)
+ax1.plot(t-T/2.0,mu/mean(mu[:20000]),'orange')
+mu = mean(FratesON_Hits[:counton,:],axis=0)
+ax1.plot(t-T/2.0,mu/mean(mu[:20000]),'green')
+mu = mean(FratesOFF_Hits[:countoff,:],axis=0)
+ax1.plot(t-T/2.0,mu/mean(mu[:19000]),'cyan')
+mu = mean(FratesNR_Hits[:countnr,:],axis=0)
+ax1.plot(t-T/2.0,mu/mean(mu[:19000]),'gray')
+ax1.set_xlim((-1,2))
+ax1.set_xticks([-1,0,1,2])
+ax1.set_ylim((0,6 ))
+ax1.set_ylabel('Norm. Rate')
+fig1.savefig(task+'StimTrigTrialAvgRates_ONOFF_Normalized.pdf',format='pdf')
+
+close()
+fig1,ax1,ax2=create_2tracefig((-.1,6))
+mu = mean(BratesON_Hits[:counton,:],axis=0)
+ax2.plot(t-T/2.0,mu,'red')
+mu = mean(FratesON_Hits[:counton,:],axis=0)
+ax1.plot(t-T/2.0,mu,'green')
+mu = mean(FratesOFF_Hits[:countoff,:],axis=0)
+ax1.plot(t-T/2.0,mu,'cyan')
+mu = mean(FratesNR_Hits[:countnr,:],axis=0)
+ax1.plot(t-T/2.0,mu,'gray')
+ax1.set_xlim((-1,2))
+ax1.set_xticks([-1,0,1,2])
+ax1.set_ylim((0,34 ))
+ax2.set_ylim((0,4))
+ax2.set_ylabel('Burst Rate [Hz]')
+ax1.set_ylabel('Firing Rate [Hz]')
+fig1.savefig(task+'StimTrigTrialAvgRates_ONOFF_Dual.pdf',format='pdf')
+
+# Saving
+downsampletime = arange(0,FratesON_Hits.shape[1],smboxsize-1)
+time = downsampletime*dt*1e3 # time in ms
+io.savemat(task+'PSTHs_3Classes_3sigOn03to04s',{'time':time,'BratesON_Hits':BratesON_Hits[:counton,downsampletime],'BratesNR_Hits':BratesNR_Hits[:countnr,downsampletime],'BratesOFF_Hits':BratesOFF_Hits[:countoff,downsampletime],'FratesON_Hits':FratesON_Hits[:counton,downsampletime],'FratesNR_Hits':FratesNR_Hits[:countnr,downsampletime],'FratesOFF_Hits':FratesOFF_Hits[:countoff,downsampletime]})
+
+
+

codes/Electrophysiology analysis/FindSigModulatedBins.m

@@ -0,0 +1,218 @@
+%%
+load('PSTH data.mat', 'Control_PSTHz','Control_bPSTHz', 'Dreadd_PSTHz','Dreadd_bPSTHz', 'Control_sorted', 'Dreadd_sorted')
+%%
+
+% bin comparison with bonferroni correction
+% any cell having at least 2 significant bins
+binNum = 46;
+BonferroniCorrection = 0.05/binNum
+threshold = abs(norminv(BonferroniCorrection/2))
+%%
+Sorted_Control_PSTHz = Control_PSTHz(:,Control_sorted(1:37));
+Sorted_Control_bPSTHz = Control_bPSTHz(:,Control_sorted(1:37));
+
+CtrlFreqSigBins = cell(37,1);
+
+CtrlBurstSigNum = zeros(37,1);
+
+for i = 1:37
+ PostPSTHz = abs(Sorted_Control_PSTHz(25:70,i));
+ CtrlFreqSigBins{i} = find(PostPSTHz > threshold);
+ 
+ 
+ 
+ PostBPSTHz = abs(Sorted_Control_bPSTHz(25:70,i));
+ CtrlBurstSigNum (i) = sum(PostBPSTHz > threshold);
+end
+
+CtrlFreqSigNum = cellfun(@length, CtrlFreqSigBins);
+
+
+CtrlAvgZscoreSigBin = zeros(37,1);
+for i = 1:37
+    CtrlAvgZscoreSigBin(i) = mean(Control_PSTHz(CtrlFreqSigBins{i},Control_sorted(i)));
+end
+
+%ControlBurstNull = find(sum(Control_bPSTHz(1:20,Control_sorted(1:37))==0));
+%CtrlBurstSigNum(find(Control_sorted==ControlBurstNull)) = NaN;
+%%
+Sorted_Dreadd_PSTHz = Dreadd_PSTHz(:,Dreadd_sorted(1:36));
+Sorted_Dreadd_bPSTHz = Dreadd_bPSTHz(:,Dreadd_sorted(1:36));
+
+ExpFreqSigBins = cell(36,1);
+
+ExpBurstSigNum = zeros(36,1);
+
+for i = 1:36
+ PostPSTHz = abs(Dreadd_PSTHz(25:70,i));
+ ExpFreqSigBins{i} = find(PostPSTHz > threshold);
+ 
+ 
+ 
+ PostBPSTHz = abs(Dreadd_bPSTHz(25:70,i));
+ ExpBurstSigNum (i) = sum(PostBPSTHz > threshold);
+end
+
+ExpFreqSigNum = cellfun(@length, ExpFreqSigBins);
+
+
+ExpAvgZscoreSigBin = zeros(36,1);
+for i = 1:36
+    ExpAvgZscoreSigBin(i) = mean(Dreadd_PSTHz(ExpFreqSigBins{i},Dreadd_sorted(i)));
+end
+
+%DreaddBurstNull = find(sum(Dreadd_bPSTHz(1:20,Dreadd_sorted(1:37))==0));
+%ExpBurstSigNum(DreaddBurstNull) = NaN;
+
+%%
+
+figure
+subplot(211)
+[~,sorted] = sort(CtrlFreqSigNum);
+plot(1:37,CtrlFreqSigNum(sorted),'ko')
+%ylim([0 20])
+xlim([0 37])
+hold on
+plot([0 40], [1 1], 'k--')
+title('Firing rate')
+ylabel('Number of significant bins')
+
+subplot(212)
+plot(1:37,CtrlBurstSigNum(sorted),'ko')
+%ylim([0 20])
+xlim([0 37])
+hold on
+plot([0 40], [1 1], 'k--')
+title('Burst rate')
+ylabel('Number of significant bins')
+
+%%
+figure
+subplot(211)
+[~,sorted] = sort(ExpFreqSigNum);
+plot(1:36,ExpFreqSigNum(sorted),'ro')
+%ylim([0 6])
+xlim([0 36])
+hold on
+plot([0 40], [1 1], 'k--')
+title('Firing rate')
+ylabel('Number of significant bins')
+xlabel('Cell Number')
+
+subplot(212)
+plot(1:36,ExpBurstSigNum(sorted),'ro')
+%ylim([0 6])
+xlim([0 36])
+hold on
+plot([0 40], [1 1], 'k--')
+title('Burst rate')
+ylabel('Number of significant bins')
+xlabel('Cell Number')
+%%
+
+CtrlModulated = find(CtrlFreqSigNum > 0);
+ExpModulated = find(ExpFreqSigNum > 0);
+
+figure
+bar([length(CtrlModulated)*100/length(CtrlFreqSigBins), (length(CtrlFreqSigBins)-length(CtrlModulated))*100/length(CtrlFreqSigBins); ...,
+    length(ExpModulated)*100/length(ExpFreqSigNum), (length(ExpFreqSigNum)-length(ExpModulated))*100/length(ExpFreqSigNum)],'stacked')
+set(gca, 'xticklabel', {'Control (n=37)', 'hM4Di (n=36)'})
+ylabel('Cell num')
+legend({'Modulated', 'Unmodulated'})
+title('Firing rate')
+
+%figure
+% subplot(121)
+% pie([length(CtrlModulated), length(CtrlFreqSigNum)-length(CtrlModulated)])
+% title('Control (n=37)')
+% %legend({'Modulated', 'Unmodulated'})
+% 
+% subplot(122)
+% pie([length(ExpModulated), length(ExpFreqSigNum)-length(ExpModulated)])
+% title('hM4Di (n=36)')
+% %legend({'Modulated', 'Unmodulated'})
+% 
+n1 = length(CtrlModulated); N1 = length(CtrlFreqSigBins);
+n2 = length(ExpModulated); N2 = length(ExpFreqSigNum);
+x1 = [repmat('a',N1,1); repmat('b',N2,1)];
+x2 = [repmat(1,n1,1); repmat(2,N1-n1,1); repmat(1,n2,1); repmat(2,N2-n2,1)];
+[tbl,chi2stat,pval] = crosstab(x1,x2)
+
+
+%%
+CtrlPosModulated = find(CtrlAvgZscoreSigBin> 0);
+CtrlNegModulated = find(CtrlAvgZscoreSigBin < 0);
+
+ExpPosModulated = find(ExpAvgZscoreSigBin> 0);
+ExpNegModulated = find(ExpAvgZscoreSigBin < 0);
+
+
+figure
+bar([length(CtrlPosModulated)*100/length(CtrlFreqSigBins), length(CtrlNegModulated)*100/length(CtrlFreqSigBins), ...,
+    (length(CtrlFreqSigBins)-(length(CtrlPosModulated)+length(CtrlNegModulated)))*100/length(CtrlFreqSigBins); 
+length(ExpPosModulated)*100/length(ExpFreqSigNum), length(ExpNegModulated)*100/length(ExpFreqSigNum), ...,
+(length(ExpFreqSigBins)-(length(ExpPosModulated)+length(ExpNegModulated)))*100/length(ExpFreqSigBins)],'stacked')
+set(gca, 'xticklabel', {'Control (n=37)', 'hM4Di (n=36)'})
+ylabel('Cell num')
+legend({'Postively Modulated', 'Negatively Modulated'})
+title('Firing rate')
+
+
+%%
+figure
+PlotAllDataPoints(CtrlFreqSigNum, ExpFreqSigNum)
+set(gca,'xticklabel', {'Control', 'hM4Di'})
+
+p = ranksum (CtrlFreqSigNum, ExpFreqSigNum)
+text(1.2, 30, ['p = ', num2str(p)])
+
+
+%%
+CtrlModulated = find(CtrlBurstSigNum > 0);
+ExpModulated = find(ExpBurstSigNum > 0);
+
+CtrlBurstCellNum = length(CtrlBurstSigNum)
+ExpBurstCellNum = length(ExpBurstSigNum)
+
+figure
+bar([length(CtrlModulated)*100/CtrlBurstCellNum, (CtrlBurstCellNum-length(CtrlModulated))*100/CtrlBurstCellNum; length(ExpModulated)*100/ExpBurstCellNum, (ExpBurstCellNum-length(ExpModulated))*100/ExpBurstCellNum],'stacked')
+set(gca, 'xticklabel', {'Control (n=37)', 'hM4Di (n=36)'})
+ylabel('Cell num')
+legend({'Modulated', 'Unmodulated'})
+title('Burst rate')
+
+% figure
+% subplot(121)
+% pie([length(CtrlModulated), CtrlBurstCellNum-length(CtrlModulated)])
+% title(['Control (n=', num2str(CtrlBurstCellNum), ')'])
+% %legend({'Modulated', 'Unmodulated'})
+% 
+% subplot(122)
+% pie([length(ExpModulated), ExpBurstCellNum-length(ExpModulated)])
+% title(['hM4Di (n=', num2str(ExpBurstCellNum), ')'])
+% %legend({'Modulated', 'Unmodulated'})
+
+n1 = length(CtrlModulated); N1 = CtrlBurstCellNum;
+n2 = length(ExpModulated); N2 = ExpBurstCellNum;
+x1 = [repmat('a',N1,1); repmat('b',N2,1)];
+x2 = [repmat(1,n1,1); repmat(2,N1-n1,1); repmat(1,n2,1); repmat(2,N2-n2,1)];
+[tbl,chi2stat,pval] = crosstab(x1,x2)
+
+%%
+figure
+PlotAllDataPoints(CtrlBurstSigNum, ExpBurstSigNum)
+set(gca,'xticklabel', {'Control', 'hM4Di'})
+
+p = ranksum (CtrlBurstSigNum, ExpBurstSigNum)
+text(1.2, 25, ['p = ', num2str(p)])
+
+
+%%
+CtrlM = Control_PSTHz(:,Control_sorted(CtrlModulated)); %Modulated cells
+ExpM = Dreadd_PSTHz(:,Dreadd_sorted(ExpModulated));
+
+figure
+boundedline(edges(1:end-1), mean(CtrlM,2), std(CtrlM,[],2)/sqrt(length(CtrlModulated)),'k') 
+hold on
+boundedline(edges(1:end-1), mean(ExpM,2), std(ExpM,[],2)/sqrt(length(ExpModulated)),'r') 
+

codes/Electrophysiology analysis/GetAllCellsPSTHnorm.m

@@ -0,0 +1,41 @@
+function [CellStruct, AllCellsPSTHnormHit] = GetAllCellsPSTHnorm(CellStruct,preStim,postStim, binsize)
+% Get PSTHs from RS cells with at least 5 trials
+
+binNum = (postStim + preStim)/binsize;
+CellNum = length(CellStruct);
+
+AllCellsPSTHnormHit = zeros(binNum,CellNum);
+
+for i = 1:length(CellStruct)
+    if CellStruct(i).CellType == 1 %Only RS cells
+        SpikeTimes = CellStruct(i).SpikeTimes;
+        
+        % Add mix trials
+        if length(CellStruct(i).Mixes) ~=0
+            Mixes = CellStruct(i).Mixes;
+            Mixes(:,3) = 10;
+            Microstims = [CellStruct(i).Microstims; Mixes];
+        else
+            Microstims = [CellStruct(i).Microstims];
+        end
+        
+        
+        % Hit trials
+        HitIdx = find(Microstims(:,2) <= 1.2 & Microstims (:,2) > 0.1);
+        HitTimes = Microstims(HitIdx,1);
+        if length(HitTimes) > 4 % at least 5 trials
+            [psthNormHit,numTrials, edges] = RasterPSTH (HitTimes, SpikeTimes, preStim, postStim, binsize,'mkplt', false );
+            CellStruct(i).psthNorm = psthNormHit;
+            CellStruct(i).HitNum = numTrials;
+            AllCellsPSTHnormHit(:,i) = psthNormHit;
+        end
+        
+%         % Miss trials
+%         MissIdx = find(Microstims(:,2) ==-0.1);
+%         MissTimes = Microstims(MissIdx,1);
+%         if length(MissTimes) > 4 % at least 5 trials
+%         [psthNormMiss,numTrials, edges] = RasterPSTH (MissTimes, SpikeTimes, preStim, postStim, binsize,'mkplt', false );
+%         AllCellsPSTHnormMiss(:,i) = psthNormMiss;
+%         end
+    end
+end

codes/Electrophysiology analysis/PlotCellAverage.m

@@ -0,0 +1,231 @@
+load('PooledCellsData.mat')
+%% Connected dots - Hit trials
+figure
+
+subplot(121)
+ConnectedDots2 (Control_hit_freq(:,1), Control_hit_freq(:,2)) % Unclear_hit_freq(:,1)])
+%set(gca,'XTickLabel',{'Control (n=41)','hM4Di (n=40)'})
+ylabel('Spike rate (Hz)')
+title('Control')
+[p h] = signrank (Control_hit_freq(:,1), Control_hit_freq(:,2)) % Unclear_hit_freq(:,1)])
+hold on
+text(1.2, 20, ['p = ', num2str(p)]) 
+
+subplot(122)
+ConnectedDots2 (newDREADD_hit_freq(:,1), newDREADD_hit_freq(:,2)) % Unclear_hit_freq(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Spike rate (Hz)')
+title('hM4Di')
+[p h] = signrank (newDREADD_hit_freq(:,1), newDREADD_hit_freq(:,2)) % Unclear_hit_freq(:,1)])
+hold on
+text(1.2, 20, ['p = ', num2str(p)])  
+ylim([0 25])
+
+figure
+
+subplot(121)
+ConnectedDots2 (Control_hit_burst(:,1), Control_hit_burst(:,2)) % Unclear_hit_burst(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Burst rate (Hz)')
+title('Control')
+[p h] = signrank (Control_hit_burst(:,1), Control_hit_burst(:,2)) % Unclear_hit_burst(:,1)])
+hold on
+text(1.2, 2, ['p = ', num2str(p)]) 
+
+subplot(122)
+ConnectedDots2 (newDREADD_hit_burst(:,1), newDREADD_hit_burst(:,2)) % Unclear_hit_burst(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Burst rate (Hz)')
+title('hM4Di')
+[p h] = signrank (newDREADD_hit_burst(:,1), newDREADD_hit_burst(:,2)) % Unclear_hit_burst(:,1)])
+hold on
+text(1.2, 2, ['p = ', num2str(p)])  
+ylim([0 2.5])
+
+%% Connected Dots - miss trials
+figure
+
+subplot(121)
+ConnectedDots2 (Control_miss_freq(:,1), Control_miss_freq(:,2)) % Unclear_miss_freq(:,1)])
+%set(gca,'XTickLabel',{'Control (n=41)','hM4Di (n=40)'})
+ylabel('Spike rate (Hz)')
+title('Control')
+ylim ([0 25])
+[p h] = signrank (Control_miss_freq(:,1), Control_miss_freq(:,2)) % Unclear_miss_freq(:,1)])
+hold on
+text(1.2, 20, ['p = ', num2str(p)]) 
+
+subplot(122)
+ConnectedDots2 (newDREADD_miss_freq(:,1), newDREADD_miss_freq(:,2)) % Unclear_miss_freq(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Spike rate (Hz)')
+title('hM4Di')
+[p h] = signrank (newDREADD_miss_freq(:,1), newDREADD_miss_freq(:,2)) % Unclear_miss_freq(:,1)])
+hold on
+text(1.2, 20, ['p = ', num2str(p)])  
+ylim([0 25])
+
+figure
+
+subplot(121)
+ConnectedDots2 (Control_miss_burst(:,1), Control_miss_burst(:,2)) % Unclear_miss_burst(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Burst rate (Hz)')
+title('Control')
+ylim([0 2.5])
+[p h] = signrank (Control_miss_burst(:,1), Control_miss_burst(:,2)) % Unclear_miss_burst(:,1)])
+hold on
+text(1.2, 2, ['p = ', num2str(p)]) 
+
+subplot(122)
+ConnectedDots2 (newDREADD_miss_burst(:,1), newDREADD_miss_burst(:,2)) % Unclear_miss_burst(:,1)])
+%set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Burst rate (Hz)')
+title('hM4Di')
+[p h] = signrank (newDREADD_miss_burst(:,1), newDREADD_miss_burst(:,2)) % Unclear_miss_burst(:,1)])
+hold on
+text(1.2, 2, ['p = ', num2str(p)])  
+ylim([0 2.5])
+
+
+%% Difference
+
+Control_hit_freq(:,3) = Control_hit_freq(:,2) - Control_hit_freq(:,1);
+newDREADD_hit_freq(:,3) = newDREADD_hit_freq(:,2) - newDREADD_hit_freq(:,1);
+Unclear_hit_freq(:,3) =Unclear_hit_freq(:,2) - Unclear_hit_freq(:,1);
+
+figure
+subplot(121)
+PlotAllDataPoints (Control_hit_freq(:,3), newDREADD_hit_freq(:,3) )
+set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Post - Pre')
+title('Spike rate change')
+
+[p h] = ranksum (Control_hit_freq(:,3), newDREADD_hit_freq(:,3))
+
+hold on
+text(1.2, 10, ['p = ', num2str(p)]) 
+
+subplot(122)
+PlotAllDataPoints (abs(Control_hit_freq(:,3)), abs(newDREADD_hit_freq(:,3)) )
+set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('|Post - Pre|')
+title('Spike rate Absolute change')
+
+[p h] = ranksum (abs(Control_hit_freq(:,3)), abs(newDREADD_hit_freq(:,3)))
+
+hold on
+text(1.2, 10, ['p = ', num2str(p)])
+
+% Bursts
+Control_hit_burst(:,3) = Control_hit_burst(:,2) - Control_hit_burst(:,1);
+newDREADD_hit_burst(:,3) = newDREADD_hit_burst(:,2) - newDREADD_hit_burst(:,1);
+Unclear_hit_burst(:,3) =Unclear_hit_burst(:,2) - Unclear_hit_burst(:,1);
+
+
+figure
+subplot(121)
+PlotAllDataPoints (Control_hit_burst(:,3), newDREADD_hit_burst(:,3) )
+set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Post - Pre')
+title('Burst rate change')
+
+[p h] = ranksum (Control_hit_burst(:,3), newDREADD_hit_burst(:,3))
+
+hold on
+text(1.2, 1, ['p = ', num2str(p)]) 
+
+
+subplot(122)
+PlotAllDataPoints (abs(Control_hit_burst(:,3)), abs(newDREADD_hit_burst(:,3)) )
+set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('|Post - Pre|')
+title('Spike rate Absolute change')
+
+[p h] = ranksum (abs(Control_hit_burst(:,3)), abs(newDREADD_hit_burst(:,3)))
+
+hold on
+text(1.2, 1, ['p = ', num2str(p)]) 
+%% Relative change
+Control_baseline = mean(Control_hit_burst(:,1));
+newDREADD_baseline = mean(newDREADD_hit_burst(:,1));
+
+Control_hit_burst(:,4) = Control_hit_burst(:,3)./Control_baseline;
+newDREADD_hit_burst(:,4) = newDREADD_hit_burst(:,3)./newDREADD_baseline;
+
+[p h] = ranksum (Control_hit_burst(:,4), newDREADD_hit_burst(:,4))
+
+figure
+PlotAllDataPoints (Control_hit_burst(:,4).*100, newDREADD_hit_burst(:,4).*100 )
+set(gca,'XTickLabel',{'Control (n=42)','hM4Di (n=26)'})
+ylabel('Relative change (%)')
+title('Burst')
+
+
+%% Post stim
+figure
+
+subplot(121)
+PlotAllDataPoints (Control_hit_freq(:,2), newDREADD_hit_freq(:,2) )
+set(gca,'XTickLabel',{'Control (n=41)','hM4Di (n=41)'})
+ylabel('Spike rate (Hz)')
+title('Hit - Post stimulation')
+
+[p h] = ranksum (Control_hit_freq(:,2), newDREADD_hit_freq(:,2))
+
+hold on
+text(1.2, 10, ['p = ', num2str(p)]) 
+
+subplot(122)
+PlotAllDataPoints (Control_miss_freq(:,2), newDREADD_miss_freq(:,2) )
+set(gca,'XTickLabel',{'Control (n=18)','hM4Di (n=17)'})
+ylabel('Spike rate (Hz)')
+title('Miss - Post stimulation')
+
+[p h] = ranksum (Control_miss_freq(:,2), newDREADD_miss_freq(:,2))
+
+hold on
+text(1.2, 8, ['p = ', num2str(p)]) 
+ylim([0 25])
+
+%% Bar wit error
+N1 = size(Control_hit_freq,1)
+Control_Avg_hitFreq = nanmean(Control_hit_freq);
+Control_SEM_hitFreq = nanstd(Control_hit_freq,1)./sqrt(N1);
+Control_Avg_hitBurst = nanmean(Control_hit_burst)  ;
+Control_SEM_hitBurst = nanstd(Control_hit_burst,1)./sqrt(N1);
+
+N2 = size(newDREADD_hit_freq,1)
+newDREADD_Avg_hitFreq = nanmean(newDREADD_hit_freq);
+newDREADD_SEM_hitFreq = nanstd(newDREADD_hit_freq,1)./sqrt(N2);
+newDREADD_Avg_hitBurst = nanmean(newDREADD_hit_burst)  ;
+newDREADD_SEM_hitBurst = nanstd(newDREADD_hit_burst,1)./sqrt(N2);
+
+
+% %% Plotting
+figure
+subplot(2,1,1)
+varargout = barwitherr([Control_SEM_hitFreq; newDREADD_SEM_hitFreq],[Control_Avg_hitFreq; newDREADD_Avg_hitFreq])  
+set(gca,'XTickLabel',{'Control (n=41)','hM4Di (n=40)'})
+ylabel('Firing rate (Hz)')
+ylim([0 5])
+
+title('hit trials')
+legend({'preStim', 'postStim'})
+
+subplot(2,1,2)
+varargout = barwitherr([Control_SEM_hitBurst; newDREADD_SEM_hitBurst],[Control_Avg_hitBurst;newDREADD_Avg_hitBurst])  
+set(gca,'XTickLabel',{'Control (n=41)','hM4Di (n=40)'})
+ylabel('Burst rate (Hz)')
+
+
+%% Violin plot
+figure
+violin(Control_hit_burst);
+set(gca, 'XTick',[1,2],'XTickLabel',{'PreStim', 'PostStim'})
+title('Control')
+
+figure
+violin(newDREADD_hit_burst);
+set(gca, 'XTick',[1,2],'XTickLabel',{'PreStim', 'PostStim'})
+title('hM4Di')

codes/Electrophysiology analysis/PreparePSTHs.m

@@ -0,0 +1,14 @@
+function [PSTH, bPSTH] = PreparePSTHs(Group, preStim,postStim,binsize,mkplt)
+
+
+
+
+Group = AllCellsBurstPSTH(Group,preStim,postStim, binsize, mkplt);
+[~, PSTH]  = GetAllCellsPSTHnorm(Group,preStim,postStim, binsize);
+
+numHits = [Group.numHits];
+trialthreshold = find(numHits>4);
+RS = find([Group.CellType] == 1);
+add = intersect(trialthreshold, RS);
+PSTH = PSTH(:,add);
+bPSTH = [Group(add).bpsthNorm];

codes/Histology analysis/Anterograde tracing/NormalizeTraces.m

@@ -0,0 +1,20 @@
+function norm_total = NormalizeTraces (totalYFP)
+
+% detrendYFP = cell(size(totalYFP,2),1);
+% for i = 1:size(totalYFP,2)
+%  lastidx= find(isnan(totalYFP(:,i))==1,1) -1;
+%  if isempty(lastidx) == 1
+%     detrendYFP{i} = detrend(totalYFP(:,i));
+%  else
+%      detrendYFP{i} = detrend(totalYFP(1:lastidx,i));
+%  end
+%  
+% end
+
+
+baseline = min(totalYFP,[],1);
+
+normfactor = max(totalYFP-baseline);
+
+norm_total = (totalYFP-baseline)./normfactor;
+

codes/Histology analysis/Anterograde tracing/PlotYFPQuantification.m

@@ -0,0 +1,43 @@
+figure
+plot(depth,PRy03)
+hold on
+plot(depth,mean(PRy03,2),'k','linew',2)
+xlim([0 1200])
+%%
+%norm_PRy03 = bsxfun(@rdivide,PRy03, max(PRy03,[],1));
+%norm_PRy02 = bsxfun(@rdivide,PRy02, max(PRy02,[],1));
+
+RatYFP = [PRy02 PRy03];
+
+
+
+%norm_total = bsxfun(@rdivide, total, max(total,[],1));
+%% Normalization
+norm_RatYFP = NormalizeTraces (RatYFP);
+
+% Mouse data start from 50 um
+startidx = dsearchn(MouseDepth,50);
+MouseYFPnew = MouseYFP;
+MouseYFPnew(1:startidx,:) = []; 
+%MouseDepthNew = MouseDepth(startidx:end);
+norm_MouseYFP = NormalizeTraces (MouseYFPnew);
+
+
+%%
+figure
+plot(RatDepth,norm_RatYFP)
+hold on
+boundedline(RatDepth, nanmean(norm_RatYFP,2), nanstd(norm_RatYFP,[],2)./sqrt(size(norm_RatYFP,2)),'k-') 
+xlim([0 1200])
+xlabel('Depth (um)')
+ylabel('Normalized intensity')
+
+
+%%
+figure
+plot(MouseDepth,norm_MouseYFP)
+hold on
+boundedline(MouseDepth, nanmean(norm_MouseYFP,2), nanstd(norm_MouseYFP,[],2)./sqrt(size(norm_MouseYFP,2)),'k-') 
+xlim([50 1000])
+xlabel('Depth (um)')
+ylabel('Normalized intensity')

codes/Histology analysis/DREADD effect quantification/DreaddEffectCorticalDepth.m

@@ -0,0 +1,73 @@
+load ('CBdata.mat', 'norm_CB', 'depth')
+load('PRh axon data.mat', 'norm_MouseYFP','MouseDepth')
+%%
+samedepth = dsearchn(MouseDepth, depth);
+
+
+%% Averaging two bins in MouseYFP
+temp1 = norm_MouseYFP(1:2:end,:);
+temp2 = norm_MouseYFP(2:2:end,:);
+
+norm_MouseYFPcorr = (temp1+temp2)./2;
+norm_MouseYFPcorr(387,:) = norm_MouseYFP(772,:);
+
+%% Effect size
+effect = nanmean(norm_MouseYFPcorr,2) .* nanmean(norm_CB,2);
+
+%%
+figure
+
+boundedline(depth, nanmean(norm_MouseYFPcorr,2), nanstd(norm_MouseYFPcorr,[],2)./sqrt(size(norm_MouseYFPcorr,2)),'g-') 
+
+hold on
+boundedline(depth, mean(norm_CB,2), std(norm_CB,[],2)./sqrt(size(norm_CB,2)),'r-') 
+plot(depth,effect,'k','linew',2)
+xlim([50 1000])
+ylim([0 1])
+xlabel('Depth (um)')
+ylabel('Normmalized Amplitude')
+
+
+% Layer boundaries (Lefort et al., 2009, Neuron)
+L1 = 128;
+L23 = 418;
+L4 = 588;
+L5 = 890;
+L6 = 1154;
+yvals = get(gca,'Ylim');
+
+hold on
+
+plot([L1 L1], yvals, 'r--')
+plot([L23 L23], yvals, 'r--')
+plot([L4 L4], yvals, 'r--')
+plot([L5 L5], yvals, 'r--')
+
+
+
+
+legend({'YFP', 'Chicago blue', 'Effect'})
+%% Quantify Area under the curve
+interval = depth(2)-depth(1);
+area = effect*interval;
+startidx = dsearchn(depth,50)+1; %MouseYFP starts from 50 um
+TotalArea = sum(area(startidx:end));
+
+LayerBoundaries = [L1; L23; L4; L5];
+LayerIdx = zeros(4,1);
+for i = 1:4
+    LayerIdx(i) = dsearchn(depth, LayerBoundaries(i));
+end
+
+
+LayerIdx = [startidx-1; LayerIdx; 387];
+
+LayerAreas = zeros(5,1);
+for i = 1:length(LayerIdx)-1
+    LayerAreas(i) = sum(area(LayerIdx(i)+1:LayerIdx(i+1)));
+end
+
+Percentage = LayerAreas./TotalArea;
+Percentage = Percentage*100
+
+

codes/Histology analysis/hM4Di quantification/PlotAvgNormDreaddExpression.m

@@ -0,0 +1,50 @@
+load('DreaddExpressionData.mat')
+%%
+AvgJS09norm = nanmean(JS09norm,2);
+
+%%
+rf = find(JS08(:,1) == 0);
+AvgJS09norm = AvgJS08norm(rf-386:rf+386);
+
+distance = JS08(rf-386:rf+386,1);
+
+%%
+TotalAvgNorm = [AvgR002norm, AvgR003norm, AvgR004norm, AvgY754norm, AvgY755norm, AvgJS09norm, ...,
+    AvgJS12norm, AvgJS13norm, AvgJS14norm, AvgJS22norm; AvgJS08norm];
+%%
+plot(distance, TotalAvgNorm)
+xlim([-1000 1000])
+hold on
+plot(distance, nanmean(TotalAvgNorm,2),'k', 'linew',2)
+yval = get(gca,'ylim')
+plot([-250 -250], yval, 'r--')
+plot([750 750], yval, 'r--')
+
+
+%% Quantification
+interval = distance(2)-distance(1);
+area = TotalAvgNorm.*interval;
+
+TotalArea = nansum(area,1);
+
+PRhBoundaries = [-250; 750];
+PRhUpperIdx = dsearchn(distance, PRhBoundaries(1));
+PRhLowerIdx = dsearchn(distance, PRhBoundaries(2));
+
+PRhArea = nansum(area(PRhUpperIdx:PRhLowerIdx,:));
+PRhAbove = nansum(area(1:PRhUpperIdx-1,:));
+PRhBelow = nansum(area(PRhLowerIdx+1:end,:));
+
+PRhPerc = 100*PRhArea./TotalArea
+PRhAbovePerc = 100*PRhAbove./TotalArea
+PRhBelowPerc = 100*PRhBelow./TotalArea
+
+%Specificity = PRhPerc./NonPRhPerc
+
+
+%%
+group = [ones(11,1); ones(11,1)+1; ones(11,1)+2];
+boxplot([PRhAbovePerc';PRhPerc';PRhBelowPerc'], group)
+ylim([0 100])
+ylabel('Expression (%)')
+set(gca,'xticklabel', {'Above PRh', 'PRh', 'Below PRh'})

codes/Nanostimulation analysis/Fig4Science_data_example_cell.m

@@ -0,0 +1,223 @@
+function Fig4Science_data_example_cell
+
+
+close all
+Fig4 = figure('Name','Fig4 Science Rebuttal 033d110511s3', 'PaperUnits', 'centimeters');
+Width = 12;
+Height = 17;
+PixelsPerCm = 85/2.54;
+set(gcf, 'PaperSize', [Width Height])
+set(gcf, 'PaperPosition', [.5 1 Width Height])
+set(gcf, 'Position', get(gcf, 'PaperPosition') * PixelsPerCm);
+set(gcf,'DefaultAxesFontSize',10, 'DefaultTextFontName', 'Helvetica','DefaultAxesFontName', 'Helvetica')
+
+PanelLetters = cell(4,1);
+PanelLetters{1,1} = 'b';
+Dimension = size(PanelLetters);
+RightMargin = .2;
+LeftMargin = .07;
+TopMargin = .07;
+BottomMargin = .05;
+InterPanelWidth = eps;
+InterPanelHeight = .05;
+LetterFontSize = 12;
+LetterX = -.05;
+LetterY = -.035;
+PanelHandles = MultiPanel(Fig4,Dimension, [LeftMargin, RightMargin, BottomMargin, TopMargin], [InterPanelWidth, InterPanelHeight], [LetterX, LetterY], PanelLetters, LetterFontSize);
+%
+PanelLetters = cell(4,1);
+PanelLetters{1,1} = 'c';
+Dimension = size(PanelLetters);
+RightMargin = .01;
+LeftMargin = .9;
+PanelHandles(:, end + 1) = MultiPanel(Fig4,Dimension, [LeftMargin, RightMargin, BottomMargin, TopMargin], [InterPanelWidth, InterPanelHeight], [LetterX, LetterY], PanelLetters, LetterFontSize);
+
+load('Fig4Science_data_example_cell_vars.mat');
+
+axes(PanelHandles(1,1))
+slen= 0.5;  %length of line representing spike
+
+HitLineWidth=3;
+
+TrialNum=length(Nano25ConIncTimes);
+TrialLim=TrialNum+0.5;
+for i= 1:length(Nano25ConIncTimes)
+   t= Nano25ConIncTimes(i,1);
+   rt= Nano25ConIncTimes(i,2);
+   s= SpikeTimesNew(SpikeTimesNew>t-1 & SpikeTimesNew<t+2);
+   for j= 1:length(s)
+       if (s(j)<t | s(j)>t+0.001) & (s(j)<t+0.199 | s(j)>t+0.2)
+           l= line([s(j)-t s(j)-t],[i-slen i+slen]);
+           set(l,'color','k')
+       end
+   end
+   if rt>0 & rt<2
+       l= line([rt rt],[i-slen i+slen]);
+       set(l,'color','r')
+       set(l,'linewidth',HitLineWidth)
+   end
+end
+axis([-1 2 0 TrialLim])
+set(gca,'xticklabel',[])
+set(gca,'ytick',[0 TrialNum])
+set(gca,'xtick',[-1 0 1])
+set(gca,'tickdir','out','FontName','Helvetica')
+myxlim=xlim;
+myylim=ylim;
+
+axes(PanelHandles(1,2))
+p= bar([Nano25IncHits]);
+set(p,'facecolor','none')
+
+ylim([0 (Nano25IncHits+Nano25IncMisses)])
+xlim([.5 1.5])
+set(gca,'xtick',[]);
+set(gca,'ytick',[0 Nano25IncHits+Nano25IncMisses],'TickLength',[0.05 0.05])
+set(gca,'ytick',[0 Nano25IncHits+Nano25IncMisses])
+set(gca,'xticklabel',[])
+set(gca,'tickdir','out','FontName','Helvetica')
+
+myxlim=xlim;
+myylim=ylim;
+text(myxlim(1)+diff(myxlim)*0.3,myylim(2)-diff(myylim)*0.05,[num2str(round(Nano25HitRate*100)), '%'],'FontName','Helvetica');
+box off
+
+
+axes(PanelHandles(2,1))
+slen= 0.5;  %length of line representing spike
+
+HitLineWidth=3;
+
+TrialNum=length(Nano100ConIncTimes);
+TrialLim=TrialNum+0.5;
+for i= 1:length(Nano100ConIncTimes)
+   t= Nano100ConIncTimes(i,1);
+   rt= Nano100ConIncTimes(i,2);
+   s= SpikeTimesNew(SpikeTimesNew>t-1 & SpikeTimesNew<t+2);
+   for j= 1:length(s)
+       if (s(j)<t | s(j)>t+0.001) & (s(j)<t+0.199 | s(j)>t+0.2)
+           l= line([s(j)-t s(j)-t],[i-slen i+slen]);
+           set(l,'color','k')
+       end
+   end
+   if rt>0 & rt<2
+       l= line([rt rt],[i-slen i+slen]);
+       set(l,'color','r')
+       set(l,'linewidth',HitLineWidth)
+   end
+end
+axis([-1 2 0 TrialLim])
+set(gca,'xticklabel',[])
+set(gca,'ytick',[0 TrialNum])
+set(gca,'xtick',[-1 0 1])
+set(gca,'tickdir','out','FontName','Helvetica')
+myxlim=xlim;
+myylim=ylim;
+
+
+axes(PanelHandles(2,2))
+p= bar([Nano100IncHits]);
+set(p,'facecolor','none')
+ylim([0 (Nano100IncHits+Nano100IncMisses)])
+xlim([.5 1.5])
+set(gca,'xtick',[]);
+set(gca,'ytick',[0 Nano100IncHits+Nano100IncMisses],'TickLength',[0.05 0.05])
+set(gca,'ytick',[0 Nano100IncHits+Nano100IncMisses])
+set(gca,'xticklabel',[])
+set(gca,'tickdir','out','FontName','Helvetica')
+
+myxlim=xlim;
+myylim=ylim;
+text(myxlim(1)+diff(myxlim)*0.3,myylim(2)-diff(myylim)*0.05,[num2str(round(Nano100HitRate*100)), '%'],'FontName','Helvetica');
+box off
+
+
+axes(PanelHandles(3,1))
+TrialNum=length(CatchIncTimes);
+TrialLim=TrialNum+0.5;
+for i= 1:length(CatchIncTimes)
+   t= CatchIncTimes(i,1);
+   rt= CatchIncTimes(i,2);
+   s= SpikeTimesNew(SpikeTimesNew>t-1 & SpikeTimesNew<t+2);
+   for j= 1:length(s)
+       l= line([s(j)-t s(j)-t],[i-slen i+slen]);
+       set(l,'color','k')
+   end
+   if rt>0 & rt<2
+       l= line([rt rt],[i-slen i+slen]);
+       set(l,'color','r')
+       set(l,'linewidth',HitLineWidth)
+   end
+end
+axis([-1 2 0 TrialLim])
+set(gca,'xticklabel',[])
+set(gca,'ytick',[0 TrialNum])
+set(gca,'xtick',[-1 0 1])
+set(gca,'tickdir','out','FontName','Helvetica')
+myxlim=xlim;
+myylim=ylim;
+
+
+axes(PanelHandles(3,2))
+p= bar([CatchIncHits]);
+set(p,'facecolor','none')
+ylim([0 CatchIncHits+CatchIncMisses])
+xlim([.5 1.5])
+set(gca,'xtick',[]);
+set(gca,'ytick',[0 CatchIncHits+CatchIncMisses],'TickLength',[0.05 0.05])
+set(gca,'xticklabel',[])
+set(gca,'tickdir','out','FontName','Helvetica')
+
+myxlim=xlim;
+myylim=ylim;
+text(myxlim(1)+diff(myxlim)*0.3,myylim(2)-diff(myylim)*0.05,[num2str(round(CatchHitRate*100)), '%'],'FontName','Helvetica');
+box off
+
+axes(PanelHandles(4,1))
+TrialNum=length(Nano25ConIncTimes);
+TrialLim=TrialNum+0.5;
+
+for i2= 1:length(microstim_sub_trials)
+   i= microstim_sub_trials(i2);
+   t= MicroIncTimes(i,1);
+   rt= MicroIncTimes(i,2);
+   s= SpikeTimesNew(SpikeTimesNew>t-1 & SpikeTimesNew<t+2);
+   for j= 1:length(s)
+       if s(j)<t | s(j)>t+0.2
+           l= line([s(j)-t s(j)-t],[i2-slen i2+slen]);
+           set(l,'color','k')
+       end 
+   end
+   if rt>0 & rt<2
+       l= line([rt rt],[i2-slen i2+slen]);
+       set(l,'color','r')
+       set(l,'linewidth',HitLineWidth)
+   end
+end
+axis([-1 2 0 TrialLim])
+set(gca,'ytick',[0 TrialNum])
+set(gca,'xtick',[-1 0 0.1 1 1.2])
+set(gca,'xticklabel',{-1;0;'';1;''})
+set(gca,'tickdir','out','FontName','Helvetica')
+myxlim=xlim;
+myylim=ylim;
+
+
+axes(PanelHandles(4,2))
+p= bar([MicroIncHits]);
+set(p,'facecolor','none')
+ylim([0 MicroIncHits+MicroIncMisses])
+xlim([.5 1.5])
+set(gca,'xtick',[]);
+set(gca,'ytick',[0 MicroIncHits+MicroIncMisses],'TickLength',[0.05 0.05])
+set(gca,'xticklabel',[])
+set(gca,'tickdir','out','FontName','Helvetica')
+
+myxlim=xlim;
+myylim=ylim;
+text(myxlim(1)+diff(myxlim)*0.3,myylim(2)-diff(myylim)*0.05,[num2str(round(MicroHitRate*100)), '%'],'FontName','Helvetica');
+box off
+
+end
+% 
+% 

codes/Nanostimulation analysis/Fig4Science_data_population_plots.m

@@ -0,0 +1,125 @@
+function Fig4Science_data_population_plots
+close all
+MyFigure = figure('Name',mfilename, 'PaperUnits', 'centimeters');
+PlotsDir='U:\Larkum\PRh paper\Matlab\MicrostimHPC\Figures\Matthew Figures\Fig 4';
+
+NumRows=3;
+NumCols=1;
+MarkSize=8;
+Width = 10*NumCols;
+Height = 6*NumRows;
+PixelsPerCm = 85/2.54;
+set(gcf, 'PaperSize', [Width Height])
+set(gcf, 'PaperPosition', [.5 1 Width Height])
+set(gcf, 'Position', get(gcf, 'PaperPosition') * PixelsPerCm);
+set(gcf,'DefaultAxesFontSize',10, 'DefaultTextFontName', 'Helvetica','DefaultAxesFontName', 'Helvetica')
+
+PanelLetters = cell(NumRows,NumCols);
+
+Dimension = size(PanelLetters);
+RightMargin = .3;
+LeftMargin = .2;
+TopMargin = .05;
+BottomMargin = .1;
+InterPanelWidth = eps;
+InterPanelHeight = .05;
+LetterFontSize = 12;
+LetterX = -.05;
+LetterY = -.035;
+PanelHandles = MultiPanel(MyFigure,Dimension, [LeftMargin, RightMargin, BottomMargin, TopMargin], [InterPanelWidth, InterPanelHeight], [LetterX, LetterY], PanelLetters, LetterFontSize);
+
+load('Fig4Science_data_population_plots_vars.mat');
+PlotsDir='U:\Larkum\PRh paper\Matlab\MicrostimHPC\Figures\Matthew Figures\Fig 4';
+
+MicroCol='Brown';
+CatchCol='Black';
+LFCol='IndianRed';
+BurstCol='DarkRed';
+
+SelectedCellInd=20;
+
+axes(PanelHandles(1,1))
+ 
+pCond1= plot(CatchHitRate*100,Nano25HitRate*100,'o');
+set(pCond1,'markersize',MarkSize,'MarkerEdgeColor','k','MarkerFaceColor','w');
+hold all;
+pCell= plot(CatchHitRate(SelectedCellInd)*100,Nano25HitRate(SelectedCellInd)*100,'o');
+set(pCell,'markersize',MarkSize,'MarkerEdgeColor','r','MarkerFaceColor','r');
+hold all;
+
+l= line([0 100], [0 100],'Color','k','LineStyle',':');
+axis([-2.5 100 -2.5 100]);
+set(gca,'xtick',0:20:100);
+set(gca,'ytick',0:20:100);
+ylabel('Stimulation hits (%)','FontName','Helvetica');
+
+set(gca,'TickDir','out','FontName','Helvetica');
+
+[H,PValue]=ttest(CatchHitRate,Nano25HitRate,0.05,'left');
+
+myxlim=xlim;
+myylim=ylim;
+
+text(myxlim(1)+diff(myxlim)*0.01,95,cat(2,'P= ',num2str(PValue,'%.2f')),'FontName','Helvetica');
+ box off;
+
+
+
+axes(PanelHandles(2,1))
+
+pCond1= plot(CatchHitRate*100,Nano100HitRate*100,'o');
+set(pCond1,'markersize',MarkSize,'MarkerEdgeColor','k','MarkerFaceColor','w');
+hold all;
+pCell= plot(CatchHitRate(SelectedCellInd)*100,Nano100HitRate(SelectedCellInd)*100,'o');
+set(pCell,'markersize',MarkSize,'MarkerEdgeColor','r','MarkerFaceColor','r');
+hold all;
+
+
+l= line([0 100], [0 100],'Color','k','LineStyle',':');
+axis([-2.5 100 -2.5 100]);
+set(gca,'xtick',0:20:100);
+set(gca,'ytick',0:20:100);
+ylabel('Stimulation hits (%)','FontName','Helvetica');
+
+set(gca,'TickDir','out','FontName','Helvetica');
+
+[H,PValue]=ttest(CatchHitRate,Nano100HitRate,0.05,'left');
+
+
+myxlim=xlim;
+myylim=ylim;
+
+text(myxlim(1)+diff(myxlim)*0.01,95,cat(2,'P= ',num2str(PValue,'%.1f')),'FontName','Helvetica');
+box off;
+
+ axes(PanelHandles(3,1))
+pCond1= plot(CatchHitRate*100,MicroHitRate*100,'o');
+set(pCond1,'markersize',MarkSize,'MarkerEdgeColor','k','MarkerFaceColor','w');
+hold all;
+pCell= plot(CatchHitRate(SelectedCellInd)*100,MicroHitRate(SelectedCellInd)*100,'o');
+set(pCell,'markersize',MarkSize,'MarkerEdgeColor','r','MarkerFaceColor','r');
+hold all;
+
+
+l= line([0 100], [0 100],'Color','k','LineStyle',':');
+axis([-2.5 100 -2.5 100]);
+set(gca,'xtick',0:20:100);
+set(gca,'ytick',0:20:100);
+xlabel('False positives (%)','FontName','Helvetica');
+ylabel('Stimulation hits (%)','FontName','Helvetica');
+
+set(gca,'TickDir','out','FontName','Helvetica');
+
+[H,PValue]=ttest(CatchHitRate,MicroHitRate,0.05,'left');
+
+
+myxlim=xlim;
+myylim=ylim;
+
+text(myxlim(1)+diff(myxlim)*0.01,5,cat(2,'P= ',num2str(PValue,'%e')),'FontName','Helvetica');
+box off;
+ 
+ 
+cd(PlotsDir);
+saveas(MyFigure,mfilename,'epsc');
+end

codes/Two-photon imaging analysis/PRh axons/PSTCaTrace.m

@@ -0,0 +1,13 @@
+function snippet = PSTCaTrace(stimframe, trace, preStim, postStim)
+
+%fs = 58.62;
+stimNum = length(stimframe);
+snippet = zeros(postStim+preStim+1, stimNum);
+
+n = length(trace); % total number of frames
+
+for i = 1:stimNum
+    if stimframe(i) + postStim <= n
+        snippet(:,i) = trace(stimframe(i)-preStim:stimframe(i)+postStim);
+    end
+end

codes/Two-photon imaging analysis/PRh axons/TwoPhotonAnalysis.m

@@ -0,0 +1,110 @@
+load('PRh axons.mat', 'Mouse16D')
+%load('200116D_200129_trace.mat')
+%%
+fs = 30; %30.3017
+preStim = 59; %frame %2 sec
+postStim = 90; %frame %3 sec
+
+% Before learning
+pre_stimframe = [Mouse16D.pre10uA_0127.stimframe];
+%pre_stimframe = CleanStimframe(pre_stimframe);
+pre_trace = session1(1).normtrace;
+
+prePSTCa = PSTCaTrace(pre_stimframe, pre_trace, preStim, postStim);
+
+% After learning
+post_stimframe = [Mouse16D.pre10uA_0128.stimframe];
+%post_stimframe = CleanStimframe(post_stimframe);
+post_trace = session2(1).normtrace;
+
+postPSTCa = PSTCaTrace(post_stimframe, post_trace, preStim, postStim);
+
+
+% Reward only
+reward_stimframe = [Mouse16D.reward_0128_2.stimframe];
+%reward_stimframe = CleanStimframe(reward_stimframe);
+
+reward_trace = session2(4).normtrace;
+
+rewardPSTCa = PSTCaTrace(reward_stimframe, reward_trace, preStim, postStim);
+
+
+%%
+
+edges = -2:1/fs:3;
+figure
+boundedline(edges(1:end-1), mean(prePSTCa,2), std(prePSTCa,[],2)/sqrt(size(prePSTCa,2)), 'k-', ...,
+    edges(1:end-1), mean(postPSTCa,2), std(postPSTCa,[],2)/sqrt(size(postPSTCa,2)), 'r-')
+xlim([-1 2])
+legend('prelearning', 'postlearning')
+ylabel('dF/F0 (%)')
+xlabel('Time (s)')
+title('Mouse 16D')
+
+
+% Reward only
+figure
+boundedline(edges(1:end-1), mean(rewardPSTCa,2), std(rewardPSTCa,[],2)/sqrt(size(rewardPSTCa,2)), 'g-')
+legend('Reward only')
+ylabel('dF/F0 (%)')
+xlabel('Time (s)')
+title('Mouse 16D')
+xlim([-1 2])
+
+
+
+%%
+
+Mouse16D.pre_trace = pre_trace;
+Mouse16D.pre_stimframe = pre_stimframe;
+Mouse16D.prePSTCa = prePSTCa;
+Mouse16D.post_trace = post_trace;
+Mouse16D.post_stimframe = post_stimframe;
+Mouse16D.postPSTCa = postPSTCa;
+Mouse16D.reward_trace = reward_trace;
+Mouse16D.reward_stimframe = reward_stimframe;
+Mouse16D.rewardPSTCa = rewardPSTCa;
+
+
+
+%% Statistics
+% Mean calcium amplitudes (between t = 0 and 1.5 s after stimulus) before and after learning were compared by Mann-Whitney U test.
+
+[p h] = ranksum(mean(Mouse15C.prePSTCa(61:105,:)), mean(Mouse15C.postPSTCa(61:105,:)))
+
+[p h] = ranksum(mean(Mouse16D.prePSTCa(61:105,:)), mean(Mouse16D.postPSTCa(61:105,:)))
+
+
+[p h] = ranksum(mean(Mouse31F.prePSTCa(61:105,:)), mean(Mouse31F.postPSTCa(61:105,:)))
+
+
+
+%% Correction fit
+prePSTCa = Mouse15C.prePSTCa;
+[fx, Corr_prePSTCa] = CorrectionFit(prePSTCa,edges);
+
+postPSTCa = Mouse15C.postPSTCa;
+[fx, Corr_postPSTCa] = CorrectionFit(postPSTCa,edges);
+
+
+figure
+boundedline(edges(1:end-1), mean(Corr_prePSTCa,2), std(Corr_prePSTCa,[],2)/sqrt(size(Corr_prePSTCa,2)), 'k-', ...,
+    edges(1:end-1), mean(Corr_postPSTCa,2), std(Corr_postPSTCa,[],2)/sqrt(size(Corr_postPSTCa,2)), 'r-')
+xlim([-1 2])
+
+
+
+[p h] = ranksum(mean(Corr_prePSTCa(61:105,:)), mean(corr_postPSTCa(61:105,:)))
+
+Mouse15C.Corr_prePSTCa = Corr_prePSTCa;
+Mouse15C.Corr_postPSTCa = Corr_postPSTCa;
+
+%% Pooled
+prePSTCa = [Mouse15C.prePSTCa Mouse16D.prePSTCa Mouse31F.prePSTCa];
+postPSTCa = [Mouse15C.postPSTCa Mouse16D.postPSTCa Mouse31F.postPSTCa];
+[p h] = ranksum(mean(prePSTCa(61:105,:)), mean(postPSTCa(61:105,:)))
+
+rewardPSTCa = [Mouse15C.rewardPSTCa Mouse16D.rewardPSTCa Mouse31F.rewardPSTCa];
+[p h] = signrank(mean(rewardPSTCa(31:60,:)), mean(rewardPSTCa(61:105,:)))
+
+