Location: BG_cAMP @ 348913e60ceb / MATLAB / type2_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://staging.physiomeproject.org/workspace/674/rawfile/348913e60cebcf090b62de195a8a13165b1512f7/MATLAB/type2_linalg_parameters.m

% This script calculates the bond graph parameters for TYPE 2 REACTIONS which are simple a+b >< c
% Solve for kappa KA KB KC using linalg methods
% require kp km kinetic parameters

clear;
% clc;
% close all;

%% Calculate stoichiometric matrix
% M = [I Nf; I Nr; 0 NcT]
M = [[1,1,1,0];[1,0,0,1];[0,1,1,-1]];

%% Set up loop
rx_names = {'PDE-IBMX','GsaGTPAC-AC','FSK-AC'};

%% Set the kinetic rate constants
% [kp km] [=] 1/fM.s, 1/s
K_M = [30e9, 0.4e9, 44e9]; %[=] fM

kp0 = 1e6;
k_all = [kp0,K_M(1)*kp0;
    kp0,K_M(2)*kp0;
    kp0,K_M(3)*kp0];

for i = 1:length(rx_names)
    disp(rx_names(i));
    %% Calculate bond graph constants from kinetic parameters
    % Note: units of kappa are fmol/s, units of K are fmol^-1
    kk = k_all(i,:);
    Kc = 1;
    k = transpose([kk 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'));

    % check that equlibrium K_A*K_B/K_C = 1
    Keq = lambda(1)*lambda(2)/lambda(3);
    disp('equilibrium Keq is');
    disp(Keq)
    %% Save bond graph parameters
    % W = [ones(18,1); W_isr; W_i; W_sr; W_i; W_i; W_i];
    % lambda = lambdaW./W;

    disp('kappa ');
    disp(lambda(1));
    disp('K');
    disp(lambda(2:end));

end