- Author:
- WeiweiAi <wai484@aucklanduni.ac.nz>
- Date:
- 2022-08-29 15:22:53+12:00
- Desc:
- Remove q and its ode in Se_BG
- Permanent Source URI:
- https://models.physiomeproject.org/workspace/6bc/rawfile/be4db1cb7bdb17d9500fe90012ba8759b0ee41e9/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));
