RatANspontact/FigsDemo.m

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