%
%script for computing both Tong pmca model (mm-formulation) and a 2-state
%simplified model of pmca (derived from higgins serca model...)
%

%flags
con_fit = 0;

%set ca concentration range

ca_low = -2;
ca_hi = 1;
ca_n = 100;

ca_cyt_soln = logspace(ca_low,ca_hi,ca_n);

ca_cyt_ss = 0.1162; % from Tong, 2011
ca_ext = 2.5*1000; %mM -> uM from Tong, 2011

%set phosphates here 
% atp_conc = 9040; %uM
% adp_conc = 10;
% p_conc = 520;
atp_conc = 3000; %uM
adp_conc = 10;
p_conc = 3000;
phosph = [atp_conc adp_conc p_conc];
%speed protein conc
sPt = [1 1]; %for normalised speed / concentration
%sPt = [15.1561 10/15.1561]; %higgins values

%computes activity of PMCA pump given cytosolic ca concentration using hill-function rep taken from
%Tong, 2011
%
%   J_pmca = Vmax / (1 + ( K_pmca/ ca_cyt)^n)
%
%Vmax    = 0.00000035*1e6; %{mM_per_msec} -> uM / sec (opencor)
%Vmax = 1e-2;
Vmax    = 1; %normalised
n       = 2;
%testing
n_t2       = 1.5;
Kpmca   = 0.0005*1e3;   %{mM} -> uM

%factor  = 1.65; %debug
%below for matching Shmygol
% factor = 1.45;
% Kpmca = 0.0003*1e3;
%below for matching Tong
factor = 1.0;


%calculate pump activity
jpmca = factor*Vmax./(1 + (Kpmca./ca_cyt_soln).^n);

%pump activity with Hill of 1.5 <-> stoichiometry of 1.5:1 Ca2+:ATP (or
%3:2)
jpmca_n2 = factor*Vmax./(1 + (Kpmca./ca_cyt_soln).^n_t2);

%modify jpmca to hit ss ca cyt at 0.1162 (used by Tong for ca_i_0)
%jpmca = jpmca - 0.0513;

%tong pump
%fig_nob_h = figure;
fig_wb_h = figure;
% set(gcf,'position',[987   708   580   637])
subplot(2,1,1)
% semilogx(ca_cyt_soln,jpmca,'b-','linewidth',2)
% hold on;grid on;

% %compare with higgins serca adaptation
% [jserc_mod, cacyt_ss] = higgins_serc_v1(ca_cyt_soln);
% 
% min_jserc_mod = min(jserc_mod);
% 
% %figure
% %higgins pump - not fit to tong...
% semilogx(ca_cyt_soln,jserc_mod + abs(min_jserc_mod),'k-')


% titlestr = ['cacyt_ss: ' num2str(cacyt_ss)];
% title(titlestr)
%title('PMCA Flux (normalised)')

%xlabel('Ca cyt (in \muM)')

%compute ss ca_cyt for this pump version
%jserc_mod_ss = ca_cyt_soln(abs(jserc_mod) <= 1e-3);

%see how lsqnonlinfit does with the higgins model and fitting to tong 'data'

%parms from Higgins: 
R = 8.314; T = 310;
%dG = -50;
dG = -5e4; %J/mol

%give initial estimates for each
%s = 15.1561; %'speed' - set to give normalised at max ca_cyt above
s = 1;
k2 = 0.6*s; %1/s
%testing
%k2 = 0.55*s; %1/s
k4 = 0.4*s;
%testing
%k4 = 0.45*s;
kn4 = 1.2e-3*s;
kn2 = 0.97*s;

K2 = kn2/k2;
K4 = kn4/k4;

K12 = 0.7; %(umol/L Cyt)^2
K32 = 1.111e-5; %(umol/L ER)^2
%K32 =  6.3634e-07;  %calculated to give ca_cyt_ss = 0.1162 from Tong, 2011 s.t. jpump = 0 at 0.1162...
                    %via K3_req = ca_cyt/(K1*sqrt(K2*K4)*ca_ext) from Eqn. in Higgins, 2006

%k2t and k4t = rates without adp/atp/p (tilde)
%used for non-equilibrium state
k2t = 2e-4;
k4t = 0.4;
kn4t = 4e-8;
%kn2t = kn2/s;
%kn2t = 3.76e-9*(k2t*k4t)/(kn4t*K12*K32); % = exp(dG/RT) = exp(-50/RT)..See Higgins, Eqn 17
kn2t = (exp(dG/(R*T))*k2t*k4t)/(kn4t*K12*K32);

K2t = kn2t/k2t;
K4t = kn4t/k4t;



% K1 = sqrt(K12);
% K3 = sqrt(K32);

%Higgins _defines_ K1^2 = kn1/k1 and K3^2 = kn3/k3 so no sqrt()
K1 = K12;
K3 = K32;

%assemble for initial try - APPARENTLY Higgins _defines_ K1^2 as kn1/n1,
%K3^2 = kn3/k3 - so all below stuff must change...there's no sqrt()
%needed... :-/
%k0 = [k2 k4 kn2 kn4 K1 K3]
k0 = [k2 k4 kn2 kn4 K1 K3];
%k02 = [k2 k4 kn2 kn4 K12 K32];

%compute steady-state ca_i_ss from above k's
ca_cyt_ss_from_k = ca_ext*(K12*K32*(kn2/k2)*(kn4/k4));

%from kt's with phosphates
ca_cyt_ss_from_kt = ca_ext*(K12*K32*(kn2t/(k2*atp_conc))*((kn4*adp_conc*p_conc)/(k4t)));

%k0tilde - for non-equilibrium
k0t = [k2t k4t kn2t kn4t K12 K32];

%compute pump activity at these baseline Higgins parms
jserc_higgins_base_nop = higgins_serc_fit_ktilde_v1(k0t,ca_cyt_soln);
jserc_higgins_base = higgins_serc_fit_ktilde_phosph_s_v1(k0t,phosph,sPt,ca_cyt_soln);

%semilogx(ca_cyt_soln,jserc_higgins_base,'ko')

%hold on; grid on;
%data to fit is the tong pump 
%kfit = lsqcurvefit('higgins_serc_fit_v1',k0,ca_cyt_soln,jpmca);

%compute jpump with this fit suite of k's
%jserc_fit = higgins_serc_fit_v1(kfit,ca_cyt_soln);

%semilogx(ca_cyt_soln,jserc_fit,'m-')

%test same thing but function using the _tilde_ k's, then calculates the
%k's internally
kfit_t = lsqcurvefit('higgins_serc_fit_ktilde_v1',k0t,ca_cyt_soln,jpmca);
%below with passing phosphates & s/pt -- won't work with lsqcurvefit args
%won't work
%kfit_t = lsqcurvefit('higgins_serc_fit_ktilde_phosph_s_v1',k0t,phosph,sPt,ca_cyt_soln,jpmca);

%same fitting but with bounds to prevent negative rrates
lower_bound = 0.001*min(k0t)*ones(1,length(k0t));
upper_bound = 1000*max(k0t)*ones(1,length(k0t));

kfit_bound_t = lsqcurvefit('higgins_serc_fit_ktilde_v1',k0t,ca_cyt_soln,jpmca,lower_bound,upper_bound);

%compute ss cacyt with these k's
%from kt's with phosphates
ca_cyt_ss_from_kfitt = ca_ext*(kfit_t(5)*kfit_t(6)*(kfit_t(3)/(kfit_t(1)*atp_conc))*((kfit_t(4)*adp_conc*p_conc)/(kfit_t(2))));


%compute jpump with this fit suite of k tildes
jserc_t_fit = higgins_serc_fit_ktilde_v1(kfit_t,ca_cyt_soln);
%jserc_t_fit = higgins_serc_fit_ktilde_phosph_s_v1(kfit_t,phosph,sPt,ca_cyt_soln);

jserc_t_bound_fit = higgins_serc_fit_ktilde_v1(kfit_bound_t,ca_cyt_soln);

%semilogx(ca_cyt_soln,jserc_t_fit,'ko')
%hold on;grid on;
% 
% %compute L2norm for this fit
% L2fit_tilde = norm(abs(jserc_t_fit - jpmca))
% 
% legend('Tong PMCA','Higgins SERCA','Fit SERCA','location','best')
% 
% ylabel('Pump flux (normalised)')
% titlestr = ['L2 fit Higgins SERCA:Tong PMCA: ' num2str(L2fit_tilde)];
% title(titlestr);

%compute PMCA conversion of Higgins: the 4-state model as applied to pmca
%with 1 Ca2+ per cycle (and electrogenic...later)
jpmca_higgins_base = higgins_pmca_v1(k0t,ca_cyt_soln);

%next compute for stoichiometry of 1.5:1
jpmca_s1p5_higgins_base = higgins_pmca_s1p5_v1(k0t,ca_cyt_soln);

%modify base pump rate for reverse action ~0.1*max
%jpmca_higgins_base_shift = jpmca_higgins_base - 0.1*max(jpmca_higgins_base);

% figure(fig_nob_h);
% %subplot(2,1,2)
% semilogx(ca_cyt_soln,jpmca,'b-','linewidth',2)
% hold on;grid on;
% semilogx(ca_cyt_soln,jpmca_higgins_base,'m-')
% % title(titlestr)
% title('PMCA Flux (Tong, 2011 normalised)')
% xlabel('Ca cyt (in \muM)')

%figure(fig_wb_h);
%subplot(2,1,2)
semilogx(ca_cyt_soln,jpmca,'b-','linewidth',2)
hold on;grid on;
semilogx(ca_cyt_soln,jpmca_higgins_base,'m-')
% title(titlestr)
title('PMCA Flux (Tong, 2011 normalised, n=2.0; Kinetic n=1.0)')
xlabel('Ca cyt (in \muM)')


% %fit the pmca conversion parms
% kfit_pmca_t = lsqcurvefit('higgins_pmca_v1',k0t,ca_cyt_soln,jpmca);
% 
% %compute jpump with this fit suite of k tildes
% jpmca_t_fit = higgins_pmca_v1(kfit_pmca_t,ca_cyt_soln);

%same but with bounded parms
kfit_pmca_bound_t = lsqcurvefit('higgins_pmca_v1',k0t,ca_cyt_soln,jpmca,lower_bound,upper_bound);

%compute jpump with this fit suite of bounded k tildes
jpmca_t_wb_fit = higgins_pmca_v1(kfit_pmca_bound_t,ca_cyt_soln);

% figure(fig_nob_h);
% semilogx(ca_cyt_soln,jpmca_t_fit,'mo')

%figure(fig_wb_h);
semilogx(ca_cyt_soln,jpmca_t_wb_fit,'mo')

legend('Tong Hill','Kinetic Higgins k','Fit Kinetic')

%compare with stoichiometry of 1.5:1 Ca2+:ATP

%same but with bounded parms
kfit_s1p5_pmca_bound_t = lsqcurvefit('higgins_pmca_s1p5_v1',k0t,ca_cyt_soln,jpmca,lower_bound,upper_bound);

%compute jpump with this fit suite of bounded k tildes
jpmca_t_s1p5_wb_fit = higgins_pmca_s1p5_v1(kfit_s1p5_pmca_bound_t,ca_cyt_soln);

subplot(2,1,2)
semilogx(ca_cyt_soln,jpmca_n2,'b-','linewidth',2)
hold on;grid on;
semilogx(ca_cyt_soln,jpmca_s1p5_higgins_base,'m-')
semilogx(ca_cyt_soln,jpmca_t_s1p5_wb_fit,'mo')
% title(titlestr)
title('PMCA Flux (Tong, 2011 normalised, n=1.5; Kinetic n=1.5)')
xlabel('Ca cyt (in \muM)')




%compute ss ca cyt for this pump version
%jserc_fit_ss = ca_cyt_soln(abs(jserc_fit) <= 4e-4);

%legend('Tong PMCA','4-state Higgins PMCA','Fit 4-state')
% textstr = ['Cacyt SS: Fit: ' num2str(jserc_fit_ss) ' 4-State: ' num2str(jserc_mod_ss)];
% text(0.0333,-0.1164,textstr)
% 
% %adjust K32 - calc to give ca_cyt_ss = 0.1162...?
% K3_adj = ca_cyt_ss/(kfit(5)*sqrt((kfit(3)/kfit(1))*kfit(4)/kfit(2))*ca_ext);
% 
% %compute pump with adjusted K3
% jserc_adj = higgins_serc_fit_v1([kfit(1) kfit(2) kfit(3) kfit(4) kfit(5) K3_adj],ca_cyt_soln);
% 
% %semilogx(ca_cyt_soln,jserc_adj,'c-')
% 
%
%try selecting points for fit of a sigmoidal:
%use the 1/2 maximal of tong pump, the max at 10 uM, but the equilibrium at ca_cyt_ss
%
%use the 'base higgins' pump to get min pump rate


% pump_ca_cyt = [min(ca_cyt_soln) ca_cyt_ss Kpmca max(ca_cyt_soln)];
% %pump_j = [min(jserc_mod) 0 0.5 max(jpmca)];
% %pump_j = [min(jserc_higgins_base) 0 0.5 max(jpmca)];
% pump_j = [min(jpmca_higgins_base_shift) 0 0.5 max(jpmca)];
% 
% %assemble parms for fit
% ksig = [Kpmca 2];
% 
% ksig_fit = lsqcurvefit('sigmoidal',ksig,pump_ca_cyt,pump_j);
% 
% %calculate pump activity with these ksig_fit
% jpmca_sig = factor*Vmax./(1 + (ksig_fit(1)./ca_cyt_soln).^ksig_fit(2));
% 
% semilogx(ca_cyt_soln,jpmca_sig,'g-')

%
% %next try with a constraint on the k's such that they satisfy thermal condition
% 
% %since we don't have k1, kn1, k3, kn3, we'll base values on assumption k1
% %>> k2, k4; and k3 >> k2, k4 then compute kn1, kn3 from there
% k1 = 10*max([k2 k4]);
% k3 = 10*max([k2 k4]);
% 
% kn1 = K1*k1;
% kn3 = K3*k3;
% 
% %define vect of these - note the k2's and k4's are modified by 's'
% k_vect = [k1 k2 k3 k4 kn1 kn2 kn3 kn4];
% %below using tildes
% k_vect_noneq = [k1 k2t k3 k4t kn1 kn2t kn3 kn4t];
% 
% %test these k's for therm constraints
% [kvect_con, kvect_eq] = higgins_therm_constraint(k_vect)
% 
% [k_vect_noneq_con, k_vect_noneq_eq] = higgins_therm_constraint(k_vect_noneq)
% 
% %compare with function using split k's -- presumed k1 prop to kn1 etc
% k_split_fit = lsqcurvefit('higgins_serc_fit_splitK_v1',k_vect,ca_cyt_soln,jpmca);
% 
% %compute jpump with this fit suite of k's
% %jserc_split_fit = higgins_serc_fit_v1(k_split_fit,ca_cyt_soln);
% jserc_split_fit = higgins_serc_fit_splitK_v1(k_split_fit,ca_cyt_soln);
% 
% semilogx(ca_cyt_soln,jserc_split_fit,'mo')
% 
% %legend('Tong PMCA','4-state Higgins','Fit 4-state','Split-K')
% %legend('Tong PMCA','4-state Higgins','Fit 4-state (Split-K)')
% 
% %quick test higgins / tong pump objective function
% 
% test_r0 = higgins_tong_serc_fit_v1(k_vect)
% 
% %compare with fit version
% test_r_fit = higgins_tong_serc_fit_Ks_v1(kfit)
% 
% test_r_presplit_fit = higgins_tong_serc_fit_v1(k_split_fit)
% 
% %compare with fit version - but with split K's
% %need to split out K1, K3
% kn1_fit = [k1 k2 k3 k4 -1 kn2 -1 kn4];
% %use the fit K1, K3
% kn1_fit = kfit(5)*k1;
% kn3_fit = kfit(6)*k3;
% 
% kfit_split_postfit = [k1 k2 k3 k4 kn1_fit kn2 kn3_fit kn4];
% 
% test_r_fit_split_postfit = higgins_tong_serc_fit_v1(kfit_split_postfit)
% 
% %quick test higgins thermal constraint
% [test_con, test_eq] = higgins_therm_constraint(kfit_split_postfit)
% 
% %zero concentration for reference
% 
% semilogx(ca_cyt_soln,zeros(1,ca_n),'r-')
% 
% % legend('Tong PMCA','4-state Higgins','Fit 4-state (Split-K)', 'Jpmca = 0','location','best')
% legend('Tong PMCA','4-state Higgins','Fit: lsq','Jpmca = 0','location','best')

if con_fit

    %
    %fit to tong pmca using the nonequilibrium kvect with the tilde k2, k4's
    %
    
    %test the therm constraint function
    [def_con, def_eq] = higgins_therm_const_tilde(k0t);
    
    
%     %compare with function using split k's -- presumed k1 prop to kn1 etc
%     k_noneq_split_fit = lsqcurvefit('higgins_serc_fit_splitK_v1',k_vect_noneq,ca_cyt_soln,jpmca);
% 
%     fprintf('fmincon test...\n');
%     % 
%     % %Solve a Constrained Nonlinear Problem
%     % 
%     % %setup options etc
%     options = optimoptions(@fmincon,...
%         'Display','iter','Algorithm','interior-point','OptimalityTolerance',1e-4,...
%         'FiniteDifferenceType','Central');
%     % % 
%     % % %solve constrained linear problem with serca fit to tong pmca
%     % [kfit_split_fit_postconst,test_r_fit_split_postconst] = fmincon(@higgins_tong_serc_fit_v1,k_split_fit,...
%     %     [],[],[],[],[],[],@higgins_therm_constraint,options);
%     % % [x,fval] = fmincon(@higgins_tong_serc_fit_v1,k_vect_noneq,...
%     % %     [],[],[],[],[],[],[],options);
%     % 
%     % %test result
%     % jserc_split_fit_const = higgins_serc_fit_splitK_v1(kfit_split_fit_postconst,ca_cyt_soln);
%     % %test_r_fit_split_postfit_postconst = higgins_tong_serc_fit_v1(kfit_split_fit_postconst)
%     % 
%     subplot(2,1,2)
%     %tong pump
%     %semilogx(ca_cyt_soln,jpmca,'b-')
%     hold on;grid on;
%     %higgins pump - not fit to tong...
%     semilogx(ca_cyt_soln,jserc_mod,'k-')
% 
%     % semilogx(ca_cyt_soln,jserc_split_fit_const,'mx')
% 
%     %use original kvect before fitting to tong
%     % %solve constrained linear problem with serca fit to tong pmca
%     %set lb to zero vector prevent k's going negative...
%     lb = zeros(1,length(k_vect));
%     [kvect_postconst,test_r_postconst] = fmincon(@higgins_tong_serc_fit_v1,k_vect,...
%         [],[],[],[],lb,[],@higgins_therm_constraint,options);
% 
%     jserc_postconst = higgins_serc_fit_splitK_v1(kvect_postconst,ca_cyt_soln);
% 
%     %semilogx(ca_cyt_soln,jserc_postconst,'go')
% 
%     %next try testing the whole thing - can fmincon fit to original higgins pump (so no need to fit) and
%     %satisfy constraint
%     [kvect_test_postconst,test_serc_r_postconst] = fmincon(@higgins_test_serc_fit_v1,k_vect,...
%         [],[],[],[],lb,[],@higgins_therm_constraint,options);
% 
% 
%     jserc_test_postconst = higgins_serc_fit_splitK_v1(kvect_test_postconst,ca_cyt_soln);
% 
%     semilogx(ca_cyt_soln,jserc_test_postconst,'rd')
% 
%     %legend('Tong','Higgins','Fit Tong Const','Fit Higgins Const','location','best')
% 
%     % %zero concentration for reference
%     % semilogx(ca_cyt_soln,zeros(1,ca_n),'r-')
%     % 
%     % % legend('Tong PMCA','4-state Higgins','Fit 4-state (Split-K)', 'Jpmca = 0','location','best')
%     % legend('Tong PMCA','4-state Higgins','Fit: lsq', 'Fig: fmin','Jpmca = 0','location','best')
    
end