cba_feasible_lp

 CBA_FEASIBLE_LP - Check flux mode for EFA feasibility and choose enzyme costs and thermodynamic forces

 [isFeasible, y, dmu] = cba_feasible(v, network, cba_constraints, cba_options)

 Test a flux vector v for feasibility (FBA, EBA, CBA) using linear programming

 Consider network with internal stoichiometric matrix N_int:
  - check if flux mode v is enzyme-optimal 
  - if so, return a possible vector 'y' of rate costs for this state
    (y = h_u * u / v), same signs as the rates

 Structure 'cba_constraints' must contains the fields:

 either    zv   : benefits for fluxes 
 or        z_ext: benefits of network.external metabolites 

 N_tot     : entire stoichiometric matrix; required for EBA test

 cba_options.test_eba : flag: should thermodynamical feasibility be tested as well?

0001 function [isFeasible, y, dmu] = cba_feasible_lp(v, network, cba_constraints, cba_options)
0002 
0003 % CBA_FEASIBLE_LP - Check flux mode for EFA feasibility and choose enzyme costs and thermodynamic forces
0004 %
0005 % [isFeasible, y, dmu] = cba_feasible(v, network, cba_constraints, cba_options)
0006 %
0007 % Test a flux vector v for feasibility (FBA, EBA, CBA) using linear programming
0008 %
0009 % Consider network with internal stoichiometric matrix N_int:
0010 %  - check if flux mode v is enzyme-optimal
0011 %  - if so, return a possible vector 'y' of rate costs for this state
0012 %    (y = h_u * u / v), same signs as the rates
0013 %
0014 % Structure 'cba_constraints' must contains the fields:
0015 %
0016 % either    zv   : benefits for fluxes
0017 % or        z_ext: benefits of network.external metabolites
0018 %
0019 % N_tot     : entire stoichiometric matrix; required for EBA test
0020 %
0021 % cba_options.test_eba : flag: should thermodynamical feasibility be tested as well?
0022 
0023 
0024 % ------------------------------------------------------------------
0025 % if necessary, reduce to active subnetwork and call the function
0026 
0027 % set all small fluxes to zero
0028 v(abs(v)/max(abs(v))<10^-5)=0;
0029 
0030 if sum(v==0),
0031   N_int = network.N(find(network.external==0),:);
0032   Es = -network.N';
0033   [v_act,N_int_act,Es_act,network_act,cba_constraints_act,ind_act,ind_met_act] = cba_reduce_to_active_subnetwork(v,N_int,Es,network,cba_constraints);
0034   [isFeasible, y_act, dmu_act] = cba_feasible_lp(v_act,network_act,cba_constraints_act,cba_options);
0035   [nm,nr]             = size(network.N);
0036   y           = nan *ones(nr,1);
0037   y(ind_act)  = y_act;  
0038   dmu            = nan *ones(nr,1);
0039   if isFeasible,
0040     dmu(ind_act)   = dmu_act;
0041   end
0042   return
0043 end
0044 
0045 
0046 % ------------------------------------------------------------------
0047 %% initialisation
0048 
0049 cba_constraints = cba_update_constraints(cba_constraints,network.N(find(network.external),:));
0050 
0051 ind_int = find(network.external==0);
0052 N_int   = network.N(ind_int,:);
0053 isFeasible  = 1;
0054 y         = [];
0055 dmu  = [];
0056 zv        = cba_constraints.zv;
0057   
0058 if ~isfield(cba_options,'test_eba'), cba_options.test_eba = 1; end
0059 
0060 %% test if mode is stationary
0061   
0062 epsilon_balance = 10^-4;
0063   
0064 if max(abs(N_int*v))>epsilon_balance,
0065   isFeasible = 0;  display('Flux mode is not stationary.');
0066 end
0067 
0068 % ------------------------------------------------------------------
0069 % first condition: mode yields positive benefit?
0070 
0071 v_w     = zv'*v;
0072   
0073 switch sign(v_w),
0074   case  0,  isFeasible = 0; display('No benefit produced.');
0075   case -1,  isFeasible = 0; display('Negative  benefit produced.');
0076 end
0077 
0078 % ------------------------------------------------------------------
0079 % second condition: existence of enzyme cost vector y
0080 
0081 if isFeasible,  
0082     
0083   %% reduce model to active reactions and remove all
0084   %% irrelevant internal metabolites
0085   
0086   nn = network_construct(N_int);
0087   cc = cba_constraints;
0088   cc.z_ext = [];
0089   
0090   %% test if sign pattern is feasible: search for a vector y such that
0091   %%
0092   %%      K' * y = K' * zv
0093   %% diag(v) * y > 0
0094   %%
0095   %% solve this by linear programming:
0096     
0097   K  = analyse_N(N_int);    
0098   nr = size(N_int,2);
0099     
0100   epsilon = 10^-4; % threshold enzyme cost
0101   c       = ones(nr,1);
0102   h       = epsilon * ones(nr,1);
0103   G       = diag(v);
0104   A       = K';
0105   b       = K' * cc.zv;
0106 
0107   %%y_act = lp236a(-c,-G,-h,A,b);
0108 
0109   if exist('cplexlp','file'),
0110     opt = optimset('Display','off');
0111     [y, fval, exitflag] = linprog(-c, -G, -h, A, b,[],[], [], opt);
0112   else
0113     opt = cplexoptimset('Display','off');
0114     [y, fval, exitflag] = cplexlp(-c, -G, -h, A, b,[],[], [], opt);
0115   end
0116 
0117   if exitflag ~= 1,
0118     exitflag
0119     isFeasible = 0;  display('Flux distribution seems to be economically infeasible.');
0120   end
0121   
0122   if cba_options.test_eba,
0123     [eba_feas,dmu] = eba_feasible_lp(v,cba_constraints.N_tot);
0124     isFeasible = isFeasible * eba_feas;
0125   end
0126 
0127 end