- Author:
- shawnm <shawnm@rocko.local>
- Date:
- 2025-06-03 15:47:26+12:00
- Desc:
- updated readme.rst
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/d28/rawfile/ee7c1be211c09cb9b1639edfc10505b89eae3510/parameter_fit_n_find/pmca_fun1p1.m
%
%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
