%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%        Victor V. Matveev, Richard Bertram and Arthur Sherman (2011)
% Ca2+ Cooperativity of Exocytosis as a Measure of Ca2+ Channel Domain Overlap.
%                      Brain Research 1398: 126-138.  
%
%  This script plots Figure 4E using data produced by CalC script Fig. 4E
%                  Victor Matveev, January 8, 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
K = 6;  % Data per line: [ config multiplicity nOpen P(R) [Ca] T_max ]

Marray    = 2:12;
po        = 0.6;
distArray = [60 100];

Btotal = 200;
Imax = 0.1;

Mica5 = @(p,M) (1-30*p+30*p.*M+75*p.^2.*M.^2-225*p.^2.*M+150*p.^2+40*p.^3.*M.^3-240*p.^3.*M.^2+440*p.^3.*M-240*p.^3+5*p.^4*M.^4-50*p.^4.*M.^3+175*p.^4.*M.^2-250*p.^4.*M+120*p.^4)...
              ./(1-15*p+15*p.*M+25*p.^2.*M.^2-75*p.^2.*M+50*p.^2+10*p.^3.*M.^3-60*p.^3.*M.^2+110*p.^3.*M-60*p.^3+p.^4.*M.^4-10*p.^4.*M.^3+35*p.^4.*M.^2-50*p.^4.*M+24*p.^4);
Mica6 = @(p,M) (1-62*p+62*p.*M+270*p.^2.*M.^2-810*p.^2.*M+540*p.^2+260*p.^3.*M.^3-1560*p.^3.*M.^2+2860*p.^3.*M-1560*p.^3+75*p.^4.*M.^4-750*p.^4.*M.^3+2625*p.^4.*M.^2-3750*p.^4.*M+1800*p.^4-90*p.^5.*M.^4+510*p.^5.*M.^3-1350*p.^5.*M.^2+1644*p.^5.*M-720*p.^5+6*p.^5.*M.^5)...
              ./(1-31*p+31*p.*M+90*p.^2.*M.^2-270*p.^2.*M+180*p.^2+65*p.^3.*M.^3-390*p.^3.*M.^2+715*p.^3.*M-390*p.^3+15*p.^4.*M.^4-150*p.^4.*M.^3+525*p.^4.*M.^2-750*p.^4.*M+360*p.^4-15*p.^5.*M.^4+85*p.^5.*M.^3-225*p.^5.*M.^2+274*p.^5.*M-120*p.^5+p.^5.*M.^5);
Mica4 = @(p,M) (1-14*p+14*p.*M+18*p.^2.*M.^2-54*p.^2.*M+36*p.^2+4*p.^3.*M.^3-24*p.^3.*M.^2+44*p.^3.*M-24*p.^3)./(1+7*p.*M-7*p+6*p.^2.*M.^2-18*p.^2.*M+12*p.^2+p.^3.*M.^3-6*p.^3.*M.^2+11*p.^3.*M-6*p.^3);
MicaMax = @(n,M) M.*(1 - (1-1./M).^n);

cs = [0.5 0.5 0.5; 0 0 0];
ls = {'*-'; 's-'};

figure(2);
titleFontSize = 12;
labelFontSize = 12;
axisFontSize  = 12;
iplot         = 1;

subplot(1, length(poArray), 1);
lstr = { };
hold on;

for iD = 1 : length(distArray);
    
    Dist = distArray(iD);   
    mica = [];
    mch  = [];

    for M = Marray

        name = ['M', num2str(M), '_ica', num2str(Imax), 'pA_B', num2str(Btotal), 'uM_', num2str(Dist), 'nm_70x50.dat'];
        f    = fopen(name,'r');
        B    = fscanf(f, '%f');
        L    = length(B) / K;
        fclose(f);

        R = '';
        R2 = '';
        mchNumer = '';

        for ind = 0 : L-1
            conf   = B(K * ind + 1);
            mult   = B(K * ind + 2);
            nOpen  = B(K * ind + 3);
            prob   = B(K * ind + 4);
            prob   = prob * (prob > 0);
            nClose = M - nOpen;

            probStr  = strcat(num2str(mult),'*',num2str(prob),  '*p^',num2str(nOpen),'*(1-p)^',num2str(nClose));
            probStr2 = strcat(num2str(mult),'*',num2str(prob^2),'*p^',num2str(nOpen),'*(1-p)^',num2str(nClose));
            R        = strcat(R,  ' + ',probStr);
            R2       = strcat(R2, ' + ',probStr2);
            mchNumer = strcat(mchNumer,' + ',probStr,'*',num2str(nOpen));
        end % configurations

        dR = diff(R);
        mica = [mica, p0 ./ subs(R,p0) .* subs(dR,p0) ];
        mch = [mch, subs(mchNumer,p0) ./ subs(R,p0) ];

    end % M loop

    marray = Marray(1):1:Marray(end);
    pmica = spline(Marray,mica,marray);
    plot(marray,pmica,ls{iD},'linewidth',2,'color',cs(iD,:));  hold on;
    lstr = [lstr, ['Dist=',num2str(Dist),'nm'] ];
    drawnow;

end % Distance loop

xlabel('M (# of channels)', 'fontsize', labelFontSize);
ylabel('m_{ICa}', 'fontsize', labelFontSize);

plot([Marray(1) Marray(end)], [ Marray(1) Marray(end)], 'color', [1 0.6 0.6], 'linewidth', 3 ); 
lstr = [lstr, 'diagonal' ];
legend(lstr);
axis([Marray(1) Marray(end) 2 7]);
drawnow;    
