Location: cellLib @ a7940fcacd45 / BG_fit / BG_paras_conv.m

Author:

WeiweiAi <wai484@aucklanduni.ac.nz>

Date:

2023-04-04 15:51:18+12:00

Desc:

Add a steady state example

Permanent Source URI:

https://models.cellml.org/workspace/6bc/rawfile/a7940fcacd455a09bdb7beac0da91eca79cf1adb/BG_fit/BG_paras_conv.m

function [kappa,K,diff]=BG_paras_conv(N_f,N_r,kf,kr,K_c,N_c,Ws)
% N_f/N_r: forward and reverse stoichiometric matrices
% kf/kr: column vectors of the forward and reverse kinetic
%rate constants; K_c: column vector of known constraints between the species
%defined in N_c; Ws: column vector of the volume of species
%% Construct M
N_fT = transpose(N_f);
N_rT = transpose(N_r);
N = N_r - N_f;
num_cols = size(N,2);
I = eye(num_cols);
zerofill = zeros(1,num_cols);
N_c_T=transpose(N_c);
if ~isempty(N_c)
M = [I N_fT; I N_rT;zerofill N_c_T;];
else
    M = [I N_fT; I N_rT];
end
%% Construct W
W = [ones(size(N,2),1);Ws;];
%% Contruct vector of kinetic paras
if ~isempty(N_c)
k = [kf;kr;K_c];
else
    k = [kf;kr];
end
%% Convert kinetic paras to BG paras
lambdaW = exp(pinv(M)*log(k));
lambda = lambdaW./W;
kappa = lambda(1:size(N,2));
K = lambda(size(N,2)+1:end);
%% Checks
N_rref = rref(N);
R_mat = null(N,'r');
K_eq = kf./kr;
zero_est = R_mat'*K_eq;
k_est = exp(M*log(lambdaW));
diff = sum(abs((k-k_est)./k));