%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%           M. Oh, S. Zhao, V. Matveev and F. Nadim (2012)
%    Neuromodulatory changes in short-term synaptic dynamics may be 
%   mediated by two distinct mechanisms of presynaptic calcium entry
%            Journal of Computational Neuroscience. 
%
%  Plot of Figure 3 using data generated by sweepGsynRandModel1.m
%                     Victor Matveev, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
global prefix plotFlag Data tArray dFit;
[ Data, dFit, tArray ] = getData(550);

Filter = 8.5;  
filter = 0.3;
Model = 1;
setNum = 2;

ms = 7;
figure('Position', [0 0 1000 600] );
ver1 = 900; ver2 = 500;

hor1 = 0.1; hor2 = 0.55; 
w = 0.42;   h = 0.065; dh = 0.08;
ver1 = 0.9; ver2 = 0.62;

hor3 = 0.1;  ver3 = 0.07;
hor4 = 0.55; ver4 = 0.07;
W = 0.3; H = 0.27;

panels = [ [hor1, ver1, w, h]; [hor1, ver1-dh, w, h]; [hor1, ver1-2*dh, w, h];  ...
           [hor1, ver2, w, h]; [hor1, ver2-dh, w, h]; [hor1, ver2-2*dh, w, h];  ...
           [hor2, ver1, w, h]; [hor2, ver1-dh, w, h]; [hor2, ver1-2*dh, w, h];  ...
           [hor2, ver2, w, h]; [hor2, ver2-dh, w, h]; [hor2, ver2-2*dh, w, h];  ...
           [hor3, ver3, W, H]; [hor4, ver4, W, H] ];


clrs = [ [0 0 1]; [1 0 0]; [0 0.5 0]; [0 0 0]; [1 0 1]; [0.5 0.5 0.5]; [1 1 0]; [0 1 1]; [0.8 0.8 0.8] ];
clri = [ [0.6 0.6 0.6]; [0 0 0]; [1 0.6 0.6]; [1 0 0] ];
Titles = {'{\bf(A)}      Control, \DeltaV = 20 mV',   '{\bf(B)}      Control, \DeltaV=40mV',...
          '{\bf(C)}      Proctolin, \DeltaV = 20 mV', '{\bf(D)}      Proctolin, \DeltaV=40mV'};

minParams    = [];  minDev       = [];
minCost      = [];  normParams   = [];
filterParams = [];  

    % tauH tauCa gsyn minfVproc minfSproc tauMmax tauMV tauMS gm hinfV hinfS 
    %   1    2     3     4         5         6      7     8    9  10    11   
    % K4 tauMcntrl minfVcnt minfScnt  gCa 
    % 12    13       14        15     16  
    
    Filter = 8.5;   % Model 1 was run with new cost function, which had larger trace dev contribution
    DataRandGsynModel1_new;
    setNum = [1];    
    [Y, I] = sort(minCost);

    
for s = setNum %length(Sets)
    
    X = minParams( I(s), : );    minfCntrl = @(V) 1./(1+exp(-(V+X(14))/X(15)));

    minfProct = @(V) 1./(1+exp(-(V+X(4))/X(5)));
    tauM = @(V) X(6) ./ cosh((V+X(7))/X(8));        %tauM = @(V) X(6) + X(7) ./ (1 + exp( (V+X(8))/X(9)));
    Vrange = linspace(-70,-10,80);
    
    subplot('Position', panels(13,:) ); hold on;

    plot(Vrange, minfCntrl(Vrange), 'k--', 'linewidth', 2);
    plot(Vrange, minfProct(Vrange), 'r-', 'linewidth', 2);
    legend('Control', 'Proctolin', 'location', 'best');
    plot([-60], minfCntrl(-60),  'o', 'markerfacecolor', [1 1 1],   'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-40], minfCntrl(-40),  'o', 'markerfacecolor', clri(1,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-20], minfCntrl(-20),  'o', 'markerfacecolor', clri(2,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-60], minfProct(-60),  'o', 'markerfacecolor', [1 1 1],   'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms);
    plot([-40], minfProct(-40),  'o', 'markerfacecolor', clri(3,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-20], minfProct(-20),  'o', 'markerfacecolor', clri(4,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-60 -60], [0 1], 'k:', [-40 -40], [0 1], 'k:', [-20 -20], [0 1], 'k:');
    title('{\bf(E)}      I_{Ca} activation in control vs proctolin', 'fontsize', 11);
    xlabel('V_{LP} (mV)'); ylabel('m_{\infty}(V)');
    axis tight;
    
    subplot('Position', panels(14,:) ); hold on;
    plot(Vrange, tauM(Vrange), 'r-', 'linewidth', 2);
    tauMin = min(tauM(Vrange)); tauMax  = max(tauM(Vrange));
    plot([-60 -60], [tauMin tauM(-60)], 'k:', [-40 -40], [tauMin tauM(-40)], 'k:', [-20 -20], [tauMin tauM(-20)], 'k:');
    plot([-60], tauM(-60),  'o', 'markerfacecolor', [1 1 1],   'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms);
    plot([-40], tauM(-40),  'o', 'markerfacecolor', clri(3,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-20], tauM(-20),  'o', 'markerfacecolor', clri(4,:), 'markeredgecolor', [0 0 0], 'linewidth', 1, 'markersize', ms );
    plot([-70 -10], X(13)*[1 1], 'k--', 'linewidth', 2);
    set(gca, 'Yscale', 'log');
    title('{\bf(F)}      I_{Ca} activation time constant in proctolin', 'fontsize', 11);
    xlabel('V_{LP} (mV)'); ylabel('\tau_m(V) (ms)');
    axis tight;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

X = abs(X);

VCa = 100;    Vsyn = -80;  VpostRest = -55;
C0 = 0;       K4 = X(12);  % Main difference compared to previous   
lambda = 0.1; 

tauH  = X(1);  tauCa = X(2);
gCa   = X(16); gm    = X(9);
gamma = X(3);
tauCa = X(2);

hinf = @(V) 1 ./ (1 + exp((V+X(10))/X(11)));
gsyn = @(C4) K4 * C4 ./ (C4 + K4);

options = odeset('RelTol', 1e-6);     
ii = 1;

for cond = [0 1]   %%%%% iterate over condition = [ control proctolin ] %%%

    if cond == 0  % control
        minf = @(V) 1./(1 + exp(-(V+X(14))/X(15)) );
        tauM = @(V) X(13);
     else
        minf = @(V) 1./(1 + exp(-(V+X(4))/X(5)) );
        tauM = @(V) X(6) ./ cosh((V+X(7))/X(8));
    end

    Vpre = -60;
    % Variables: m, h, Vpost, Ca; 
    f = @(t,Y) [ (minf(Vpre) - Y(1)) / tauM(Vpre); (hinf(Vpre) - Y(2)) / tauH; ...
                 -gamma*gsyn(Y(4).^4).*(Y(3)-Vsyn)-gm*(Y(3)-VpostRest); -(Y(4)-C0)/tauCa-lambda*gCa*Y(1).^2.*Y(2).*(Vpre-VCa); ];
    IC = [minf(Vpre), hinf(Vpre), VpostRest, C0];
    IC0 = fsolve( @(x) f(0,x), IC, optimset('TolFun',1e-8,'MaxFunEvals',5000,'Display','off') );
    
    
    for DV = [20 40] %%%%%%%% iterate over pulse amplitude %%%%%%%%%%%%%%%%
        
        T = []; Y = []; ICA = []; VLP = []; MM=[]; HH=[];
        T0 = 0; IC = IC0;
        
        for i = 0 : 10;  %%%%%%%% iterate over pulse number %%%%%%%%%%%%
           DT = 100 + 400 * (i > 0);
           Vpre = -60 + DV * mod(i,2);
           
           % Variables: m, h, Vpost, Ca; 
           f = @(t,Y) [ (minf(Vpre) - Y(1)) / tauM(Vpre); (hinf(Vpre) - Y(2)) / tauH; ...
                     -gamma*gsyn(Y(4).^4).*(Y(3)-Vsyn)-gm*(Y(3)-VpostRest); -(Y(4)-C0)/tauCa-lambda*gCa*Y(1).^2.*Y(2).*(Vpre-VCa); ];        
                 
           T1 = T0 + DT;
           [Tn, Yn] = ode15s(f, [T0 T1], IC, options);
           T = [T; Tn]; Y = [Y; Yn];  VLP = [VLP; Vpre*ones(size(Tn)) ];
           ICA = [ICA; lambda*gCa*Yn(:,1).^2.*Yn(:,2).*(Vpre-VCa) ];
           MM = [MM; Yn(:,1)]; HH = [HH; Yn(:,2) ];
           T0 = T1;
           IC = Yn(end,:);

        end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over pulses  
               
            dFit(ii, :) = squeeze( extract(T, Y(:,3)-VpostRest, tArray) );
 
            
            subplot('Position', panels(1+(ii-1)*3, :));
            set(gca,'box','off', 'XTick', []);
            clr = [1 0 1] + [0 0.5 0]*(ii==3);
            plot(T, ICA, 'color', clri(ii,:), 'linewidth', 2 ); hold on;
            if ii > 2
                mm = min(ICA);
                axis([0 5000 mm 0]);
                set(gca,'Ytick',[], 'XTick', [], 'box','off'); 
                subplot('Position', panels(1+(ii-3)*3, :));
                axis([0 5000 mm 0]);
            else
                ylabel('I_{Ca}');
            end
            set(gca,'Xtick',[], 'box','off'); 
            
            subplot('Position', panels(2+(ii-1)*3,:));
            plot(T, Y(:,4), 'color', clri(ii,:), 'linewidth', 2 ); 
            set(gca,'Xtick',[], 'box','off');
            if ii > 2
                mm = max(Y(:,4));
                axis([0 5000 0 mm]);
                set(gca,'Ytick',[]);
                subplot('Position', panels(2+(ii-3)*3,:));
                axis([0 5000 0 mm]);
            else
                ylabel('Ca');
            end
            ii = ii + 1;
                
    end %%%%%%%%%%%%%%  loop over DV = [20 40]   %%%%%%%%%%%%%%%%%%%%%%%%%

end %%%%%%%%% loop over condition = [ control proctolin ] %%%%%%%%%%%%%%%%

clr = [0 0.6 0]; lw = 1.5;
scale = sum(sum(Data)) / sum(sum(dFit));

for ii = 1 : 4
    subplot('Position', panels(1+(ii-1)*3, :));
    title(Titles{ii}, 'fontsize', 11);
    subplot('Position', panels(3+(ii-1)*3, :)); hold on;
    plot(tArray, -55+Data(ii,:)/scale, '-', 'color', clr, 'linewidth', lw);
    plot(tArray, -55+dFit(ii,:), 'color', clri(ii,:), 'linewidth', 2);
    set(gca,'box','off');
    if ii == 1 || ii == 3
        set(gca, 'XTick', [], 'box', 'off');
    else
        xlabel('Time (ms)');
    end
    if ii > 2
                mm = min(min([dFit(ii,:), Data(ii,:)/scale]));
                axis([0, 5000, mm-55, -55]);
                set(gca,'Ytick',[], 'box','off');
                subplot('Position', panels(3+(ii-3)*3, :)); 
                axis([0, 5000, mm-55, -55]);
    else
        ylabel('V_{PD}');
    end
    if ii == 1 || ii == 3
        legend('Data', 'Model', 'Location', 'Best');
    end
end

