%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%        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 Fig. 4A, using data produced by equiM5_vs_Dist
%                  Victor Matveev, January 8, 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Imax   = 0.05;          % Single-channel current in pA
Btotal = 200;           % Total buffer concentration
R0     = 1.112090e-09;  % Background release rate at [Ca]bgr = 50 nM
M      = 5;             % Number of channels in a ring
ppA    = [0.6];         % Open channel fraction

special = [29.9863   48.7956   65.3833   79.5306   92.2105];
cs = [0 0 0; 0.6 0.6 1; 0 1 0; 1 0 0; 1 1 0 ];

titleFontSize = 15;
labelFontSize = 12;
axisFontSize  = 12;

name = ['equi_M5_Imax', num2str(Imax), 'pA_B', num2str(Btotal), 'uM_70x50.dat'];

f = fopen(name,'r');
B = fscanf(f, '%f');
fclose(f);

if B(1) == B(2)
    B = B(2:end);
end;

K     = 22;
iDist = length(B)/K;
A     = reshape(B,K,iDist);
Dist  = squeeze( A(1,:,:) );

for ip = 1:length(ppA)

    pp = ppA(ip);
    mica = [];
    mch = [];
    rarray = [];
    ls = 'ABCDEF';

    for iD = 1 : iDist

        R = [num2str(R0),'*(1-p)^',num2str(M)];
        mchNumer = '';
        iProx = 0;

        base = K * (iD - 1);

        nOpenArray   = [1 2 2 3 3 4 5];  % Number of open channels in each of the 7 configurations
        multiplicity = [5 5 5 5 5 5 1];  % Multiplicity of each of the 7 configurations

        for ind = 1 : 7 
                        prob   = B( 1 + ind + base);
                        prob   = prob * (prob > 0);

                        nOpen  = nOpenArray(ind);
                        nClose = M - nOpen;
                        Factor = multiplicity(ind);

                        probStr = strcat(num2str(Factor),'*',num2str(prob),'*p^',num2str(nOpen),'*(1-p)^',num2str(nClose));
                        R       = strcat(R,  ' + ',probStr);
                        if nOpen > 0
                            mchNumer = strcat(mchNumer,' + ',num2str(nOpen),'*',probStr);
                        end;          
        end

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

    end;

    figure(1);
    subplot(1,length(ppA),ip);
    plot(Dist*1000,mica,'-k','linewidth',2); 
    hold on;
    xlabel('Distance (nm)', 'fontsize', labelFontSize);
    ylabel('m_{ICa}', 'fontsize', labelFontSize);
    if ip ~= 1
        legend('m_{ICa}','location','best');
    end;
    axis([20 Dist(end)*1000 micaMin micaMax]);
    set(gca,'box','off');

    for j = 1 : length(special)
        yy = spline(Dist*1000, mica, special(j)); 
        plot(special(j), yy, 'o', 'markeredgecolor', [0 0 0], 'markerfacecolor', cs(j,:), 'markersize', 10 );
    end;

    figure(3);
    subplot(1,length(ppA),ip); 
    hold on;
    plot(rarray, mica, 'om', 'linewidth',2); 
    xlabel('Release', 'fontsize', labelFontSize);
    ylabel('m_{ICa}', 'fontsize', labelFontSize);

end

distR    = rarray;
distD    = Dist*1000;
distMica = mica;
distMch  = mch;

