function [jpump] = higgins_pmca_v1(kt,ca_cyt)
%
%function computes steady-state flux for PMCA adaptation of higgins serca
%=> only 1x Ca2+ transported per cycle
%=> is also Electrogenic
%
%this version excludes the electrogenic component for testing
%
%parms are as given by Higgins:
%
%   kt = [k2t k4t kn2t kn4t K1 K3]
%
%the 'tilde' versions:
%
%  
%   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
%
%below from Sage
%Jserc = (Pt*ci*k1*k2*k3*k4 - Pt*ce*kn1*kn2*kn3*kn4) /
% (ce*kn2*kn3*kn4 + ((ce*kn3 + k3 + k4)*k2 + k3*k4 + 
%   (ce*kn3 + k4)*kn2)*ci*k1 + (ce*kn3*kn4 + k3*(k4 + kn4))*k2 + 
%       (ce*kn3*kn4 + k3*(k4 + kn4) + (ce*kn3 + k4 + kn4)*kn2)*kn1)
%
%with 'ci' = cytosolic ca and 'ce' = extracellular; k's as before with K1 = kn1/k1 ,
%K3 = kn3/k3

%
%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.314; 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);

%split K12 and K32 into k1, kn1, k3, kn3 assuming k1 >> kn1, k3 >> kn3
k1 = 1e3;
k3 = 1e3;
kn1 = K12*k1;
kn3 = K32*k3;

%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 terms
temp_num1 = Pt*(ca_cyt*k1*k2*k3*k4 - ca_ext*kn1*kn2*kn3*kn4);
temp_den1 = ca_ext*kn2*kn3*kn4;
temp_den2 = ((ca_ext*kn3 + k3 + k4)*k2 + k3*k4 + (ca_ext*kn3 + k4)*kn2)*ca_cyt*k1;
temp_den3 = (ca_ext*kn3*kn4 + k3*(k4 + kn4))*k2;
temp_den4 = (ca_ext*kn3*kn4 + k3*(k4 + kn4) + (ca_ext*kn3 + k4 + kn4)*kn2)*kn1;

jpump = temp_num1./(temp_den1 + temp_den2 + temp_den3 + temp_den4);