CBA_ECONOMIC_POTENTIALS - Compute the economic potentials and their sensitivities in a kinetic model [w, v, dw_du, dv_du] = cba_economic_potentials(network, c, cba_constraints,ind_controllable) w contains only the internal economic potentials Compute the economic potentials for a kinetic model in a given stationary state as well as their derivatives dw_du with respect to the controllable enzyme levels The calculation works only for linear benefit functions (defined by weights zv, zc)
0001 function [w, v, dw_du, dv_du, c] = cba_economic_potentials(network, c, zv, zc, ind_controllable) 0002 0003 % CBA_ECONOMIC_POTENTIALS - Compute the economic potentials and their sensitivities in a kinetic model 0004 % 0005 % [w, v, dw_du, dv_du] = cba_economic_potentials(network, c, cba_constraints,ind_controllable) 0006 % 0007 % w contains only the internal economic potentials 0008 % Compute the economic potentials for a kinetic model in a given stationary state 0009 % as well as their derivatives dw_du with respect to the controllable enzyme levels 0010 % 0011 % The calculation works only for linear benefit functions (defined by weights zv, zc) 0012 0013 0014 % --------------------------------------------------------- 0015 % to be sure ... put system into steady state 0016 0017 [c,v] = network_steady_state(network,c,200); 0018 [nm,nr] = size(network.N); 0019 ind_int = find(network.external==0); 0020 0021 0022 % --------------------------------------------------------- 0023 % Compute economic potentials w 0024 % phi: virtual supply fluxes for independent metabolites 0025 0026 0027 Ec = elasticities(network,c); 0028 [CJ, CS, L_int, NR_int, M, Madj, indp_among_internal] = control_coefficients(network.N, Ec, network.external); 0029 0030 CS_int_phi = - L_int * pinv(M); 0031 CS_phi = zeros(nm,length(indp_among_internal)); 0032 CS_phi(ind_int,:) = CS_int_phi; 0033 CJ_phi = Ec(:,ind_int) * CS_int_phi; 0034 w_ind = [zv' * CJ_phi + zc' * CS_phi]'; 0035 w = zeros(nm,1); 0036 w(ind_int(indp_among_internal)) = w_ind; 0037 0038 0039 % --------------------------------------------------------- 0040 % Compute sensitivities dw_du and dv_du 0041 0042 if nargout > 2, 0043 0044 %% NUMERICAL DERIVATIVES 0045 uu = network.kinetics.u; 0046 dd = 10^-8; 0047 0048 for itt = 1:nr, 0049 my_du1 = 0 * uu; 0050 my_du1(itt) = uu(itt) * dd; 0051 my_du2 = -my_du1; 0052 my_du1(uu==0) = dd; 0053 my_du2(uu==0) = 0; 0054 network.kinetics.u = uu + my_du1; 0055 [my_w1, my_v1] = cba_economic_potentials(network, c, zv, zc, ind_controllable); 0056 network.kinetics.u = uu + my_du2; 0057 [my_w2, my_v2] = cba_economic_potentials(network, c, zv, zc, ind_controllable); 0058 dw_du(:,itt) = [my_w1 - my_w2] / [my_du1(itt) - my_du2(itt)]; 0059 dv_du(:,itt) = [my_v1 - my_v2] / [my_du1(itt) - my_du2(itt)]; 0060 end 0061 0062 %% ANALYTICAL CALCULATION: 0063 %% R is not actually control matrix, but "response" to enzyme levels 0064 % Ec = elasticities(network,c); 0065 % [Ec,Ep,parameters,Ecc,Ecp,Epp,p] = elasticities(network,c); 0066 % Ecu = Ecp(:,:,1:nr); 0067 % v = network_velocities(c,network); 0068 % rho_2 = diag(1./v) * Ec; 0069 %% VORSICHT: GAMMA-FORMEL IST eventuell noch FALSCH!!!! 0070 %% es wird keine term ber�ksichtigt, der den direkten einfluss 0071 %von phi mit dem indirekten von u koppelt. nochmal nachrechnen! 0072 % Gamma = tensor_product(tensor_product(Ecc,CS_phi),CS,2,1) ... 0073 % + permute(tensor_product(Ecu,CS_phi,2,1),[1,3,2]); 0074 % RS_phi_u = tensor_product(CS,Gamma); 0075 % RJ_phi_u = tensor_product(Ec,RS_phi_u); 0076 % 0077 % dw_du_ind = squeeze(tensor_product(zv', RJ_phi_u) + tensor_product(zc', RS_phi_u) ); 0078 % dw_du = zeros(nm,nr); 0079 % dw_du(ind_int(indp_among_internal),:) = dw_du_ind; 0080 0081 end