Location: BG_cAMP @ 348913e60ceb / MATLAB / all_linalg_parameters.m

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://models.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)