cba_check_model

 CBA_CHECK_MODEL - Check kinetic model for being in an enzyme-optimal state

 [isFeasible,res] = cba_check_model(network,c,u,v,cba_constraints)

 Check a kinetic model for a feasible stationary state

 'res' contains various kinds of information about the state: 

   res.c_adjusted
   res.v_adjusted
   res.u_adjusted
   res.c_mismatch_norm
   res.v_mismatch_norm
   res.mismatch_v_constraints
   res.stationarity_mismatch_norm
   res.flux_sign_mistakes
   res.dynamically_instable
   res.negative_control
   res.R_f_u
   res.economically_unstable
   res.fitness_gradient
   res.fitness_hessian
   res.benefit_hessian
   res.cost_hessian

0001 function [isFeasible, res] = cba_check_model(network, c, u, v, cba_constraints)
0002 
0003 % CBA_CHECK_MODEL - Check kinetic model for being in an enzyme-optimal state
0004 %
0005 % [isFeasible,res] = cba_check_model(network,c,u,v,cba_constraints)
0006 %
0007 % Check a kinetic model for a feasible stationary state
0008 %
0009 % 'res' contains various kinds of information about the state:
0010 %
0011 %   res.c_adjusted
0012 %   res.v_adjusted
0013 %   res.u_adjusted
0014 %   res.c_mismatch_norm
0015 %   res.v_mismatch_norm
0016 %   res.mismatch_v_constraints
0017 %   res.stationarity_mismatch_norm
0018 %   res.flux_sign_mistakes
0019 %   res.dynamically_instable
0020 %   res.negative_control
0021 %   res.R_f_u
0022 %   res.economically_unstable
0023 %   res.fitness_gradient
0024 %   res.fitness_hessian
0025 %   res.benefit_hessian
0026 %   res.cost_hessian
0027 
0028 
0029 % stationarity
0030 
0031 % set small enzyme levels to zero
0032 
0033 u(u/median(u) < 0.001) = 0;
0034 
0035 [nm,nr] = size(network.N);
0036 
0037 network.kinetics.u = u;
0038 
0039 [c_st,v_st]    = network_steady_state(network,c);
0040 
0041 res.c_adjusted = c_st;
0042 res.v_adjusted = v_st;
0043 res.u_adjusted = u;
0044 
0045 res.c_mismatch_norm = norm(c_st - c)/norm(c);
0046 res.v_mismatch_norm = norm(v_st - v)/norm(v);
0047 
0048 res.mismatch_v_constraints     = sum(v<cba_constraints.v_min) + sum(v>cba_constraints.v_max);
0049 res.stationarity_mismatch_norm = norm(network.N(network.external==0,:) * v_st);
0050 res.flux_sign_mistakes         = find([v_st~=0] .*  [sign(v_st .* [ log(network.kinetics.Keq) - network.N'*log(c_st)]) ==-1]);
0051 
0052 R   = basic_control_analysis(network,c_st);
0053 RJu = R.RJ(:,1:nr);
0054 
0055 res.dynamically_instable =  sum(real(eig(R.M))>0)>0;
0056 
0057 R_f_u =  cba_constraints.zv' * RJu;
0058 R_f_u(abs(R_f_u)<0.001*median(abs(R_f_u))) = 0;
0059 
0060 res.negative_control = find([u~=0].*R_f_u' <0);
0061 
0062 if res.negative_control,
0063   res.R_f_u = R_f_u;
0064 end
0065 
0066 if length(res.negative_control),
0067   figure(10); netgraph_concentrations(network,[], [u~=0].*R_f_u',1); title('Fitness/enzyme response coefficients');
0068 end
0069 
0070 [f_cost, grad_cost, H_cost] = feval(cba_constraints.fitness_options.cost_function, u, cba_constraints.fitness_options);
0071 
0072 fitness_gradient = [cba_constraints.zv' * RJu] - grad_cost';
0073 res.fitness_gradient_active_norm = norm(fitness_gradient .* [u ~=0]');
0074 
0075 [Ec,Ep,parameters,Ecc,Ecp,Epp,p] = elasticities(network,c);
0076 [CJ, CS] = control_coefficients(network.N, Ec, network.external);
0077 [RS,RJ,RS2,RJ2] = response_coefficients(CS,Ec,Ep,Ecc,Ecp,Epp);
0078 
0079 H_benefit  = squeeze(tensor_product(cba_constraints.zv', RJ2(:,1:nr,1:nr)));
0080 H_fitness  = H_benefit - H_cost;
0081 
0082 res.economically_unstable = sum(eig(H_fitness(find(u>0),find(u>0)))>0)>0;
0083 
0084 res.fitness_gradient = fitness_gradient;
0085 res.fitness_hessian  = H_fitness;
0086 res.benefit_hessian  = H_benefit;
0087 res.cost_hessian     = H_cost;
0088 
0089 isFeasible = ...
0090     [res.c_mismatch_norm^2/nm <10^-5] ...
0091     * [res.v_mismatch_norm^2/nr <10^-5] ...
0092     * [res.stationarity_mismatch_norm^2/nr <10^-5] ...
0093     * [res.mismatch_v_constraints == 0] ...
0094     * [length(res.flux_sign_mistakes) == 0] ...
0095     * [res.dynamically_instable == 0] ...
0096     * [length(res.negative_control) == 0] ...
0097     * [res.fitness_gradient_active_norm^2/sum(u>0) < 10^-5] ...
0098     * [res.economically_unstable == 0];