- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2022-01-28 09:32:45+13:00
- Desc:
- Adding combined crossbridge_TRPN module
- Permanent Source URI:
- https://models.cellml.org/workspace/75f/rawfile/3bc405a7c9215edb27f452c85b8b8e939dd50dbc/composite_whole_B1AR/MATLAB/find_BG_parameters_whole.m
% find bond-graph parameters whole whole cell with multiple modules
% present:
% cAMP; LRGbinding; B1AR; PKA; PLB; Inhib1; GsProtein;
% based on Pan cardiac AP 2018 code
clear;
% clc;
close all;
%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
parent_dir = current_dir(1:Idx_backslash(end-1));
main_dir = current_dir(1:Idx_backslash(end));
output_dir = [current_dir filesep 'output' filesep];
%% Define constants
R = 8.314; % unit J/mol/K
T = 310;
F = 96485;
cNao_st = 140; % unit mM
cNai_st = 10; % unit mM
cKo_st = 5.4;
cKi_st = 145;
N_A = 6.022e23;
A_cap = 1.534e-4; % Unit cm^2
%% Define volumes (unit pL)
V_myo = 34.4;
V_e = 5.182; % external volume
V_SR = V_myo*0.035; % SR volume
V_di = V_myo*0.0539; % diadic volume
V = struct();
V.V_myo = V_myo;
V.V_SR = V_SR;
V.V_di = V_di;
V.V_e = V_e;
%% Load stoichiometric matrices and kinetic rate constants
array_names = {'cAMP'; 'LRGbinding'; 'beta1AR'; 'PKA'; 'PLB'; 'Inhib1'; 'GsProtein'};
num_subsystems = length(array_names);
sys_struct = cell(num_subsystems,1);
rxnIDs = {};
Knames = {};
keepKname = {};
for i_system = 1:num_subsystems
sys_name = array_names{i_system};
sys_dir = [main_dir sys_name '\MATLAB\'];
cd(sys_dir)
sys_struct{i_system}.name = sys_name;
forward_mat_path = ['data\all_forward_matrix.txt'];
reverse_mat_path = ['data\all_reverse_matrix.txt'];
N_f = load_matrix(forward_mat_path);
N_r = load_matrix(reverse_mat_path);
sys_struct{i_system}.N_f = N_f;
sys_struct{i_system}.N_r = N_r;
% disp(array_names{i_system});
dims = struct();
dims.num_rows = size(N_f,1);
dims.num_cols = size(N_f,2);
I = eye(dims.num_cols);
M = [I transpose(N_f); I transpose(N_r)];
[k_kinetic, N_cT, K_C, W] = kinetic_parameters(M, true, dims, V);
sys_struct{i_system}.kfkr = k_kinetic;
sys_struct{i_system}.Kc = K_C;
sys_struct{i_system}.N_cT = N_cT;
% sys_struct{i_system}.W = W(dims.num_cols+1:end);
fID = fopen('data\rxnID.txt', 'r');
rxnID = textscan(fID,'%s', 'delimiter','\n');
fclose(fID);
rxnIDs = [rxnIDs; rxnID{1}];
fID = fopen('data\Kname.txt', 'r');
Kname = textscan(fID,'%s', 'delimiter','\n');
fclose(fID);
Knames = [Knames; Kname{1}];
end
cd(current_dir)
% relations between submodule to whole module
sys_struct{1}.I_vec = 1:15; % cAMP
sys_struct{2}.I_vec = 16:21; % LRGbinding
sys_struct{3}.I_vec = [19 21 22:29]; % B1AR
sys_struct{4}.I_vec = [2 30:42]; % PKA
sys_struct{5}.I_vec = [43 28 44:49]; % PLB
sys_struct{6}.I_vec = [28 48 50:53]; % Inhib1
sys_struct{7}.I_vec = [17 18 20 19 21 14 54:55]; % GsProtein
num_rows = sys_struct{size(sys_struct,1)}.I_vec(end);
N_f = [];
N_r = [];
for i_system = 1:num_subsystems
% disp(array_names{i_system});
T = calcT(sys_struct{i_system}.I_vec,num_rows);
sys_struct{i_system}.T = T;
N_f = [N_f T*sys_struct{i_system}.N_f];
N_r = [N_r T*sys_struct{i_system}.N_r];
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);
%% Set up the vectors for kinetic rate constants
% cAMP; LRGbinding; B1AR; PKA; PLB; Inhib1; GsProtein
kf = [];
kr = [];
W = ones(size(N,2),1);
for i_system=1:num_subsystems
nrx = length(sys_struct{i_system}.kfkr)/2;
% disp(array_names{i_system});
% disp(nrx)
kf = [kf, sys_struct{i_system}.kfkr(1:nrx)'];
kr = [kr, sys_struct{i_system}.kfkr(nrx+1:end)'];
% W = [W; sys_struct{i_system}.W];
end
k_kinetic = [kf, kr]';
W = [ones(size(N,2),1);V_myo*ones(15,1);V_myo*ones(14,1);V_myo*ones(13,1);
V_myo*ones(7,1);V_myo*ones(4,1);V_myo*ones(2,1)];
lambdaW = exp(pinv(M)*log(k_kinetic));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);
save([output_dir 'whole_saucerman03_BG_parameters.mat'],'kappa','K','k_kinetic');
%% Checks
N_rref = rref(N);
R_mat = null(N,'r');
k_est = exp(M*log(lambdaW));
diff = (k_kinetic - k_est)./k_kinetic;
error = sum(abs(diff));
K_eq = kf./kr;
zero_est = R_mat'*K_eq';
%% print outputs
for ik=1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g, pub: out};\n',rxnIDs{ik},kappa(ik));
end
Kunique = {};
for ik=1:length(Knames)
if ~any(strcmp(Kunique,Knames{ik}))
Kunique = [Kunique; Knames{ik}];
end
end
for ik = 1:length(Kunique)
fprintf('var K_%s: per_fmol {init: %g, pub: out};\n',Kunique{ik},K(ik));
end
%% FUNCTION SCRIPTS
function matrix = load_matrix(stoich_path)
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;
matrix = 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
matrix(i_row,i_col) = str2double(line_split{i_col});
end
end
end
function T = calcT(I_vec,num_rows)
num_cols = length(I_vec);
T = zeros(num_rows,num_cols);
for i = 1:num_cols
T(I_vec(i),i) = 1;
end
end