Upload files to ''

FigsDemo.m

@@ -0,0 +1,327 @@
+function FigsDemo
+% FIGSDEMO - demonstration of plotting contents of a processed data file to
+% make many of the figures in the manuscript. Runs in Matlab.
+% This assumes that dD, dEN files containing data from experiment SW421 are
+% in the Matlab workspace.
+
+gray = [0.5 0.5 0.5];
+
+% Get the data from SW421
+load('SW421_dD_dEN_C.mat')
+
+% FIGURE 1A. Plot the raw waveform (blue), fit waveform (red), and
+% deconvolved event signal (green) for the examples in this figure.The
+% EPSCs are located at the following times, s.:
+EPSClocs = [337.0708    1.0000;
+  		337.1683    1.0000;
+  		337.2380    7.0000;
+  		337.2835    5.0000;
+  		337.3419    4.0000;
+  		337.3556    1.0000;
+  		337.4153    3.0000;
+  		337.4400    2.0000;
+  		337.4828    1.0000;
+  		337.6855    6.0000;
+  		337.7531    3.0000;
+  		337.8212    1.0000;
+ 		337.9187    1.0000];
+figure(1); clf
+% Signals are sampled at dD.ut_s s. Convert files to indices of the sigs.
+indices = round(EPSClocs(:,1)/dD.ut_s); fivems = round(0.005/dD.ut_s);
+rw = []; fw = []; ds = []; % To hold the 3 plotted signals
+for jn=1:size(EPSClocs,1)
+    rw = [rw; dD.DataDC(indices(jn)-fivems:indices(jn)+fivems)];
+    fw = [fw; dEN.dataOIL(indices(jn)-fivems:indices(jn)+fivems)];
+    ds = [ds; dEN.xsigL(indices(jn)-fivems:indices(jn)+fivems)];
+end
+tim = 0:dD.ut_s:(length(rw)-1)*dD.ut_s;
+plot(tim, fw, 'r-'); hold on
+plot(tim, rw, 'b-'); plot(tim, ds, 'g-')
+xlabel('Seconds'); ylabel('Currents, pA')
+legend('Fitted waveshape', 'Data', 'Deconv. signal')
+title('Figure 1A')
+
+
+% FIGURE 1-supplementary can be obtained as 
+%     plot(dD.timK, dD.AvgKernel)
+
+
+% FIGURE 2A - histograms of dEN.cluspks, which are the amplitudes of the
+% EPSCS, which contain dEN.nsigpks events.
+thpks = -dEN.cluspks;   % EPSC amplitudes, sign changed to +
+% Bin edges estimated with Scott's rule:
+bw = 3.5*std(thpks)/numel(thpks)^1/3;
+bw = 10*ceil(bw/10); 
+edges = (-50:bw:max(thpks)+bw)';
+fprintf('bw=%g pA.\n', bw)
+% Plot the amplitude histograms
+figure(2); clf
+%i1 = find(dEN.nsigpks==1);   % Indices of monos in the data
+histogram(thpks(dEN.nsigpks==1), edges)
+hold on
+for jn=2:4
+    if jn<4
+        iI = find(dEN.nsigpks==jn);
+    else
+        iI = find(dEN.nsigpks>=4);
+    end
+    histogram(thpks(iI), edges, 'DisplayStyle','stairs', 'LineWidth', 2)       % multis
+end
+histogram(thpks, edges,'DisplayStyle','stairs', 'EdgeColor','k', ...
+                     'LineWidth',2)
+legend('Monophasic','2-events', '3-events', '>=4 events','all')
+xlabel('Peak current, pA'); ylabel('Number')
+title('Figure 2A')
+
+
+% FIGURES 2C, 2D. The data are in Figure234data.mat.
+% Expt H is in row 6; expt T is in row 2.
+load Figure234data; aa = Figure234data;
+figure(25); clf
+nex = size(aa.Nevents,1);
+aaS = sum(aa.Nevents,2);
+aaPL = (ones(1,size(aa.Nevents,2))./aaS).*aa.Nevents; % Compute fractions
+semilogy(aa.colmns, mean(aaPL), 'k-', 'Linewidth', 2)  % Mean of fractions
+hold on
+for jr=1:nex; semilogy(aa.colmns, aaPL(jr,:), 'b-'); end
+xlabel('Number of events'); ylabel('Fraction of EPSCs')
+legend('Mean', 'Indiv. expts.')
+title('Figure 2C')
+
+figure(26); clf
+MA = mean(aa.Mampl)';  % Mean amplitudes across expts
+plot(aa.colmns', MA, 'k-', 'Linewidth', 2); hold on 
+mdl1 = fitnlm(aa.colmns', MA, @(b,xv)b(1)+b(2)*xv+b(3)*xv.^2, ...
+            [MA(1) (MA(end)-MA(1))/length(MA) 1]); % Fit with a quadratic
+plot(aa.colmns', mdl1.Fitted, 'g--', 'LineWidth', 2) % Plot the fit
+for jr=1:nex; plot(aa.colmns, aa.Mampl(jr,:), 'b-'); end % Plot the data
+R2 =  1-mdl1.SSE/mdl1.SST; % Coefficient of determination
+fprintf('Fig 2D Model: Ampl = %g + %g*NEv + %g*NEv^2. R^2=%g, P=%g\n', ...
+      mdl1.Coefficients.Estimate, R2, 1-tcdf(sqrt(R2*(nex-2)/(1-R2)),nex-2))
+xlabel('Number of events'); ylabel('Mean amplitude, pA')
+legend('Mean', 'Quad. fit', 'Indiv. expts')
+title('Figure 2D')
+
+
+% FIGURE 3A - histogram of dEN.clusArea, the areas of EPSCs.
+thare = -dEN.clusArea*dEN.ut_s;   % EPSC amplitudes, sign changed to +,
+                                 % converted to pC = pA*s
+% Bin edges set for 100 bins:
+upperlimit = 0.1*ceil(prctile(thare, 99)/0.1); bw = upperlimit/100;
+edges = (0:bw:upperlimit)';
+fprintf('bw=%g pA.\n', bw)
+% Plot the amplitude histograms
+figure(3); clf
+for jn=1:4
+    if jn<4
+        iI = find(dEN.nsigpks==jn);
+    else
+        iI = find(dEN.nsigpks>=4);
+    end
+    histogram(thare(iI), edges, 'DisplayStyle','stairs', 'LineWidth', 2)
+    hold on
+end
+legend('Monophasic','2-events', '3-events', '>=4 events')
+xlabel('EPSC area, pC'); ylabel('Number')
+title('Figure 3A')
+
+
+% FIGURE 3D. The data are in Figure234data.mat.
+figure(36); clf
+MAR = mean(aa.Marea);  % Mean area across expts
+mdl2 = fitnlm(aa.colmns, MAR, @(b,xv)b(1)+b(2)*xv, ...
+            [MAR(1) (MAR(end)-MAR(1))/length(MAR)]); % Fit with a line
+plot(aa.colmns, MAR, 'k-', 'LineWidth', 2); hold on % Plot avg. data
+plot(aa.colmns', mdl2.Fitted, 'g--', 'LineWidth', 2) % Plot the fit
+for jr=1:nex; plot(aa.colmns,aa.Marea(jr,:), 'b-'); end  % Plot data
+R2A =  1-mdl2.SSE/mdl2.SST; % Coefficient of determination
+fprintf('Fig 3D Model: Area = %g + %g*Nev, R^2=%g, P=%g\n', ...
+ mdl2.Coefficients.Estimate, R2A, 1-tcdf(sqrt(R2A*(nex-2)/(1-R2A)),nex-2))
+xlabel('Number of events'); ylabel('EPSCarea, pC')
+legend('Mean', 'Indiv. expts.')
+title('Figure 3D')
+axv = axis; axis([1 axv(2) 0 axv(4)])
+
+% Expt H is in row 6; expt T is in row 2.
+
+
+% FIGURE 4A. The data are in Figure234data.cmNsig, not normalized by expt.
+% normcm is the normalizer. There are NaNs in the data, where the analysis
+% yielded too few samples.
+figure(4); clf
+nex = size(aa.cmNsig, 2); nec = size(aa.cmNsig, 1);
+cmNmean = zeros(nec,1);  % Compute mean across expts
+for j1=1:nec
+    cmNmean(j1) = mean(aa.cmNsig(j1,~isnan(aa.cmNsig(j1,:))));
+end
+Mcmn = cmNmean/mean(aa.normcm);  % Normalize the mean
+plot((1:nec)', Mcmn, 'k-', 'LineWidth', 3) % Plot normalized mean across expts
+hold on
+mdl4 = fitnlm((1:nec)', Mcmn, @(b,xv)b(1)+b(2)*xv, ...
+            [Mcmn(1) (Mcmn(end)-Mcmn(1))/nec]); % Fit with a line 
+plot((1:nec)',mdl4.Fitted,'g--','Linewidth', 2)  % Plot fit line
+for ju=1:nex
+    plot((1:nec)', aa.cmNsig(:,ju)/aa.normcm(ju), 'b-')  % Plot data
+end
+R2cm =  1-mdl4.SSE/mdl4.SST; % Coefficient of determination of fit
+fprintf('Fig 4A Model: M(nth) = %g + %g*N(N-1th), R^2=%g, P=%g\n', ...
+  mdl4.Coefficients.Estimate, R2cm, 1-tcdf(sqrt(R2cm*(nex-2)/(1-R2cm)), ...
+  nec-2))
+axv = axis; axis([0 axv(2) 0 axv(4)])
+xlabel('Event count in n-1''th cluster')
+ylabel('Mean count in n''th cluster')
+legend('Mean', 'Linear model', 'Data')
+title('Figure 4A')                
+
+% FIGURE 4B. 
+figure(45); clf
+iplB = find(Figure234data.intNN(:,1)<=0.1);  % Only for intervals <=0.1 s
+[xxB,yyB,~] = slideMean(Figure234data.intNN(iplB,1), ...
+     Figure234data.intNN(iplB,2), 9, 1); % Smooth data in sets of 9
+                                              % (~9*33 original EPSCs)
+plot(xxB, yyB, 'k-', 'LineWidth', 2); hold on
+mnxxB = min(xxB); iplB2 = find(xxB<=0.02);  % Model fit only 1st 0.02 s
+mdlB = fitnlm(xxB(iplB2)-mnxxB, yyB(iplB2), ...
+               @(b,xv)b(1)-b(2)*exp(-xv/b(3)), [1 0.3 0.001]);
+R2B =  1-mdlB.SSE/mdlB.SST; % Coefficient of determination of fit
+fprintf('Fig 4B Model: NN = %g - %g*exp(-Int/%g), R^2=%g, P=%g\n', ...
+  mdlB.Coefficients.Estimate, R2B, ...
+  1-tcdf(sqrt(R2B*(length(xxB(iplB2))-2)/(1-R2B)),nec-2))
+plot(xxB(iplB2), mdlB.Fitted, 'g--', 'LineWidth', 2)
+plot(Figure234data.intNN(iplB,1), Figure234data.intNN(iplB,2), 'o', ...
+     'Color', gray, 'MarkerSize',4, 'MarkerFaceColor', gray)
+axv = axis; axis([0 0.02 axv(3:4)])
+line([0 0.1], [1, 1], 'Color', 'k', 'LineStyle', '--')
+xlabel('Interval to prev cluster, s.')
+ylabel('Normalized number of events/EPSC')
+legend('Mean', 'Exponential model', 'Data')
+title('Figure 4B')
+mdlB      % Print the model
+
+
+% FIGURE 4C. Data are in Figure234data.allArea. These are all the data
+% points plotted in Fig 4C, not separated by experiment. Similar data for
+% normalized amplitude versus interval are in Figure234data.allAmpl
+figure(46); clf
+ipl = find(Figure234data.intArea(:,1)<=0.1);  % Only for intervals <=100 ms
+[xx,yy,~] = slideMean(Figure234data.intArea(ipl,1), ...
+     Figure234data.intArea(ipl,2), 21, 0.5); % Smooth data in sets in sets
+                                              % of 21 (~1000 original EPSCs)
+plot(xx, yy, 'b-', 'LineWidth', 2); hold on
+mnxx = min(xx);
+mdlR = fitnlm(xx-mnxx, yy, @(b,xv)b(1)-b(2)*exp(-xv/b(3)), [1 0.3 0.01]);
+R2R =  1-mdlR.SSE/mdlR.SST; % Coefficient of determination of fit
+fprintf('Fig 4C Model: Area = %g - %g*exp(-Int/%g), R^2=%g, P=%g\n', ...
+  mdlR.Coefficients.Estimate, R2R, ...
+  1-tcdf(sqrt(R2R*(length(xx)-2)/(1-R2R)),nec-2))
+plot(xx, mdlR.Fitted, 'g--', 'LineWidth', 2)
+plot(Figure234data.intArea(ipl,1), Figure234data.intArea(ipl,2), '.', ...
+              'Color', gray)
+line([0 0.1], [1, 1], 'Color', 'k', 'LineStyle', '--')
+xlabel('Interval to prev cluster, s.')
+ylabel('Normalized area of cluster')
+legend('Mean', 'Exponential model', 'Data')
+title('Figure 4C')
+mdlR      % Print the model
+% For the analysis of amplitude, mentioned in the ms, repeat the above
+% using Figure234data.intAmpl
+[xxM,yyM,~] = slideMean(Figure234data.intAmpl(ipl,1), ...
+     Figure234data.intAmpl(ipl,2), 21, 0.5); % Smooth data in sets in sets
+mnxxM = min(xxM);
+mdlM = fitnlm(xxM-mnxxM, yyM, @(b,xv)b(1)-b(2)*exp(-xv/b(3)), [1 0.3 0.01]);
+R2M =  1-mdlM.SSE/mdlM.SST; % Coefficient of determination of fit
+fprintf('Not plotted, but same analysis on EPSC amplitude gives:\n')
+fprintf('Fig 4C Model: Ampl = %g - %g*exp(-Int/%g), R^2=%g, P=%g\n', ...
+  mdlM.Coefficients.Estimate, R2M, ...
+  1-tcdf(sqrt(R2M*(length(xxM)-2)/(1-R2M)),nec-2))
+mdlM    % Print the model
+
+% FIGURE 5A. Plot a randomly selected 10 examples of EPSCs with NEv events.
+% EPSCs are from the example dataset in dD, dEN.
+isdone = 'x';
+while isdone~='q' && isdone~='Q'
+    maxntsh = 7;  % Max number of examples to show
+    axtim = 0.005; % Length of x-axis
+    axleng = round(axtim/dEN.ut_s);   % Axis length in indices for axtim s
+    onems = round(0.001/dEN.ut_s);  % 1 ms in indices
+
+    NEv = input('Show EPSCs with ?? events? ');
+    iE = find(dEN.nsigpks==NEv);    % All the EPSCs with NEv events
+    ntsh = min(maxntsh, length(iE));    % The number to show
+    ir = randperm(length(iE)); iEP = iE(ir(1:ntsh)); % randomly choose ntsh exmpls
+    tim = 1000*(0:dEN.ut_s:(axleng-1)*dEN.ut_s)'; % time axis, in ms
+    figure(5); clf; ah = 0.9/ntsh;
+    wpos = get(5,'Position'); set(5,'Position',[wpos(1) wpos(2)+wpos(4)-840  ...
+                                             wpos(3) 840])
+    for je=1:length(iEP)
+        ithEst = dEN.iClust(iEP(je),2) - onems; % Index, start of EPSC - 1 ms
+        ithEen = dEN.iClust(iEP(je),3) + onems; % and end of the EPSC + 1 ms
+        ipad = length(tim) - (ithEen - ithEst + 1); % Zero padding to fill axtim s
+        if ipad<0; ithEen = ithEen + ipad; ipad = 0; end  % If long long one
+        axes('Position', [0.1 0.96-ah*je 0.8 ah-0.02]);
+        plot(tim, [dEN.dataOIL(ithEst:ithEen); zeros(ipad,1)], 'r-', ...
+            'LineWidth', 1.5); hold on
+        plot(tim, [dD.DataDC(ithEst:ithEen); zeros(ipad,1)], 'b-', ...
+            'LineWidth', 1.5)
+        plot(tim, [dEN.xsigL(ithEst:ithEen); zeros(ipad,1)], 'g--', ...
+            'LineWidth', 1)
+        tm0 = ithEst*dEN.ut_s*1000 - 0.1;
+        plot(dEN.iEvClustL{iEP(je)}(:,1)*dEN.ut_s*1000-tm0, dEN.PeakAreas{iEP(je)}, ...
+                 'go', 'MarkerSize', 4, 'MarkerFaceColor', 'g')
+        if je==1
+            legend('Fit waveform','Data','Deconv. signal','Event ampls', ...
+                'AutoUpdate','off')
+        end
+        line([0 axtim], [0 0], 'Color', 'k', 'LineWidth', 0.8)
+        if je==length(iEP)
+            xlabel('Time, ms'); ylabel('EPSC, pA')
+        end
+    end
+    isdone = input('Quit, Again? ','s'); isdone = isdone(1);
+end
+
+
+return
+
+
+
+function [xxo, yyo, Nadv] = slideMean(xxe,yye,N,Fadv)
+% SLIDEMEAN - for each N (xx,yy) points in order sort(xx), compute 
+% xxo=mean(xx) and yyo=mean(yy), advance by (Fadv*N) points between computes.
+%       [xxo, yyo, Nadv] = slideMean(xx,yy,N,Fadv);
+% Fadv =0.5 is 50% overlap, =1 is independent samples, =0 is a bad idea.
+[xx,iso] = sort(xxe); yy = yye(iso);  % Sort data along x axis
+Nadv = round(max(Fadv*N, N/10));  % Number of points to advance between bins
+xxo = []; yyo = [];
+i1 = 1; i2 = N;
+while i2<=length(xx)
+    xxo = [xxo; mean(xx(i1:i2))];
+    yyo = [yyo; mean(yy(i1:i2))];
+    i1 = i1 + Nadv; i2 = i2 + Nadv;
+end
+return
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+