cba_reconstruct_model

 CBA_RECONSTRUCT_MODEL - Build model from economical flux profile

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