- Author:
- Shelley Fong <s.fong@auckland.ac.nz>
- Date:
- 2021-05-31 18:20:34+12:00
- Desc:
- increased guess of forward kinetic rate constant in all reactions
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/674/rawfile/b22d4b553a41f8ed56c3e4a790af928d6bbe8213/MATLAB/type1_linalg_parameters.m
% This script calculates the bond graph parameters for all reactions of the
% cAMP AC-basal 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.
% coupled rxn?
clear;
clc;
close all;
%% 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
stoich_path = [data_dir filesep 'type1_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);
%% Load reverse matrix
stoich_path = [data_dir filesep 'type1_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);
%% Calculate stoichiometric matrix
nr = num_cols;
N = N_r - N_f;
N_T = N_rT - N_fT;
I = eye(nr);
M = [I N_fT; I N_rT];
% constraint of equilibrium where q_substrate = q_product
% eqn 2.36 of Pan thesis
N_cT = zeros(1,size(M,2)); % one row each per constraint
N_cT(1,3) = 1;
N_cT(1,4) = -1;
M = [M; N_cT];
%% set up loop
rx_names = {'basalAC','GsaGTPAC','FSKAC','PDE'};
%% Set the kinetic rate constants
% [kp km(kcat)] [=] 1/fM.s, 1/s
% K_M [=] fM
Km_all = [1.03e12, 3.15e11, 8.6e11, 1.3e9];
kcat_all = [0.2, 8.5, 0, 5];
kp0 = 1e-12;
for i=1:length(rx_names)
disp(rx_names(i));
Km = Km_all(i);
kcat = kcat_all(i);
% R1a/b
kap = kp0*kcat*1000; % Guess
kam = kap*Km - kcat;
kbp = kcat;
% kbm = 0; % assume reaction is irreversible
% Calculate remaining parameter using detailed balance for biochemical
% cycles e.g. single enzyme reaction
kbm = (kap*kbp)/kam;
%% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
K_C = 1;
k_kinetic = [kap kbp kam kbm]';
k = transpose([k_kinetic' K_C]);
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));
k_predicted = k_sub(3:end);
% Check that there is a detailed balance constraint
Z = transpose(null(transpose(M),'r'));
%% Save bond graph parameters
% W = [ones(18,1); W_isr; W_i; W_sr; W_i; W_i; W_i];
% lambda = lambdaW./W;
kappa = lambda(1:nr);
K = lambda(nr+1:end);
save([output_dir 'cAMP_',rx_names{i},'_params.mat'],'kappa','K','k_kinetic');
disp('kappa ');
disp(kappa);
disp('K');
disp(K);
disp(newline)
end