cba_reconstruct_model_OLD

 CBA_RECONSTRUCT_MODEL - Build model from economical flux profile

 %% OLD VERSION IN WHICH THE equality pc = G' pcm - zc is not yet used

 [network, res, cba_constraints] = cba_reconstruct_model(network,v,mu,cba_constraints,cba_options,y,w,c)

 The function determines all model quantitaties and the rate laws 
 based on an economical flux mode. 

 Important: it is assumed that all enzymes in the model are controllable


 Input
   network  - Network structure (as in Metabolic Network Toolbox)
   v        - Precalculated metabolic fluxes
   w        - Precalculated economic potentials
   y        - Precalculated enzyme costs
   mu       - Precalculated chemical potentials
   c        - Precalculated metabolite levels
              
   For the inputs cba_constraints and cba_options, see cba_default_options
   y and w together must satisfy the economic balance equation

 Output
   network          - Model with rate laws (in field 'kinetics');
   res              - All results in matlab struct
   cba_constraints  - Updated constraints data structure

0001 function [network, res, cba_constraints] = cba_reconstruct_model(network,v,mu,cba_constraints,cba_options,y,w,c)
0002 
0003 % CBA_RECONSTRUCT_MODEL - Build model from economical flux profile
0004 %
0005 % %% OLD VERSION IN WHICH THE equality pc = G' pcm - zc is not yet used
0006 %
0007 % [network, res, cba_constraints] = cba_reconstruct_model(network,v,mu,cba_constraints,cba_options,y,w,c)
0008 %
0009 % The function determines all model quantitaties and the rate laws
0010 % based on an economical flux mode.
0011 %
0012 % Important: it is assumed that all enzymes in the model are controllable
0013 %
0014 %
0015 % Input
0016 %   network  - Network structure (as in Metabolic Network Toolbox)
0017 %   v        - Precalculated metabolic fluxes
0018 %   w        - Precalculated economic potentials
0019 %   y        - Precalculated enzyme costs
0020 %   mu       - Precalculated chemical potentials
0021 %   c        - Precalculated metabolite levels
0022 %
0023 %   For the inputs cba_constraints and cba_options, see cba_default_options
0024 %   y and w together must satisfy the economic balance equation
0025 %
0026 % Output
0027 %   network          - Model with rate laws (in field 'kinetics');
0028 %   res              - All results in matlab struct
0029 %   cba_constraints  - Updated constraints data structure
0030 
0031 
0032 % ----------------------------------------------------------------------
0033 % Restrict model to active submodel (variable names ..._act) and initialise variables
0034 
0035 % Q_ext is upposed to have length = #external metab  (or 0, if omitted)
0036 
0037 if isempty(cba_constraints.Q_ext), 
0038   Q_ext = nan * ones(sum(network.external),1); 
0039 else, 
0040   Q_ext = cba_constraints.Q_ext;
0041 end 
0042 
0043 if isfield(network,'kinetics'), network = rmfield(network,'kinetics'); end
0044 
0045 [nm,nr] = size(network.N);
0046 ind_int = find(network.external==0);
0047 ind_ext = find(network.external==1);
0048 
0049 [v_act, N_int_act, Es_act, network_act, cba_constraints_act, ind_active, ind_met_active] = cba_reduce_to_active_subnetwork(v,network.N(ind_int,:),[],network,cba_constraints);
0050 
0051 ind_ext_act = find(network_act.external);
0052 ind_int_act = find(network_act.external==0);
0053 mu_act      = mu(ind_met_active);
0054 y_act       = y(ind_active);
0055 w_act       = w(ind_met_active);
0056 c_act       = c(ind_met_active);
0057 N_act       = network_act.N;
0058 W_act       = network_act.regulation_matrix;
0059 
0060 % assume hill coefficients = 1 (already implicitly assumed in 'make_structure_matrices')
0061 [Mplus_act, Mminus_act, Wplus_act, Wminus_act, nm_act, nr_act] = make_structure_matrices(N_act,W_act,ind_ext_act);
0062 M_act       = Mplus_act + Mminus_act;
0063 h_act       = ones(nr_act,1); 
0064 n_ext_act   = sum(network_act.external);
0065 A_act       = - N_act'* mu_act;
0066 zeta_act    = exp(h_act .* A_act / RT);
0067 hu_act      = cba_constraints.hu(ind_active);
0068 
0069 % note that here, the vector zc refers to all (not only internal)
0070 % metabolites, but the entries for external metabolites are ignored
0071 
0072 if length(cba_constraints_act.zc),
0073   zc_act = cba_constraints_act.zc(ind_met_active);
0074 else,
0075   zc_act = zeros(size(c_act));
0076 end
0077 
0078 Q_act_predefined = - c_act .* zc_act;
0079 dum              = nan * ones(nm,1);
0080 dum(find(network.external)) = Q_ext;
0081 Q_act_predefined(ind_ext_act) = dum(ind_met_active(ind_ext_act)); 
0082 
0083 
0084 % ------------------------------------
0085 % find blocks in link matrix
0086 
0087 Nint_act = N_act(find(network_act.external ==0),:); 
0088 [L_act, NR_act, ind_ind_met_act] = reduce_N(Nint_act);
0089 L_act_scaled = diag(1./c_act(ind_int_act)) * L_act * diag(c_act(ind_int_act(ind_ind_met_act)));
0090 L_blocks = matrix_find_blocks(L_act);
0091 ind_dep_met_act = setdiff(1:length(ind_int_act),ind_ind_met_act);
0092 
0093 [ni,nj] = size(L_act); 
0094 if ni==nj, 
0095   display(sprintf('\n  The model does not contain conserved moieties')); 
0096 else
0097   display('  The model contains conserved moieties'); 
0098 end
0099 
0100 
0101 % -------------------------------------------------------------
0102 % Investment balance equations; precompute necessary information for each metabolite:
0103 %   my_Q
0104 %   my_beta_M
0105 %   my_alpha_A
0106 %   my_beta_I
0107 
0108 
0109 for it = 1:nm_act,
0110   %% reactions affected by the metabolite
0111   ind_rea = find(abs(M_act(:,it))+abs(Wplus_act(:,it))+abs(Wminus_act(:,it)));
0112   my_M      =       M_act(ind_rea,it);
0113   my_Mplus  =  Mplus_act(ind_rea,it);
0114   my_Mminus = Mminus_act(ind_rea,it);
0115   my_Wplus  =  Wplus_act(ind_rea,it);
0116   my_Wminus = Wminus_act(ind_rea,it);
0117   my_y      =      y_act(ind_rea);
0118   
0119   % sign constraint for Q // set sign for external metabolites
0120   Q_min = -inf;  
0121   Q_max =  inf;  
0122   
0123   %% reactions must not be completely switched off or saturated: 0.05 < beta < 0.95
0124   %% CHECK THE SAME FOR ALLOSTERIC REGLUATION:: FIX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
0125   %% ALSO FURTHER BELOW!!!
0126 
0127   x_min = [ Q_min; 0.05 * ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0128   x_max = [ Q_max; 0.95 * ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0129   
0130   if network_act.external(it),
0131     %% external metabolites
0132     if isfinite(Q_act_predefined(it)),
0133       my_Q_ext     = Q_act_predefined(it);
0134       my_Q_ext_std = 0.1*abs(my_Q_ext);
0135     else
0136       net_prod  = N_act(it,:) * v_act;
0137       my_Q_ext = 0;
0138       if net_prod < 0, my_Q_ext =  1; end
0139       if net_prod > 0, my_Q_ext = -1; end
0140       my_Q_ext_std = 1; 
0141     end
0142     x_prior_mean = [my_Q_ext;     1/2 * ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0143     x_prior_std  = [my_Q_ext_std;       ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0144   else
0145     %% internal metabolites
0146     my_Q_given = Q_act_predefined(it);
0147     x_prior_mean = [my_Q_given;           1/2 * ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0148     x_prior_std  = [0.00001 * [1 + my_Q_given]; ones(length(find([my_M; my_Wplus; my_Wminus])),1)];
0149   end
0150   
0151   % thermodynamic part of elasticities
0152   my_E_T = [zeta_act(ind_rea) .* my_Mplus - my_Mminus] ./ [zeta_act(ind_rea)-1];
0153 
0154   my_ind_M      = find(my_M);       my_n_M      = length( my_ind_M);
0155   my_ind_Wplus  = find(my_Wplus);   my_n_Wplus  = length( my_ind_Wplus);
0156   my_ind_Wminus = find(my_Wminus);  my_n_Wminus = length( my_ind_Wminus);
0157   
0158   Aeq = full([1, ...
0159               my_y(my_ind_M)'      .*  my_M(my_ind_M)',   ...
0160             - my_y(my_ind_Wplus)'  .*  my_Wplus(my_ind_Wplus)',   ...
0161               my_y(my_ind_Wminus)' .* my_Wminus(my_ind_Wminus)']);
0162   
0163   Beq = my_y' * my_E_T;
0164   
0165   rr{it}.Aeq          = Aeq;
0166   rr{it}.Beq          = Beq;
0167   rr{it}.x_min        = x_min;
0168   rr{it}.x_max        = x_max;
0169   rr{it}.x_prior_mean = x_prior_mean;
0170   rr{it}.x_prior_std  = x_prior_std;
0171   rr{it}.my_n_M       = my_n_M  ;
0172   rr{it}.my_n_Wplus   = my_n_Wplus  ;
0173   rr{it}.my_n_Wminus  = my_n_Wminus ;
0174   
0175 end
0176 
0177 
0178 % -------------------------------------------------------------
0179 % Initialise more
0180 
0181 Q_act       = nan * ones(nm_act,1);
0182 beta_M_act  = zeros(nr_act,nm_act); beta_M_act(find(M_act)) = nan;
0183 alpha_A_act = zeros(nr_act,nm_act); alpha_A_act(find(Wplus_act)) = nan;
0184 beta_I_act  = zeros(nr_act,nm_act);  beta_I_act(find(Wminus_act)) = nan;
0185 
0186 
0187 % -------------------------------------------------------------
0188 % Solve investment balance equation for external metabolites
0189 
0190 opt = optimset('Display','off','Algorithm','interior-point-convex');
0191 
0192 for it = 1:n_ext_act,
0193   ind   = ind_ext_act(it);
0194   [x_opt,dum,exitflag] = quadprog(diag(rr{ind}.x_prior_std.^-2), -diag(rr{ind}.x_prior_std.^-2)* rr{ind}.x_prior_mean, [],[],rr{ind}.Aeq,rr{ind}.Beq,rr{ind}.x_min,rr{ind}.x_max,[],opt);
0195   if exitflag ~=1, exitflag 
0196     error('Error in optimisation'); 
0197   end 
0198   
0199   my_Q               = x_opt(1);                     x_opt = x_opt(2:end);
0200   my_beta_M_act      = x_opt(1:rr{ind}.my_n_M);      x_opt = x_opt(rr{ind}.my_n_M+1:end);
0201   my_alpha_A_act     = x_opt(1:rr{ind}.my_n_Wplus);  x_opt = x_opt(rr{ind}.my_n_Wplus+1:end);
0202   my_beta_I_act      = x_opt(1:rr{ind}.my_n_Wminus);
0203   
0204   Q_act(ind,1)                            = my_Q; 
0205   beta_M_act(find(M_act(:,ind)),ind)      = my_beta_M_act;
0206   alpha_A_act(find(Wplus_act(:,ind)),ind) = my_alpha_A_act;
0207   beta_I_act(find(Wminus_act(:,ind)),ind) = my_beta_I_act;
0208 end
0209 
0210 
0211 % -------------------------------------------------------------
0212 % Solve investment balance equation for internal metabolites
0213 % go through blocks of link matrix and solve equations for each block
0214 % (for instance, blocks with only one entry: independend metabolites
0215 %  on which no other metabolite depends)
0216 
0217 for it = 1:length(L_blocks),
0218   ind_L_col       = L_blocks{it}.columns;
0219   ind_L_row       = L_blocks{it}.rows;
0220   my_L            = L_act(ind_L_row,ind_L_col);
0221   my_L_scaled     = L_act_scaled(ind_L_row,ind_L_col);
0222 
0223   my_x_prior_std  = [];
0224   my_x_prior_mean = [];
0225   my_Aeq          = [];
0226   my_Beq          = [];
0227   my_x_min        = [];
0228   my_x_max        = [];
0229 
0230   for itt = 1:length(ind_L_row)
0231     it_int = ind_int_act(ind_L_row(itt));
0232     my_x_min        = [my_x_min;        rr{it_int}.x_min          ];
0233     my_x_max        = [my_x_max;        rr{it_int}.x_max          ];
0234     my_x_prior_std  = [my_x_prior_std;  rr{it_int}.x_prior_std;   ];
0235     my_x_prior_mean = [my_x_prior_mean; rr{it_int}.x_prior_mean   ];
0236     my_Aeq          = matrix_add_block(my_Aeq,rr{it_int}.Aeq);
0237     my_Beq          = [my_Beq;          rr{it_int}.Beq            ];    
0238   end
0239 
0240   my_Aeq = my_L_scaled' * my_Aeq;
0241   my_Beq = my_L_scaled' * my_Beq;
0242 
0243   [x_opt,fval,exitflag] = quadprog(diag(my_x_prior_std.^-2), -diag(my_x_prior_std.^-2) * my_x_prior_mean, [],[],my_Aeq,my_Beq,my_x_min,my_x_max,[],opt);
0244   if exitflag ~=1, error('Error in optimisation'); end 
0245   
0246   for itt = 1:length(ind_L_row)
0247     it_int = ind_int_act(ind_L_row(itt));
0248     my_Q               = x_opt(1);                        x_opt = x_opt(2:end);
0249     my_beta_M_act      = x_opt(1:rr{it_int}.my_n_M);      x_opt = x_opt(rr{it_int}.my_n_M+1:end);
0250     my_alpha_A_act     = x_opt(1:rr{it_int}.my_n_Wplus);  x_opt = x_opt(rr{it_int}.my_n_Wplus+1:end);
0251     my_beta_I_act      = x_opt(1:rr{it_int}.my_n_Wminus); x_opt = x_opt(rr{it_int}.my_n_Wminus+1:end);
0252   
0253     Q_act(it_int,1)                               = my_Q; 
0254     beta_M_act(find(M_act(:,it_int)),it_int)          = my_beta_M_act;
0255     alpha_A_act(find(Wplus_act(:,it_int)),it_int) = my_alpha_A_act;
0256     beta_I_act(find(Wminus_act(:,it_int)),it_int) = my_beta_I_act;
0257   end
0258 
0259 end
0260 
0261 
0262 % --------------------------------------------------------------------------
0263 % update Q values (in active subnetwork) if necessary
0264 
0265 Q_mismatch = norm(Q_act(ind_int_act) - Q_act_predefined(ind_int_act));
0266 
0267 if Q_mismatch/length(Q_mismatch) < 10^-8,
0268   display('  Predefined internal economic loads have been realised');
0269   Q_act(ind_int_act) = Q_act_predefined(ind_int_act);
0270 else,
0271   display('  Feasible solution requires change of internal economic loads');
0272   [Q_act(ind_int_act), Q_act_predefined(ind_int_act)]
0273 end
0274 
0275 Q_act(abs(Q_act)<10^-8) = 0;
0276 
0277 
0278 % --------------------------------------
0279 % Convert saturation values to kinetic constants and build kinetics data struture for active submodel
0280 
0281 alpha_M_act = alpha_to_betagamma(beta_M_act);
0282 alpha_I_act = alpha_to_betagamma(beta_I_act);
0283 
0284 KM_act = alpha_to_k(alpha_M_act,c_act,h_act);
0285 KA_act = alpha_to_k(alpha_A_act,c_act,h_act);
0286 KI_act = alpha_to_k(alpha_I_act,c_act,h_act);
0287 
0288 % set u values y/hu, adjust KV values to yield the right flux
0289 
0290 network_act.kinetics.type = 'ms';
0291 network_act.kinetics.u    = y_act ./ hu_act;
0292 network_act.kinetics.c    = c_act;
0293 network_act.kinetics.KA   = KA_act;
0294 network_act.kinetics.KI   = KI_act;
0295 network_act.kinetics.KM   = KM_act;
0296 network_act.kinetics.KV   = ones(nr_act,1);
0297 network_act.kinetics.Keq  = exp(N_act' * [log(c_act)-1/RT*mu_act]);
0298 network_act.kinetics.h    = h_act;
0299 
0300 vv = network_velocities(c_act, network_act);
0301 
0302 if find(vv ==0), error('zero flux encountered'); end
0303 if find(vv.*v_act<0), error('wrong flux direction'); end
0304 
0305 network_act.kinetics.KV = v_act ./ vv .* network_act.kinetics.KV;
0306 
0307 % check
0308 % v_act_kinetic = network_velocities(c_act, network_act); [v_act_kinetic,v_act]
0309 
0310 u_act = network_act.kinetics.u;
0311 
0312 
0313 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0314 % check flux benefit response coefficients (should be equal to hu)
0315 % R_act    = basic_control_analysis(network_act,c_act,struct('used',ones(length(v_act),1)));
0316 % Rfvu_act = [ cba_constraints.zv(ind_active)' * R_act.RJ ]';
0317 % Rfcu_act = [ -[Q_act(ind_int_act)./c_act(ind_int_act)]' * R_act.RS(ind_int_act,:) ]';
0318 % ES_act1 = diag(1./v_act) * R_act.epsilon_1 * diag(c_act);
0319 % figure(18); plot(cba_constraints.hu(ind_active), [Rfvu_act + Rfcu_act] .* [v_act~=0],'.');
0320 
0321 % check investment condition:
0322 % zc_int_act_updated = -Q_act(ind_int_act)./c_act(ind_int_act);
0323 % E_T_act = diag(1./[zeta_act-1]) * [diag(zeta_act) * Mplus_act - Mminus_act];
0324 % ES_act  = E_T_act - beta_M_act .* M_act + alpha_A_act .* Wplus_act - beta_I_act .* Wminus_act;
0325 % ES_act_unscaled = diag(v_act) * ES_act  * diag(1./c_act);
0326 % LHS     = - L_act' *  zc_int_act_updated;
0327 % RHS = L_act_scaled' * ES_act(:,ind_int_act)' * y_act
0328 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0329 
0330 
0331 % ----------------------------------------------------------------------
0332 % Return to entire model with inactive reactions -> compute w and u
0333 % construct the kinetics for the entire model ...
0334 
0335 % set all metabolites from the inactive subnetwork = external
0336 
0337 u = 0 * v;    u(ind_active)     = u_act;  
0338 Q = nan * c;  Q(ind_met_active) = Q_act;  
0339 
0340 em = sparse(nr,nm);
0341 
0342 % for metabolites in inactive part of the network:
0343 % remove allosteric interactions, but keep them as reactants (KM values)
0344 
0345 network.kinetics                 = set_kinetics(network,'ms');
0346 network.kinetics.u               = u;  
0347 network.kinetics.c               = c;
0348 network.kinetics.KA              = em;
0349 network.kinetics.KI              = em;
0350 network.kinetics.KA(ind_active, ind_met_active) = network_act.kinetics.KA ;  
0351 network.kinetics.KI(ind_active, ind_met_active) = network_act.kinetics.KI ;  
0352 network.kinetics.KM(ind_active, ind_met_active) = network_act.kinetics.KM ;
0353 network.kinetics.KV(ind_active)  = network_act.kinetics.KV;
0354 network.kinetics.Keq(ind_active) = network_act.kinetics.Keq;  
0355 network.kinetics.h(ind_active)   = network_act.kinetics.h;  
0356 
0357 
0358 % ----------------------------------------------------------------------
0359 % set all metabolites from inactive subnetwork external
0360 
0361 ind_met_nonactive = setdiff(1:nm,ind_met_active);
0362 network.external(ind_met_nonactive) = 1;
0363 network.regulation_matrix(:,ind_met_nonactive) = 0;
0364 
0365 all_zext = zeros(nm,1);
0366 all_zext(ind_ext) = cba_constraints.z_ext;
0367 cba_constraints.z_ext = all_zext(find(network.external));
0368 
0369 % check: v_kinetic = network_velocities(c,network); [v_kinetic,v]
0370 
0371 
0372 % ----------------------------------------------------------------------
0373 % .. compute the enzyme response coefficients and ..
0374 
0375 R    = basic_control_analysis(network,c);
0376 E_sc = diag(1./v) * R.epsilon_1 * diag(c); % scaled elasticities
0377 
0378 % % check elasticities:
0379 % % 1. from reconstructed ms rate law
0380 % E_sc_act_1 = E_sc(ind_active,ind_met_active);
0381 % % 2. from original beta values
0382 % E_T_act    = diag(1./[zeta_act-1]) * [diag(zeta_act) * Mplus_act - Mminus_act];
0383 % E_sc_act_2 = E_T_act - beta_M_act .* M_act + alpha_A_act .* Wplus_act - beta_I_act .* Wminus_act;
0384 
0385 
0386 % ----------------------------------------------------------------------
0387 % check the balance equations again
0388 
0389 % economic balance equation (left hand side, right hand side)
0390 LHS = v.*[network.N'*w + cba_constraints.z_int];
0391 RHS = y;
0392 mismatch = norm(LHS-RHS);
0393 display(sprintf('  Economic balance equation:   mismatch %f',mismatch))
0394 if mismatch>10^-5,  [ LHS, RHS, LHS-RHS]
0395 end
0396 
0397 % investment balance equation (left hand side, right hand side)
0398 % first external, then internal metabolites
0399 % (only metabolites in active subnetwork are considered)
0400 
0401 LHS = [Q(ind_met_active(ind_ext_act)); ...
0402        L_act' * Q(ind_met_active(ind_int_act));];
0403 
0404 RHS = [E_sc(ind_active,ind_met_active(ind_ext_act))' * y(ind_active); ...
0405        L_act_scaled' * E_sc(ind_active,ind_met_active(ind_int_act))' * y(ind_active)];
0406 
0407 mismatch = norm(LHS-RHS);
0408 display(sprintf('  Investment balance equation: mismatch %f',mismatch))
0409 if mismatch>10^-5, 
0410   [ LHS, RHS, LHS-RHS]
0411 end
0412 
0413 
0414 % ----------------------------------------------------------------------
0415 % .. check the steady state for stability;
0416 % - eigenvalues must not have a positive real part
0417 % - all complex eigenvalues must have a negative real part
0418 
0419 M_eigenvalues = eig(full(R.M));
0420 is_stable = [max(real(M_eigenvalues))<=0] * [sum(real(M_eigenvalues(find(imag(M_eigenvalues))))==0)==0];
0421 
0422 switch is_stable,
0423   case 1, display('  o The steady state is stable');
0424   case 0, warning('  o The steady state is unstable; please sample again!');
0425 end
0426 
0427 
0428 % ----------------------------------------------------------------------
0429 % .. check if all active enzymes have a positive influence on the benefit
0430 
0431 zc_updated = - Q ./ c;
0432 zc_updated(find(network.external)) = 0;
0433 zc_updated(isnan(zc_updated)) = 0;
0434 
0435 Rfvu = [ cba_constraints.zv' * R.RJ ]';
0436 Rfcu = [         zc_updated' * R.RS ]';
0437 Rfu  = Rfvu + Rfcu;
0438 
0439 switch min(sign(Rfu(ind_active))),
0440   case 1,    display('  o All active enzymes have positive marginal benefits');
0441   otherwise, error('  Active enzyme with non-positive marginal benefit encountered');
0442 end
0443 
0444 
0445 % --------------------------------------------------------------------------------
0446 % check eigenvalues of fitness curvature matrix
0447 % - for an enzyme-econonic state, the matrix needs to be negative semidefinite
0448 % - to get unique solutions for differential expression prediction,
0449 %   it needs to be negative definite
0450 
0451 if cba_options.check_curvatures,
0452   %% FIXES NEEDED:
0453   %% HERE THE FIRST ORDER IS COMPUTED FOR THE SECOND TIME .. MAYBE COMPUTE ONLY ONCE?????
0454   %% COMPUTING THE TENSOR TAKES LONG .. JUST COMPUTE THE RELEVANT PART!!!
0455   [Ec,Eu,parameters,Ecc,Ecu,Euu] = elasticities(network,c,struct('only_enzyme_levels',1));
0456   [CJ, CS]     = control_coefficients(network.N, Ec,network.external);
0457   [RSu, RJu, RSuu, RJuu] = response_coefficients(CS, Ec, Eu, Ecc, Ecu, Euu);
0458   f_u_active = [cba_constraints.zv' * RJu(1:nr,ind_active)]';
0459   f_uu_active = squeeze(tensor_product(cba_constraints.zv', RJuu(1:nr,ind_active,ind_active)));
0460 
0461   P_orth = eye(length(ind_active)) - 1/[f_u_active'*f_u_active] * [f_u_active*f_u_active'];
0462   
0463   %% check eigenvalues in subspace orthogonal on f_u??  IS THAT CORRECT AT ALL??
0464   [f_uu_active_eigenvectors, f_uu_active_eigenvalues] = eig(full(  P_orth' * f_uu_active * P_orth));
0465   [f_uu_active_eigenvalues,order] = sort(diag(f_uu_active_eigenvalues));
0466   f_uu_eigenvectors = f_uu_active_eigenvectors(:,order);
0467   display(sprintf('  %d directions with positive benefit curvature',sum(f_uu_active_eigenvalues>0)));
0468   
0469   figure(1); netgraph_concentrations(network,[],network.kinetics.u,1,struct('arrowstyle','none'));
0470   title('Enzyme levels');
0471   show_cc = nan * ones(nr,1); show_cc(ind_active) = f_u_active;
0472   figure(2); netgraph_concentrations(network,[],network.kinetics.u,1,struct('arrowstyle','none'));
0473   title('Enzyme benefit gradient');
0474   show_cc = nan * ones(nr,1); show_cc(ind_active) = f_uu_eigenvectors(:,end);
0475   figure(3); netgraph_concentrations(network,[],  show_cc,1,struct('arrowstyle','none'));
0476   title('Curvature eigenvector');
0477 end
0478 
0479 
0480 % ---------------------------------------------
0481 % output data structure 'res'
0482 
0483 u(v==0) = 0;
0484 
0485 res.w   = w;
0486 res.u   = u;
0487 res.Q   = Q;
0488 res.fc_updated = zc_updated;
0489 res.kinetics = network.kinetics;
0490 res.R   = R;
0491 res.Rfu = Rfu;
0492 res.Rfvu = Rfvu;
0493 res.Rfcu = Rfcu;
0494 
0495 cba_constraints.zc    = fc_updated;
0496 cba_constraints.Q_ext = Q(find(network.external));