spectral_response_coefficients.m

来自「％The Metabolic Networks Toolbox contains」· M 代码 · 共 54 行

%[RS,RJ,RS_2_2omega,RJ_2_2omega,RS_2_2omega0,RV_2_2omega0,CS,CS_2_omega,CJ,CJ_2_omega] = ...
%  spectral_response_coefficients(N0, L, M0, delta_par,epsilon_1,pi_1,epsilon_2,rho_2,pi_2)
%
%calculate frequency-dependent response coefficients
%
% L                   link matrix
% N0                  reduced stoichiometric matrix
% delta_par           structure with field 'omega': angular frequency 
% epsilon, pi, rho... elasticities of 1st and 2nd order
% M0                  Jacobian matrix (optional)

function [RS_omega,RJ_omega,RS_2_2omega,RJ_2_2omega,RS_2_2omega0,RV_2_2omega0,CS_omega,CS_2_omega,CJ_omega,CJ_2_omega] = spectral_response_coefficients(N0,L,M0,delta_par,epsilon_1,pi_1,epsilon_2,rho_2,pi_2)

n_react = size(N0,2);
n_indep = size(N0,1);
clear i;

if ~exist('M0','var'), 
  M0 = N0 * epsilon_1 * L; 
end

CS          = - L * inv(M0                                    ) * N0;
CS_omega    = - L * inv(M0 -i *   delta_par.omega*eye(n_indep)) * N0;

CJ          =  epsilon_1 * CS         + eye(n_react);
CJ_omega    =  epsilon_1 * CS_omega   + eye(n_react);
 
RS_omega    = CS_omega * pi_1;
RJ_omega    = CJ_omega * pi_1;

if nargout>2,

CS_2_omega  = - L * inv(M0 -i * 2*delta_par.omega*eye(n_indep)) * N0;
CJ_2_omega  =  epsilon_1 * CS_2_omega + eye(n_react);

Gamma_omega_omega =  ...
    tensor_product(tensor_product(epsilon_2,RS_omega),RS_omega,2,1) ...
    + tensor_product(rho_2,RS_omega,2,1) ...
    + permute(tensor_product(rho_2,RS_omega,2,1),[1,3,2])...
    + pi_2;

Gamma_omega_minusomega =  ...
    tensor_product(tensor_product(epsilon_2,RS_omega.''),RS_omega,2,1) ...
    + tensor_product(rho_2,RS_omega.'',2,1)  ...
    + permute(tensor_product(rho_2,RS_omega,2,1),[1,3,2])...
    + pi_2;

RS_2_2omega  =   1/sqrt(2*pi) *  tensor_product(CS_2_omega,Gamma_omega_omega);
RJ_2_2omega  =   1/sqrt(2*pi) *  tensor_product(CJ_2_omega,Gamma_omega_omega);

RS_2_2omega0 =   1/sqrt(2*pi) *  tensor_product(CS,Gamma_omega_minusomega);
RV_2_2omega0 =   1/sqrt(2*pi) *  tensor_product(CJ,Gamma_omega_minusomega);

end