- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2022-03-03 09:58:16+13:00
- Desc:
- Changing Xi to X_i, and adding python param finder code
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/82c/rawfile/73d6bc1016edb1be3ceadf6cccc4aa05ef251e58/BG_fit_parameters_CC/compare_PSO_alpha_beta_kappa.m
% check fit of krkf output (alpha_fit) for kf zf kr zr by comparing against
% kinetic (alpha)
clear;
V = (-120:60)'; % unit mV
gate_state_name = {'C3','c2','c1','O','IF','IS'};
R = 8.314;
T = 310;
RT = R*T;
F = 96500;
vol = 2.58E-17; % L
Ko = 4.5;
P_O_Na = 1.78e-6;
P_C1 = 1.119e-4;
P_C2 = 0.0071021;
P_C3 = 0.8294071;
P_IF = 0.10199797;
P_IS = 0.0613826;
alpha = 55.5*exp(0.05547153*(V-12));
beta = 2.357*exp(-0.036588*V);
alpha_in = 2.172*ones(size(V));
beta_in = 1.077*ones(size(V));
alpha_alpha = 65.5*exp(0.05547153*(V-36));
beta_beta = 2.9357*exp(-0.02158*V);
beta_i = 439*exp(-0.02352*(V+25))*4.5/Ko;
alpha_i = 656*exp(0.000942*V).*((4.5/Ko).^ 0.3);
mu = alpha_i.*beta_beta.*alpha_alpha./(alpha_alpha.*beta_i);
alpha = [alpha alpha_in alpha_alpha alpha_i alpha_alpha];
beta = [beta beta_in beta_beta beta_i mu];
% some BG parameters: from fit_alpha_beta_kappa_constant_z.m
krkf = [68.438114 2.1524652
2.172 1.077
21.333536 0.55120363
638.49363 55.998552
21.333536 9414.4161
];
krkf = [krkf(:,1) [1;0;1;0;1] krkf(:,2) [-1;0;-1;-1;4]];
%% test kf, kr, zf, zr for fits against alpha and beta
if true
% % BG
% zf_exponent = exp(bsxfun(@times, krkf(:,2)', 1e-3*F*V/RT));
% zr_exponent = exp(bsxfun(@times, krkf(:,4)', 1e-3*F*V/RT));
% alpha_fit = bsxfun(@times, krkf(:,1)', zf_exponent);
% beta_fit = bsxfun(@times, krkf(:,3)', zr_exponent);
for j = 1:5
zf_exponent = exp(krkf(j,2) * 1e-3*F*V/RT);
zr_exponent = exp(krkf(j,4) * 1e-3*F*V/RT);
alpha_fit(:,j) = krkf(j,1) * zf_exponent;
beta_fit(:,j) = krkf(j,3) * zr_exponent;
end
% figure, subplot(1,2,1)
% plot(V, alpha, '.');
% title('alpha_{OG}'); legend('1','2','3','4','5','6');
% subplot(1,2,2)
% plot(V, alpha_fit, '.'); ylim([0 7e4]);
% title('alpha_{fit}'); legend('1','2','3','4','5','6');
%
% figure, subplot(1,2,1)
% plot(V, beta, '.'); legend('1','2','3','4','5','6');
% title('beta_{OG}');
% subplot(1,2,2)
% plot(V, beta_fit, '.'); ylim([0 14e4]); legend('1','2','3','4','5','6');
% title('beta_{fit}');
figure,
for ik = 1:5
subplot(2,3,ik)
plot(V, alpha(:,ik), V, alpha_fit(:,ik));
legend('OG','fit');
xlabel('V'); ylabel('rate'); title(num2str(ik));
end
suptitle('alpha');
figure,
for ik = 1:5
subplot(2,3,ik)
plot(V, beta(:,ik), V, beta_fit(:,ik));
legend('OG','fit');
xlabel('V'); ylabel('rate'); title(num2str(ik));
end
suptitle('beta');
% % plot Open state only
% figure,
% subplot(1,2,1)
% semilogy(V, alpha(:,4), '.', V, alpha_fit(:,4),'.');
% title('alpha P_{Open}'); legend('OG','fit'); xlabel('V');
%
% subplot(1,2,2)
% plot(V, beta(:,4), '.', V, beta_fit(:,4),'.');
% title('beta P_{Open}'); legend('OG','fit'); xlabel('V');
end
%% test v (dPdt) using kappa, K, zf, zr.
% BEWARE: not comparing like and like (OG kinetic: dim'less but BG: mol/s)
if false
alpha = alpha'; beta = beta';
% declare K_c, q_Sc
xgg = 2e-19; % mol
% q_sC3 = 1e-9;
% q_sC2 = 1e-11;
% q_sC1 = 1e-11;
% q_sO = 1e-11;
% q_sIF = 1e-11;
% q_sIS = 1e-11;
q_sC3 = xgg*P_C3;
q_sC2 = xgg*P_C2;
q_sC1 = xgg*P_C1;
q_sO = xgg*P_O_Na;
q_sIF = xgg*P_IF;
q_sIS = xgg*P_IS;
q0 = [q_sC3 q_sC2 q_sC1 q_sO q_sIF q_sIS];
% one K per gate variable
% K [=] per_(f)mol
K_C3 = 535.0649335;
K_C2 = 9.771415327;
K_C1 = 0.268890314;
K_O = 579.0861639;
K_IF = 116.818253;
K_IS = 4.82E-07;
% one kappa per reaction
% kappa [=] (f)mol_per_sec
kappa_rC3 = 19.61852542;
kappa_rC2 = 951.5336816;
kappa_rC1 = 139.9921881;
kappa_rC1IF = 3.17E-04;
kappa_rO = 0.07040671;
kappa_rIF = 0.785664583;
zf = krkf(:,2);
zr = krkf(:,4);
% state variable potentials
mu_sC3 = RT*log(K_C3*q_sC3);
mu_sC2 = RT*log(K_C2*q_sC2);
mu_sC1 = RT*log(K_C1*q_sC1);
mu_sO = RT*log(K_O*q_sO);
mu_sIF = RT*log(K_IF*q_sIF);
mu_sIS = RT*log(K_IS*q_sIS);
Volt = V*1e-3;
Af_rC3 = mu_sC3+zf(1)*F*Volt;
Ar_rC3 = mu_sC2+zr(1)*F*Volt;
Af_rC2 = mu_sC2+zf(2)*F*Volt;
Ar_rC2 = mu_sC1+zr(2)*F*Volt;
Af_rC1 = mu_sC1+zf(3)*F*Volt;
Ar_rC1 = mu_sO+zr(3)*F*Volt;
Ar_rC1IF = mu_sIF+zr(4)*F*Volt;
Af_rC1IF = mu_sC1+zf(4)*F*Volt;
Af_rO = mu_sO+zf(5)*F*Volt;
Ar_rO = mu_sIF+zr(5)*F*Volt;
Af_rIF = mu_sIF+zf(6)*F*Volt;
Ar_rIF = mu_sIS+zr(6)*F*Volt;
v_rC3 = kappa_rC3*(exp(Af_rC3/RT)-exp(Ar_rC3/RT));
v_rC2 = kappa_rC2*(exp(Af_rC2/RT)-exp(Ar_rC2/RT));
v_rC1 = kappa_rC1*(exp(Af_rC1/RT)-exp(Ar_rC1/RT));
v_rC1IF = kappa_rC1IF*(exp(Af_rC1IF/RT)-exp(Ar_rC1IF/RT));
v_rO = kappa_rO*(exp(Af_rO/RT)-exp(Ar_rO/RT));
v_rIF = kappa_rIF*(exp(Af_rIF/RT)-exp(Ar_rIF/RT));
% state variable velocities [=] mol/s
%v = [ sc3 sc2 sc1 o if is]
v_bg(:,1) = -v_rC3;
v_bg(:,2) = v_rC3-v_rC2;
v_bg(:,3) = v_rC2-v_rC1-v_rC1IF;
v_bg(:,4)= v_rC1-v_rO;
v_bg(:,5) = v_rC1IF+v_rO-v_rIF;
v_bg(:,6) = v_rIF;
% v_bg = v_bg * 1e3 / vol; % unit: mM/s
v_bg = 1e-3*v_bg ./ q0; % unit: 1/s
% plot rates
figure,
for i = 1:6
subplot(2,3,i)
plot(V, v_CC(i,:),V, v_bg(:,i));
xlabel('V'); ylabel('rate mM/s'); title('CC');
legend('CC','BG');
title(gate_state_name{i});
end
end
