Delete 'FigsDemo.m'

FigsDemo.m

@@ -1,327 +0,0 @@
-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
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-