- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-10-22 15:42:06+13:00
- Desc:
- Adding parameter finding in Python and tidied cellml code
- Permanent Source URI:
- https://models.cellml.org/workspace/6d7/rawfile/fc4ffea6efa99d4c793d0b6d2b8ea4559540734b/matlab/LCConly_cardiac_AP_bg_parameters.m
% Based on Pan cardiac action potential BG modelling 2018
% only LCC channel is retained.
clear;
clc;
close all;
%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
data_dir = [current_dir filesep 'data' filesep];
code_dir = [current_dir filesep 'code' filesep];
output_dir = [current_dir filesep 'output' filesep];
storage_dir = [current_dir filesep 'storage' filesep];
%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;
%
N_A = 6.022e23;
A_cap = 1.534e-4; % Unit cm^2
x_LCC = 50000/N_A*1e15; % unit fmol
%% Define volumes (unit pL)
W_i = 38;
W_e = 5.182;
% %% Import IV parameters
% L-type calcium channel
% Load P_GHK (unit pL/s)
load([storage_dir 'LCC_P_GHK_Ca.mat']);
P_LCC_Ca = P_GHK; % Unit pL/s
load([storage_dir 'LCC_P_GHK_K.mat']);
P_LCC_K = P_GHK; % Unit pL/s
%
load([storage_dir 'LCC_d_parameters.mat']);
params_d = params_vec;
load([storage_dir 'LCC_f_parameters.mat']);
params_f = params_vec;
%% Load stoichiometric matrices
stoich_path = [data_dir 'LCC_forward_matrix.txt'];
stoich_file_id = fopen(stoich_path,'r');
stoich_file_data = textscan(stoich_file_id,'%s','delimiter','\n');
fclose(stoich_file_id);
num_rows = length(stoich_file_data{1});
num_cols = sum(stoich_file_data{1}{1} == ',')+1;
N_f = zeros(num_rows,num_cols);
for i_row = 1:num_rows
line_str = stoich_file_data{1}{i_row};
line_split = regexp(line_str,',','split');
for i_col = 1:num_cols
N_f(i_row,i_col) = str2double(line_split{i_col});
end
end
% Load reverse matrix
stoich_path = [data_dir 'LCC_reverse_matrix.txt'];
stoich_file_id = fopen(stoich_path,'r');
stoich_file_data = textscan(stoich_file_id,'%s','delimiter','\n');
fclose(stoich_file_id);
num_rows = length(stoich_file_data{1}); % num of kappa + K
num_cols = sum(stoich_file_data{1}{1} == ',')+1; % num of reactions
N_r = zeros(num_rows,num_cols);
for i_row = 1:num_rows
line_str = stoich_file_data{1}{i_row};
line_split = regexp(line_str,',','split');
for i_col = 1:num_cols
N_r(i_row,i_col) = str2double(line_split{i_col});
end
end
N_fT = transpose(N_f);
N_rT = transpose(N_r);
N = N_r - N_f;
N_T = N_rT - N_fT;
num_cols = size(N,2);
I = eye(num_cols);
M = [I N_fT; I N_rT];
M_rref = rref(M);
M_T = transpose(M);
M_T_rref = rref(M_T);
%% Set up the vectors
% L-type calcium channel
alpha_d0_bg = params_d(1)*1e3; % unit s^-1
beta_d0_bg = params_d(3)*1e3; % unit s^-1
alpha_f1_0_bg = params_f(1)*1e3; % unit s^-1
beta_f1_0_bg = params_f(3)*1e3; % unit s^-1
alpha_f2_0_bg = params_f(5)*1e3; % unit s^-1
beta_f2_0_bg = params_f(7)*1e3; % unit s^-1
K_f3_0 = beta_f1_0_bg*alpha_f2_0_bg/alpha_f1_0_bg/beta_f2_0_bg; % dimensionless
rate_f3 = 1e5; % unit s^-1
KmCa = 0.6e-3; % Unit mM
rate_fCa = 1e5; % unit s^-1
kf_LCC = [P_LCC_Ca/x_LCC; ... % R_GHK_Ca1
P_LCC_Ca/x_LCC; ... % R_GHK_Ca2
P_LCC_K/x_LCC; ... % R_GHK_K1
P_LCC_K/x_LCC; ... % R_GHK_K2
alpha_d0_bg; ... % Rd000
alpha_d0_bg; ... % Rd010
alpha_d0_bg; ... % Rd020
alpha_d0_bg; ... % Rd001
alpha_d0_bg; ... % Rd011
alpha_d0_bg; ... % Rd021
alpha_f1_0_bg; ... % Rf1_000
alpha_f1_0_bg; ... % Rf1_100
alpha_f1_0_bg; ... % Rf1_001
alpha_f1_0_bg; ... % Rf1_101
alpha_f2_0_bg; ... % Rf2_000
alpha_f2_0_bg; ... % Rf2_100
alpha_f2_0_bg; ... % Rf2_001
alpha_f2_0_bg; ... % Rf2_101
K_f3_0*rate_f3; ... % Rf3_010
K_f3_0*rate_f3; ... % Rf3_110
K_f3_0*rate_f3; ... % Rf3_011
K_f3_0*rate_f3; ... % Rf3_111
rate_fCa; ... % RfCa000
rate_fCa; ... % RfCa100
rate_fCa; ... % RfCa010
rate_fCa; ... % RfCa110
rate_fCa; ... % RfCa020
rate_fCa]; % RfCa120
kr_LCC = [P_LCC_Ca/x_LCC; ... % R_GHK_Ca1
P_LCC_Ca/x_LCC; ... % R_GHK_Ca2
P_LCC_K/x_LCC; ... % R_GHK_K1
P_LCC_K/x_LCC; ... % R_GHK_K2
beta_d0_bg; ... % Rd000
beta_d0_bg; ... % Rd010
beta_d0_bg; ... % Rd020
beta_d0_bg; ... % Rd001
beta_d0_bg; ... % Rd011
beta_d0_bg; ... % Rd021
beta_f1_0_bg; ... % Rf1_000
beta_f1_0_bg; ... % Rf1_100
beta_f1_0_bg; ... % Rf1_001
beta_f1_0_bg; ... % Rf1_101
beta_f2_0_bg; ... % Rf2_000
beta_f2_0_bg; ... % Rf2_100
beta_f2_0_bg; ... % Rf2_001
beta_f2_0_bg; ... % Rf2_101
rate_f3; ... % Rf3_010
rate_f3; ... % Rf3_110
rate_f3; ... % Rf3_011
rate_f3; ... % Rf3_111
rate_fCa/KmCa; ... % RfCa000
rate_fCa/KmCa; ... % RfCa100
rate_fCa/KmCa; ... % RfCa010
rate_fCa/KmCa; ... % RfCa110
rate_fCa/KmCa; ... % RfCa020
rate_fCa/KmCa]; % RfCa120
% Overall
k = [kf_LCC;kr_LCC];
W = [ones(size(N,2),1);W_i;W_e;W_i;W_e;ones(12,1)]; % kappa Cai Cae Ki Ke [gates]
lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);
save([output_dir 'LCC_parameters.mat'],'kappa','K','params_f');
%% Checks
N_rref = rref(N);
R_mat = null(N,'r');
K_eq = kf_LCC./kr_LCC;
zero_est = R_mat'*K_eq;
k_est = exp(M*log(lambdaW));
diff = sum(abs((k-k_est)./k));
%% output for cellML
if false
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);
fID = fopen([data_dir 'rxnID.txt'], 'r');
rxnID = textscan(fID,'%s', 'delimiter','\n');
fclose(fID);
fID = fopen([data_dir 'Kname.txt'], 'r');
Kname = textscan(fID,'%s');
fclose(fID);
for j = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{1}{j},kappa(j));
end
for j = 1:length(K)
fprintf('var K_%s: per_fmol {init: %g};\n',Kname{1}{j},K(j));
end
if false
qname = {'L','LR','LRG','LR_BARK','LRG_BARK','R','RG','B1d','B1p','PKACI_B1p',...
'B1tot','PKACI','BARK','G','GaT','GBy','GaD','GsGDP','Pi'};
qval = [0.0001672850,0.0000016970,0.0001046605,0.0000016970,0.0001046605,0.0002197247,...
0.0000254602,0.0000000000,...
0.0000219260,.0000219260,0.0004579000,0.0022120882,0.0228000000,0.1453052187,...
0.0009519000,0.0009762200,0.0000244948,0.0000000000,570]';
for j=1:length(qname)
fprintf('var q_%s: fmol {init: %g};\n',qname{j}, qval(j));
end
end
for j=1:length(K)
fprintf('var q_%s: fmol {init: 1e-13};\n',Kname{1}{j});
end
for j=1:length(kappa)
fprintf('var v%s: fmol_per_sec;\n',rxnID{1}{j})
end
for j=1:length(K)
fprintf('var mu_%s: J_per_mol;\n',Kname{1}{j});
end
for j=1:length(K)
fprintf('mu_%s = R*T*ln(K_%s*q_%s);\n',Kname{1}{j},Kname{1}{j},Kname{1}{j});
end
for j=1:length(kappa)
fprintf('v%s = kappa_%s*exp(mu_a/(R*T));\n',rxnID{1}{j}, rxnID{1}{j})
end
for j=1:length(K)
fprintf('ode(q_%s, time) = p;\n',Kname{1}{j});
end
end
