Location: BG_cAMP @ d29401fb021b / MATLAB / type1_linalg_parameters.m

Author:

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

Date:

2021-07-22 12:36:57+12:00

Desc:

Removing unnecessary files

Permanent Source URI:

https://staging.physiomeproject.org/workspace/674/rawfile/d29401fb021b702a4deac586d4d7ec0eb7717f24/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