- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-06-21 16:30:43+12:00
- Desc:
- Codes to reproduce Saucerman Figs 2 and 3a
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/674/rawfile/348913e60cebcf090b62de195a8a13165b1512f7/MATLAB/all_linalg_parameters.m
% This script calculates the bond graph parameters for all reactions of the
% cAMP pathway, based on SERCA model of Pan et al, which is based
% on Tran et al. (2009). Parameters calculated by using the kinetic
% parameters and stoichiometric matrix.
clear;
clc;
close all;
%% booleans
include_FSK = false;
include_type2_reactions = true;
use_type2_constraints = false;
use_type2_constraints = use_type2_constraints * include_type2_reactions;
%% Set directories
current_dir = cd;
Idx_backslash = find(current_dir == filesep);
data_dir = [current_dir filesep 'data' filesep];
output_dir = [current_dir filesep 'output' filesep];
%% Load forward matrix
if include_type2_reactions
stoich_path = [data_dir filesep 'all_forward_matrix.txt'];
else
stoich_path = [data_dir filesep 'all_noType2_forward_matrix.txt'];
end
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
if include_type2_reactions
stoich_path = [data_dir filesep 'all_reverse_matrix.txt'];
else
stoich_path = [data_dir filesep 'all_noType2_reverse_matrix.txt'];
end
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
if include_FSK == false
fsk_rxn_ids = [5,6,11];
N_f(:,fsk_rxn_ids) = [];
N_r(:,fsk_rxn_ids) = [];
N_f([7,8,15],:) = [];
N_r([7,8,15],:) = [];
num_cols = num_cols - length(fsk_rxn_ids);
end
N_fT = transpose(N_f);
N_rT = transpose(N_r);
%% Calculate stoichiometric matrix
N = N_r - N_f;
N_T = N_rT - N_fT;
I = eye(num_cols);
M = [I N_fT; I N_rT];
% show substrate ATP is in eqlm with product cAMP
if use_type2_constraints
N_cT = zeros(5,size(M,2)); % one row each per constraint
else
N_cT = zeros(2,size(M,2));
end
N_cT(1,num_cols + 1) = 1;
N_cT(1,num_cols + 2) = -1;
% another constraint for equilibrium between cAMP and 5AMP
N_cT(2,num_cols + 2) = 1;
N_cT(2,num_cols + 11) = -1;
if use_type2_constraints
% equilibrium between PDEinh and PDE/IBMX
N_cT(3,num_cols + 9) = 1;
N_cT(3,num_cols + 12) = 1;
N_cT(3,num_cols + 13) = -1;
% equilibrium between GsAC and Gs/AC
N_cT(4,num_cols + 3) = 1;
N_cT(4,num_cols + 14) = 1;
N_cT(4,num_cols + 5) = -1;
% equilibrium between FSKAC and FSK/AC
N_cT(5,end) = 1;
N_cT(5,num_cols + 3) = 1;
N_cT(5,num_cols + 7) = -1;
end
M = [M; N_cT];
%% Set the kinetic rate constants
% knowns: K_m [=] fmol/L
% kcat [=] per_sec
K1_m = 1.03e12;
k1cat = 0.2;
K2_m = 3.15e+11;
k2cat = 8.5;
K3_m = 8.60e+11;
k3cat = 1e-16;
K4_m = 1.30e+09;
k4cat = 5;
K5_m = 3.00e+10;
K6_m = 4.00e+08;
K7_m = 4.40e+10;
if false
% use k1m == k2p (GUESS. Not a new equilibrium constraint)
% k1p = 2kcat/Km
% use detailed balance to find remaining k2m
% parameter
% for j = 1:4
% end
const = 1e-50;
t1mult = 1e0;
t2mult = 1e8;
% R1 a/b
k1ap = 2*k1cat/K1_m/t1mult;
k1am =const; %k1cat;
k1bp =const; %k1cat;
% R2 a/b
k2ap = 2*k2cat/K2_m/t1mult;
k2am =const; %k2cat;
k2bp =const; %k2cat;
% R3 a/b
k3ap = 2*k3cat/K3_m/t1mult;
k3am =const; %k3cat;
k3bp =const; %k3cat;
% R4 a/b
k4ap = 2*k4cat/K4_m/t1mult;
k4am =const; %k4cat;
k4bp =const; %k4cat;
if include_type2_reactions
% R5
k5p =const; %K5_m/t2mult;
k5m =const; %K5_m/t2mult;
% R6
k6p =const; %K6_m/t2mult;
k6m =const; %K6_m/t2mult;
% R7
k7p =const; %K7_m/t2mult;
k7m =const; %K7_m/t2mult;
end
else
kappamult = 1e-12;
% R1 a/b
k1ap = 0.93*kappamult; % 59 dessauer_1997 (no kappamult)
k1am = k1ap*K1_m - k1cat;
k1bp = k1cat;
% R2 a/b
k2ap = 0.1*kappamult; % 1 iancu
k2am = k2ap*K2_m - k2cat;
k2bp = k2cat;
% R3 a/b
k3ap = 1; % Guess - like Gsa
k3am = k3ap*K3_m - k3cat;
k3bp = k3cat;
% R4 a/b
k4ap = 1e0; % Guess 1e6
k4am = k4ap*K4_m - k4cat;
k4bp = k4cat;
if include_type2_reactions
const = 1;
% R5
k5p =100; %K5_m/t2mult;
k5m =100; %K5_m/t2mult;
% R6
k6p =const; %K6_m/t2mult;
k6m =const; %K6_m/t2mult;
% R7
k7p =const; %K7_m/t2mult;
k7m =const; %K7_m/t2mult;
% k5p = 1;
% k5m = K5_m*k5p;
% % R6
% k6p = 1;
% k6m = K6_m*k6p;
% % R7
% k7p = 1;
% k7m = K7_m*k7p;
end
end
% Calculate remaining parameter using detailed balance
k1bm = (k1ap*k1bp)/k1am;
k2bm = (k2ap*k2bp)/k2am;
k3bm = (k3ap*k3bp)/k3am;
k4bm = (k4ap*k4bp)/k4am;
if include_type2_reactions == false
k5p = [];
k6p = [];
k7p = [];
k5m = [];
k6m = [];
k7m = [];
end
if include_FSK == false
k3ap = [];
k3bp = [];
k3am = [];
k3bm = [];
k7p = [];
k7m = [];
end
%% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
Kc = ones(1,size(N_cT,1));
% Kc = [];
k_kinetic = [k1ap k1bp ...
k2ap k2bp ...
k3ap k3bp ...
k4ap k4bp ...
k5p ...
k6p ...
k7p ...
k1am k1bm ...
k2am k2bm ...
k3am k3bm ...
k4am k4bm ...
k5m ...
k6m ...
k7m]';
k = transpose([k_kinetic' Kc]);
lambda = exp(pinv(M) * log(k));
% Check that kinetic parameters are reproduced by bond graph parameters
k_sub = exp(M*log(lambda));
diff = (k - k_sub)./k;
error = sum(abs(diff));
% Check that there is a detailed balance constraint
Z = transpose(null(transpose(M),'r'));
%% Save bond graph parameters
kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);
if include_type2_reactions == false
kappa = [kappa;0;0;0];
K = [K;100;100;100;100];
end
save([output_dir 'all_params.mat'],'kappa','K','k_kinetic');
if include_FSK
rxnID = {'1a','1b','2a','2b','3a','3b','4a','4b','5','6','7'};
Kname = {'ATP','cAMP','AC','AC_ATP','Gs_AC','Gs_AC_ATP','FSK_AC','FSK_AC_ATP','PDE', ...
'PDE_cAMP','5AMP','IBMX','PDEinh','Gs','FSK'};
else
rxnID = {'1a','1b','2a','2b','4a','4b','5','6'};
Kname = {'ATP','cAMP','AC','AC_ATP','Gs_AC','Gs_AC_ATP','PDE', ...
'PDE_cAMP','5AMP','IBMX','PDEinh','Gs'};
end
for j = 1:length(kappa)
fprintf('var kappa_%s: fmol_per_sec {init: %g};\n',rxnID{j},kappa(j));
end
for j = 1:length(K)
fprintf('var K_%s: per_fmol {init: %g};\n',Kname{j},K(j));
end
if include_FSK == false
fprintf('var kappa_3a: fmol_per_sec {init: 0};\n');
fprintf('var kappa_3b: fmol_per_sec {init: 0};\n');
fprintf('var kappa_7: fmol_per_sec {init: 0};\n');
fprintf('var K_FSK: per_fmol {init: 1e3};\n');
fprintf('var K_FSK_AC: per_fmol {init: 1e3};\n');
fprintf('var K_FSK_AC_ATP: per_fmol {init: 1e3};\n');
end
disp(newline)
