Location: BG_cAMP @ b22d4b553a41 / MATLAB / cAMP_basalAC_parameters.m

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/cAMP_basalAC_parameters.m

% This script calculates the bond graph parameters for 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.

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 'basalAC_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 'basalAC_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
N = N_r - N_f;
N_T = N_rT - N_fT;

I = eye(num_cols);

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,2),1,-1,0,0]; 

M = [M; N_cT];

%% Set the kinetic rate constants
K_m = 1.03e12; % [=] fmol/L
kcat = 0.2; % [=] per_sec
% R1
k1p = 1e-12; % Guess
k1m = k1p*K_m - kcat;

% R2
k2p = kcat;
% k2m = 0; % assume reaction is  irreversible

% Calculate remaining parameter using detailed balance
k2m = (k1p*k2p)/k1m;

%% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
% k = transpose([k_12p k_24p k_22ap k_45p k_56p k_68p k_89p k_910p k_101p k_12m k_24m k_22am k_45m k_56m k_68m k_89m k_910m k_101m K_MgATP 1]);
k = transpose([k1p k2p k1m k2m 1]);
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
% W = [ones(18,1); W_isr; W_i; W_sr; W_i; W_i; W_i];
% lambda = lambdaW./W;

kappa = lambda(1:num_cols);
K = lambda(num_cols+1:end);

save([output_dir 'cAMP_basalAC_params.mat'],'kappa','K');