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

Imax = 0.05;            % Single-channel Ca2+ current (in pA)
Dist = 30;              % Radius of channel ring (in nm)
R0   = 1.112090e-09;    % Equilibrium release rate at [Ca]bgr = 50 nM
M    = 5;               % Number of channels in a ring
ppA  = [0.6];           % Open channel fraction

micaMax = 3.1;
micaMin = 2.0;

special = [200 1000 2000 3000 4000 ];
cs      = [0 0 0; 0.6 0.6 1; 0 1 0; 1 0 0; 1 1 0 ];

titleFontSize = 14;
labelFontSize = 12;
axisFontSize  = 12;

name = ['equi_M5_Imax', num2str(Imax), 'pA_Dist', num2str(Dist), 'nm_70x50.dat'];

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

K       = 22;
iBtotal = length(B)/K;
A       = reshape(B,K,iBtotal);
Btot    = squeeze( A(1,:,:) );

for ip = 1:length(ppA)

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

    for iB = 1 : iBtotal

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

    base = K * (iB - 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   = A( 1 + ind + base);
                    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

    if subs(R,pp) > 10*R0
        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) ];
        btot = [btot, Btot(iB)];
    end

    end;

    figure(2);
    subplot(1,length(ppA),ip);
    plot(btot,spline(btot,mica,btot),'-k','linewidth',2); 
    hold on;
    xlabel('B_{total} (\muM)', 'fontsize', labelFontSize);
    ylabel('m_{ICa}', 'fontsize', labelFontSize);
    if ip ~= 1
       legend('m_{ICa}','m_{CH}','Location','best');
    end
    axis([Btot(1) 5000 micaMin micaMax]);
    set(gca,'box','off');

    for j = 1 : length(special)
        yy = spline(btot, 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, 'sr', 'linewidth',2); 
    xlabel('Release', 'fontsize', labelFontSize);
    ylabel('m_{ICa}', 'fontsize', labelFontSize);

end

text(200,micaMax-0.5,['i_{Ca}=', num2str(Imax), 'pA;  Dist=', num2str(Dist), ' nm']);
PRmap = spline(distR, distD, rarray);

figure(5); 
plot(btot, PRmap, 'k-', 'linewidth', 2); 
hold on;
xlabel('B_{total} (\muM) at fixed distance', 'fontsize', labelFontSize); 
ylabel('Distance (nm), at fixed B_{total}', 'fontsize', labelFontSize);
title('Effective distance as a function of buffer concentration', 'fontsize', titleFontSize);
axis([0 5000 10 100]);

Dspecial = spline(btot,PRmap,special) % enter into M5equi_Dist

for j = 1 : length(special)
    plot(special(j),    Dspecial(j),    'o', 'markeredgecolor', [0 0 1], 'markerfacecolor', cs(j,:), 'markersize', 10 );
    text(special(j)+70, Dspecial(j)-9,  sprintf('m_{ICa}=%.3g', 0.5 * (spline(btot,  mica, special(j)) + spline(distD, distMica, Dspecial(j)) )));
    text(special(j)+70, Dspecial(j)-4,  sprintf('m_{CH}=%.3g', 0.5 * (spline(btot,  mch, special(j)) + spline(distD, distMch, Dspecial(j)) )));
    text(special(j)+70, Dspecial(j)-14, sprintf('R=%.1e', (spline(btot, rarray, special(j))/spline(btot, rarray, special(1)) ) ) );
end



