cba_economic_state

 CBA_ECONOMIC_STATE - Find fluxes, chemical and economic potentials satisfying FBA/EBA/CBA constraints

 result = cba_economic_state(N, external, cba_constraints, cba_options)

 Determine feasible flux vector v, chemical potential differences dmu, 
 and rate values y under EBA and CBA constraints with either
  - linear FBA objective (zv' * v) 
  - sum of squared residuals objective sum([[v-v_mean]/v_std].^2)
 plus a small Euclidean norm term for regularisation

 Output:
   result.v      - Flux mode
   result.dmu    - Chemical potential differences
   result.y      - (only exists if cba_options.cba_conditions == 'y')
   result.w      - (only exists if cba_options.cba_conditions == 'w')
   result.value  - Objective function

0001 function result = cba_economic_state(N, external, cba_constraints, cba_options)
0002 
0003 % CBA_ECONOMIC_STATE - Find fluxes, chemical and economic potentials satisfying FBA/EBA/CBA constraints
0004 %
0005 % result = cba_economic_state(N, external, cba_constraints, cba_options)
0006 %
0007 % Determine feasible flux vector v, chemical potential differences dmu,
0008 % and rate values y under EBA and CBA constraints with either
0009 %  - linear FBA objective (zv' * v)
0010 %  - sum of squared residuals objective sum([[v-v_mean]/v_std].^2)
0011 % plus a small Euclidean norm term for regularisation
0012 %
0013 % Output:
0014 %   result.v      - Flux mode
0015 %   result.dmu    - Chemical potential differences
0016 %   result.y      - (only exists if cba_options.cba_conditions == 'y')
0017 %   result.w      - (only exists if cba_options.cba_conditions == 'w')
0018 %   result.value  - Objective function
0019 
0020 
0021 % The constraints given by cba_constraints.mu_fix, cba_constraints.dmu_fix are ignored
0022 
0023 
0024 %  ------------------------------------------------------------
0025 % initialisation
0026 
0027 [nm,nr]      = size(N);
0028 
0029 epsilon = 10^-3;
0030 
0031 cba_constraints.v_min(find(cba_constraints.v_sign==1))  = epsilon;
0032 cba_constraints.v_max(find(cba_constraints.v_sign==-1)) = -epsilon;
0033 
0034 if isempty(cba_constraints.dmu_min),  cba_constraints.dmu_min = -10*ones(nr,1); end
0035 if isempty(cba_constraints.dmu_max),  cba_constraints.dmu_max =  10*ones(nr,1); end
0036 if isempty(cba_constraints.y_min),    cba_constraints.y_min   = -10*ones(nr,1); end
0037 if isempty(cba_constraints.y_max),    cba_constraints.y_max   =  10*ones(nr,1); end
0038 if isempty(cba_constraints.w_min),    cba_constraints.w_min   = -10*ones(nm,1); end
0039 if isempty(cba_constraints.w_max),    cba_constraints.w_max   =  10*ones(nm,1); end
0040 
0041 zv           = cba_constraints.zv;
0042 ind_v_fix    = find(isfinite(cba_constraints.v_fix));
0043 ind_extsigns = find(isfinite(cba_constraints.ext_sign));
0044 n_int        = sum(external==0); 
0045 N_int        = N(find(external==0),:);
0046 Ktot         = analyse_N(N);
0047 Kint         = analyse_N(N(find(external==0),:));
0048 nKint        = size(Kint,2);
0049 nKtot        = size(Ktot,2);
0050 my_eye       = eye(nr);
0051 regularisation_weight  = 0; % weights for regularisation term
0052 
0053 
0054 % ------------------------------------------------------------
0055 % vector x contains all (v, dmu, y) or (v, dmu, w)
0056 % indices within vector (either ind_y or ind_w is unnecessary)
0057 
0058 ind_v   =           1:nr;
0059 ind_dmu =     nr + (1:nr);
0060 ind_y   = 2 * nr + (1:nr);
0061 ind_w   = 2 * nr + (1:nm);
0062 
0063 
0064 % ------------------------------------------------------
0065 % cba_constraints:
0066 %
0067 % (where v = x(ind_v), dmu = x(ind_dmu), y = x(ind_y), w = x(ind_w) )
0068 %
0069 % limits           : cba_constraints.v_min   <=  x(ind_v)   <= cba_constraints.v_max
0070 %                  : cba_constraints.dmu_min <=  x(ind_dmu) <= cba_constraints.dmu_max
0071 %                  : cba_constraints.y_min   <=  x(ind_y)   <= cba_constraints.y_max
0072 % regularisation   : sum(x.^2) = min !
0073 %
0074 % fixed fluxes     : x(ind_v_fix)          = cba_constraints.v_fix(ind_v_fix)
0075 % stationary fluxes: N_int * x(ind_v )   = 0
0076 % positive objective: zv' * x(ind_v )   > 0
0077 %
0078 % eba, wegscheider : Ktot' * x(ind_dmu)  = 0
0079 % eba signs        : diag(x(ind_v)) * x(ind_dmu) <= - epsilon < 0   where v == 0
0080 %
0081 % variant 1: use rate values y
0082 %   y balance      : Kint' * x(ind_y)    = Kint' * zv
0083 %   y signs        : diag(x(ind_v)) * x(ind_y) >= epsilon > 0   where v == 0
0084 %
0085 % variant 2: use metabolite values w
0086 %   w signs        : diag(x(ind_v)) * [ N' * x(ind_w) + zv] >= epsilon > 0 if v == 0
0087 
0088 % ------------------------------------------------------
0089 % rewrite this in the form
0090 %
0091 % Aeq * x  = beq
0092 %   A * x <= b
0093 %       C  = cbaf2(x,N)   (for signs)
0094 %     min != cbaf1(x)     (regularisation)
0095 
0096 
0097 switch cba_options.cba_conditions, 
0098   
0099   case 'y', 
0100 
0101     display('Optimisation with rate values y');
0102 
0103     %% Equality cba_constraints:
0104     %%  v_fix
0105     %%  stationarity
0106     %%  wegscheider for dmu
0107     %%  balance for y
0108     
0109     Aeq  = double([my_eye(ind_v_fix,:), zeros(length(ind_v_fix),2*nr)      ; ...
0110                    N_int,               zeros(n_int,2*nr)                  ; ...
0111                    zeros(nKtot,nr),     Ktot',            zeros(nKtot,nr)  ; ...
0112                    zeros(nKint,2*nr),                     Kint'             ]);
0113 
0114     beq  = double([cba_constraints.v_fix(ind_v_fix);... 
0115                    zeros(n_int,1)              ;...
0116                    zeros(nKtot,1)              ;...
0117                    Kint' * zv ]);
0118 
0119     xmin   = [cba_constraints.v_min; cba_constraints.dmu_min; cba_constraints.y_min];
0120     xmax   = [cba_constraints.v_max; cba_constraints.dmu_max; cba_constraints.y_max];
0121     zvfull = [cba_constraints.zv ; zeros(nr,1); zeros(nr,1)];
0122 
0123     %% inequality cba_constraints:
0124     %% positive objective
0125     %% external signs
0126     
0127     A = [-zv', zeros(1,2*nr) ; ...
0128          double([-diag(cba_constraints.ext_sign(ind_extsigns)) * N(ind_extsigns,:), ...
0129                  zeros(length(ind_extsigns),2*nr)])];
0130     
0131     b = [0; zeros(length(ind_extsigns),1)];
0132     
0133   case 'w',  
0134 
0135     display('Optimisation with metabolite values w');
0136 
0137     %% Equality cba_constraints:
0138     %%  v_fix
0139     %%  stationarity
0140     %%  wegscheider for dmu
0141     
0142     Aeq  = double([my_eye(ind_v_fix,:),    zeros(length(ind_v_fix),nr+nm) ; ...
0143                    N_int,                 zeros(n_int,nr+nm)              ; ...
0144                    zeros(nKtot,nr),       Ktot',          zeros(nKtot,nm)]);
0145     
0146     beq  = double([cba_constraints.v_fix(ind_v_fix);... 
0147                    zeros(n_int,1)           ;...
0148                    zeros(nKtot,1)          ]);
0149 
0150     xmin   = [cba_constraints.v_min; cba_constraints.dmu_min; cba_constraints.w_min];
0151     xmax   = [cba_constraints.v_max; cba_constraints.dmu_max; cba_constraints.w_max];
0152     zvfull = [cba_constraints.zv; zeros(nr,1); zeros(nm,1)];
0153 
0154     
0155     %% inequality cba_constraints:
0156     %% positive objective
0157     %% external signs
0158     
0159     A = [zv', zeros(1,nr+nm) ; ...
0160          double([-diag(cba_constraints.ext_sign(ind_extsigns)) * N(ind_extsigns,:), ...
0161                  zeros(length(ind_extsigns),nr+nm)])];
0162     
0163     b = [0; zeros(length(ind_extsigns),1)];
0164     
0165     %% inequality cba_constraints:
0166     %% positive objective
0167     %% external signs
0168     
0169     A = [-zv', zeros(1,nr+nm) ; ...
0170          double([-diag(cba_constraints.ext_sign(ind_extsigns)) * N(ind_extsigns,:), ...
0171                  zeros(length(ind_extsigns),nr+nm)])];
0172     
0173     b = [0; zeros(length(ind_extsigns),1)];
0174     
0175 end
0176 
0177 
0178 % ----------------------------------------------------------
0179 % optimisation
0180 
0181 ppp.Ntrans  = full(N');
0182 ppp.ind_v   = ind_v;
0183 ppp.ind_dmu = ind_dmu;
0184 ppp.ind_y   = ind_y;
0185 ppp.ind_w   = ind_w;
0186 ppp.zv      = zv;
0187 
0188 opt = optimset('MaxFunEvals',10000);
0189 
0190 switch cba_options.objective
0191   
0192   case 'fba',
0193 
0194     switch cba_options.cba_conditions
0195       case 'y', %% use constraint function cba2f
0196         x0 = rand(3*nr,1)-.5;
0197         weights = regularisation_weight * ones(size(x0));
0198         ppp.epsilon = -0.1;
0199         [xtry,negvalue,exitflag] = fmincon(@(xx)cbaf1(xx,zvfull,weights), x0, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2(xx,ppp), opt);
0200         %% run, with slightly relaxed feasibility criterion
0201         ppp.epsilon = -0.001;
0202         [x,negvalue,exitflag] = fmincon(@(x)cbaf1(x,zvfull,weights), xtry, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2(xx,ppp), opt);
0203       case 'w', %% use constraint function cba2fw
0204         x0 = rand(2*nr+nm,1)-.5;
0205         weights = regularisation_weight * ones(size(x0));
0206         ppp.epsilon = -0.1;
0207         [xtry,negvalue,exitflag] = fmincon(@(xx)cbaf1(xx,zvfull,weights), x0, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2w(xx,ppp), opt);
0208         %% run, with slightly relaxed feasibility criterion
0209         ppp.epsilon = -0.001;
0210         [x,negvalue,exitflag] = fmincon(@(x)cbaf1(x,zvfull,weights), xtry, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2w(xx,ppp), opt);
0211     end
0212     
0213   case 'fit',
0214 
0215     v_mean = cba_constraints.v_mean;
0216     v_std  = cba_constraints.v_std;
0217     switch cba_options.cba_conditions
0218       case 'y', %% use constraint function cba2f
0219         x0 = rand(3*nr,1)-.5;
0220         ppp.epsilon = -0.1;
0221         [xtry,value,exitflag] = fmincon(@(xx)cbaf3(xx,v_mean,v_std), x0, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2(xx,ppp), opt);
0222         %% run, with slightly relaxed feasibility criterion
0223         ppp.epsilon = -0.001;
0224         [x,value,exitflag] = fmincon(@(x)cbaf3(x,v_mean,v_std), xtry, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2(xx,ppp), opt);
0225       case 'w', %% use constraint function cba2fw
0226         x0 = rand(2*nr+nm,1)-.5;
0227         weights = regularisation_weight * ones(size(x0));
0228         ppp.epsilon = -0.1;
0229         [xtry,value,exitflag] = fmincon(@(xx)cbaf3(xx,v_mean,v_std), x0, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2w(xx,ppp), opt)
0230         %% run, with slightly relaxed feasibility criterion
0231         ppp.epsilon = -0.001;
0232         [x,value,exitflag] = fmincon(@(x)cbaf3(x,v_mean,v_std), xtry, A, b, Aeq, beq, xmin, xmax, @(xx)cbaf2w(xx,ppp), opt)
0233     end
0234     
0235 end
0236 
0237 % ----------------------------------------------------------
0238 % format the output
0239 
0240 v     = nan * ones(nr,1);
0241 dmu   = nan * ones(nr,1);
0242 y     = nan * ones(nr,1);
0243 w     = nan * ones(nm,1);
0244 value = nan;
0245 
0246 if exitflag<=0,
0247   display('CBA failed');
0248 else,
0249   v   = x(ind_v);
0250   dmu = x(ind_dmu);
0251   switch cba_options.objective
0252     case 'fba',
0253       value = -negvalue; 
0254   end
0255   switch cba_options.cba_conditions,
0256     case 'y', y = x(ind_y);
0257     case 'w', w = x(ind_w);
0258   end
0259 
0260   if value < 0, warning('Objective function is negative'); end
0261   if sum(v.*dmu > 0), warning('Solution violates EBA constraint'); end
0262   switch cba_options.cba_conditions,
0263     case 'y', 
0264       if sum(v.*y < 0), warning('Solution violates CBA constraint'); end
0265     case 'w',
0266       if sum(v.*[N'*w+zv] < 0), warning('Solution violates CBA constraint'); end
0267   end
0268   
0269 end
0270 
0271 result.value = value; 
0272 result.v     = v;
0273 result.dmu   = dmu;
0274 
0275 switch cba_options.cba_conditions,
0276   case 'y',    result.y  = y;
0277   case 'w',    result.w  = w;
0278 end