Location: BG_cAMP @ a51cfe492772 / MATLAB / cAMP_PDE_IBMX_parametersm.m

Author:

Shelley Fong <s.fong@auckland.ac.nz>

Date:

2021-06-04 16:32:29+12:00

Desc:

New BG parameters

Permanent Source URI:

https://models.physiomeproject.org/workspace/674/rawfile/a51cfe49277237dc75e18131fa08ca3089c799c1/MATLAB/cAMP_PDE_IBMX_parametersm.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.

% UNFINISHED: WILL NOT WORK WITH MATRIX ROUTINE

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 'IBMX_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 'IBMX_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];

N_cT = [0,1,1,-1]; % constraint of equilibrium where q_sub = q_prod. Take Ksub1*Ksub2 = Kprod

M = [M; N_cT];

%% Set the kinetic rate constants
K_m = 30e9; % [=] fmol/L
% R1
k1p = 1e-12; % Guess
k1m = k1p*K_m;

% detailed balance not possible for this type of reaction

%% Calculate bond graph constants from kinetic parameters
% Note: units of kappa are fmol/s, units of K are fmol^-1
k = transpose([k1p k1m 1]);
lambda = exp(pinv(M) * log(k));

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

% 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));
% Km_kinetic = Ksub1Ksub2/Kprod; % equilibrium constraint
Km_kinetic = K(2)*K(1)/K(3);
Km_diff = (K_m - Km_kinetic)/K_m;

% 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;

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