function [jpump] = higgins_serc_fit_ktilde_v1(kt,ca_cyt)
%
%function for computing pump rate for Higgins serca model; this takes
%cytosolic calcium as well as 6 kinetic rates for k2t, k4t, kn2t, kn4t and K1 K3 from 
%multi-state model of Higgins, 2006:
%
%   kt = [k2t k4t kn2t kn4t K1 K3]
%
%Parms are taken as the "tilde" versions
%for k2 and k4 as defined by Higgins:
%
%   k2 = k2t*[ATP]*s
%   k4 = k4t*s
%   kn2 prop to exp(DeltaG / RT)  <- gibbs free energy constraint
%                           (from K12*K2t*K32*K4t = exp(DeltaG / RT), see
%                           Higgins Eqn (17)
%   kn2 = exp(DeltaG / RT)*kn2t / (K4t*K12*K32) * s
%       = 3.76e-9 * kn2t / (K4t*K12*K32) * s
%   kn4 = kn4t*[ADP]*[P]*s
%
%with 's' speed, ATP, ADP, P all set internally here.
%
%K12, and K32 are the _nonsquared_ K1 = kn1/k1, K3 = kn3/k3 (should be
%constant)
%
%with above conversions from the tilde k's then
%
%It computes pump rate as given by:
%   
%   Jserc = 2*(-kn2*K3^2*Casr^2*kn4*Kn1^2 + k2*Cacyt^2*k4)*(Pc || Psd)
%           ---------------------------------------------------------------
%           (Casr*Cacyt)^2*K3^2*(k2 + kn2) + Cacyt^2*(k4 + k2) +
%           Casr^2*K1^2*K3^2*(kn2 + kn4) + K1^2*(k4 + kn4)
%


% %compute tong pump activity - NOT THIS VERSION
% %set ca concentration range
% 
% ca_low = -2;
% ca_hi = 1;
% ca_n = 100;
% 
% ca_cyt_soln = logspace(ca_low,ca_hi,ca_n);
% 
% %calculate pump activity via Tong formulation
% factor = 1.0;
% Vmax = 1.0;
% n       = 2;
% Kpmca   = 0.0005*1e3;   %{mM} -> uM
% jpmca = factor*Vmax./(1 + (Kpmca./ca_cyt_soln).^n);

%
%set parms here  
%
%below speed for fitting pump data should not be normalised
%s = 15.1561; %'speed' - set to give normalised at max ca_cyt above
s = 1;

%pump concentration
%Pt = 10.0;  %umol/L Cyt
%Pt = 10/s;  %umol/L Cyt
Pt = 1.0;

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

%ATP/ADP/P0 concentrations (from Higgins)
%atp_conc = 3000; %uM
%adp_conc = 10;
%p_conc = 3000;
%below for 0.1162 uM cai ss
 atp_conc = 9040; %uM
 adp_conc = 10;
 p_conc = 520;

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

%pull out final k's K12, K32 (NOT squared)
K12 = kt(5);
K32 = kt(6);

%convert from passed tilde k's
k2 = kt(1)*atp_conc*s;
k4 = kt(2)*s;
%below kn2 simply = kn2t
%kn2 = kt(3);
%kn2 = s*(exp(dG/(R*T))*kt(3)*kt(2))/(kt(4)*K12*K32)
%kn2 = s*(exp(dG/R*T)*kt(3)*kt(2))/(kt(4)*K12*K32);
%below compute kn2t from thermal constraint and use as kn2
kn2t = s*(exp(dG/R*T)*kt(1)*kt(2))/(kt(4)*K12*K32);
kn2 = kn2t;
kn4 = kt(4)*adp_conc*p_conc*s;



%compute actual serca pump flux
temp_num1 = -kn2*K32*ca_ext^2*kn4*K12;
temp_num2 = k2*k4*ca_cyt.^2;

temp_den1 = (ca_cyt*ca_ext).^2*K32*(k2 + kn2);
temp_den2 = (k4 + k2)*ca_cyt.^2;
temp_den3 = K12*K32*(kn2 + kn4)*ca_ext^2;
temp_den4 = K12*(k4 + kn4);

%assemble terms
jpump = 2*Pt*(temp_num1 + temp_num2)./(temp_den1 + temp_den2 + temp_den3 + temp_den4);

%shift
%jpump = jpump + 0.036091160584956;