%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%           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. 
%
%          Objective (fitness) function for Model 2
%                     Victor Matveev, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Cost, Factors, Status] = FitCostGsynModel2(X)           
% Parameters: 
% gMI tauCa gamma kpV kpS gm tauH hinfV hinfS 
%  1   2     3     4   5   6   7   8    9   
% K4  tauM  minfV minfS gCa  km 
% 10   11     12    13  14   15

global prefix plotFlag Data tArray dFit;
X = abs(X);

VCa = 100;    Vsyn = -80;  VpostRest = -55;
C0 = 0;       K4 = X(10);  % Main difference compared to previous   
lambda = 0.1; 
 
gCa   = lambda*X(14);  gamma = X(3);  gm   = X(6);
tauCa = X(2);          tauH  = X(7);  tauM = X(11);

minf = @(V) 1./(1+exp(-(V+X(12))/X(13)));
hinf = @(V) 1 ./ (1 + exp((V+X(8))/X(9)));
gsyn = @(C4) K4 * C4 ./ (C4 + K4);
kp = @(V) 1. / (1 + exp(-(V+X(4))/X(5)) );
km = 1e-4 + X(15);

options = odeset('RelTol', 1e-6);
Vmax = zeros(4,11);
Vend = zeros(4,11);
Dev  = zeros(4,1);
ii = 1;       

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

    gMI = cond * lambda * X(1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
    Vpre = -60;
    minfty = minf(Vpre); hinfty = hinf(Vpre);
    kpV = kp(Vpre);  rate = kpV + km;
    % Variables: m, h, Vpost, Ca, x; 
    f = @(t,Y) [ (minfty - Y(1)) / tauM; (hinfty - Y(2)) / tauH; ...
                 -gamma*gsyn(Y(4).^4).*(Y(3)-Vsyn)-gm*(Y(3)-VpostRest); -(Y(4)-C0)/tauCa-(gCa*Y(1).^2.*Y(2)+gMI*Y(5)).*(Vpre-VCa); ...
                 kpV - rate*Y(5)]; 
    IC = [minf(Vpre), hinf(Vpre), VpostRest, C0, 0];
    %IC0 = fsolve( @(x) f(0,x), IC, optimset('Display','off') );
    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 = [];
        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, x;   
            minfty = minf(Vpre); hinfty = hinf(Vpre);
            kpV = kp(Vpre);  rate = kpV + km;
            % Variables: m, h, Vpost, Ca, x; 
            f = @(t,Y) [ (minfty - Y(1)) / tauM; (hinfty - Y(2)) / tauH; ...
                         -gamma*gsyn(Y(4).^4).*(Y(3)-Vsyn)-gm*(Y(3)-VpostRest); -(Y(4)-C0)/tauCa-(gCa*Y(1).^2.*Y(2)+gMI*Y(5)).*(Vpre-VCa); ...
                         kpV - rate*Y(5)]; 
         
           T1 = T0 + DT;
           [Tn, Yn] = ode15s(f, [T0 T1], IC, options);
           T = [T; Tn]; Y = [Y; Yn];
           T0 = T1;
           IC = Yn(end,:);
 
           Vend(ii, i+1) = VpostRest - Yn(end,3); 
           Vmax(ii, i+1) = max( VpostRest - Yn(:,3) );

        end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Loop over pulses  
        
        dFit(ii, :) = squeeze( extract(T, Y(:,3)-VpostRest, tArray) );
        ii = ii + 1;
       
    end %%%%%%%%%%%%%%  loop over DV = [20 40]   %%%%%%%%%%%%%%%%%%%%%%%%%

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

scale = sum(sum(Data)) / sum(sum(dFit)); 
        
Dev(1,1) = sum( abs(dFit(1,:)*scale - Data(1,:)) );
Dev(2,1) = sum( abs(dFit(2,:)*scale - Data(2,:)) );
Dev(3,1) = sum( abs(dFit(3,:)*scale - Data(3,:)) );
Dev(4,1) = sum( abs(dFit(4,:)*scale - Data(4,:)) );

[Cost, Factors] = costBalancedInfo(Dev, Vend, Vmax);     

if  isnan(Cost) || Cost == -Inf
        Cost = Inf;
end;

if ~plotFlag
    return
end

 Params = sprintf('%.3g ', X');
 Status = sprintf('Cost=%.5g Dev=%.3g D_1=%.3g D_2=%.3g F=%.3g [%.1e %.1e %.1g]', ...
                    Cost, Factors);
 figure(1);
 for ii = 1 :4 
            subplot(4,1,ii);
            plot(tArray, dFit(ii,:)*scale, tArray, Data(ii,:) );
            if ii == 2               
                title(Params);
            elseif ii == 1
                title(prefix);
            elseif ii == 4
                title(Status);
            end;
            axis tight;
 end;
 drawnow; 

end % function