- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2022-04-14 15:36:40+12:00
- Desc:
- Scaling number of gating variables by new SA
- Permanent Source URI:
- https://models.cellml.org/workspace/6d7/rawfile/c091245cb6b6feb4858b16355c1dac7803d068ca/matlab_parameter_fitting/LCC_bg_parameters.m
clear;
clc;
%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
main_dir = current_dir(1:Idx_backslash(end));
data_dir = ['data' filesep];
code_dir = ['code' filesep];
output_dir = ['output' filesep];
storage_dir = ['storage' filesep];
%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;
N_A = 6.022e23;
x_LCC = 50000/N_A*1e15; % unit fmol
%% Define volumes (unit pL)
W_i = 3.71;
W_e = 0.5062;
%% Import IV parameters
A_geo = 1.1e-5; %1.534e-4; % Unit cm^2
% Load P_GHK (unit pL/s)
load([storage_dir 'LCC_P_GHK_Ca.mat']);
P_Ca = P_GHK;
load([storage_dir 'LCC_P_GHK_K.mat']);
P_K = P_GHK;
%% Import gate transition parameters
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
N_fT = transpose(N_f);
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_cols = sum(stoich_file_data{1}{1} == ',')+1;
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_rT = transpose(N_r);
N = N_r - N_f;
N_T = N_rT - N_fT;
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
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 = [P_Ca/x_LCC; ... % R_GHK_Ca1
P_Ca/x_LCC; ... % R_GHK_Ca2
P_K/x_LCC; ... % R_GHK_K1
P_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 = [P_Ca/x_LCC; ... % R_GHK_Ca1
P_Ca/x_LCC; ... % R_GHK_Ca2
P_K/x_LCC; ... % R_GHK_K1
P_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
k = [kf;kr];
W = [ones(28,1);W_i;W_e;W_i;W_e;ones(12,1)];
lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:28);
K = lambda(29:end);
save('LR_LCC_ch_parameters.mat','params_d','params_f',...
'kappa','K');
%% Checks
N_rref = rref(N);
R_mat = null(N,'r');
K_eq = kf./kr;
zero_est = R_mat'*K_eq;
k_est = exp(M*log(lambdaW));
diff = sum(abs((k-k_est)./k));
