network_steady_state.m

％The Metabolic Networks Toolbox contains functions to create, ％modify, display, and simulate bioche

% [S,J,Sdot] = network_steady_state(network,s,integrate,L_int,NR_int,indep_met_int)
%
% calculate steady state concentrations S, fluxes J, and elasticities epsilon
% given network and parameters initial metabolite concentrations s (column vector)

function  [S,J,Sdot] = network_steady_state(network,s,integrate,L_int,NR_int,indep_met_int)

 if ~exist('s','var'), s=ones(size(network.metabolites)); end
 if ~exist('integrate','var'), integrate=0; end

 [n_S,n_A] = size(network.N);

 external  = find(network.external);
 internal  = setdiff(1:n_S,external);
 
 s_ext       = s(external);
 s_int       = s(internal);
 N_int       = network.N(internal,:);

 if ~exist('L_int','var'), 
 % fprintf('Computing conservation relations...\n');
 [L_int, NR_int, indep_met_int] = reduce_N(N_int);
 end

% let the network run for some time to get a better initial estimate
 
 if integrate,
   [t, s_int_t,s_t,met_int] = network_integrate(network, s,1);
   plot(t,s_t);
   s_int = s_int_t(:,end);
 end

 s_ind = s_int(indep_met_int);
 T = s_int-L_int*s_ind;

% fprintf('Computing steady state...\n');
 
 S_ind = abs(fminsearch(@calc_sdot_sq,s_ind,...
	   optimset('TolFun',0.005,'MaxFunEvals',1000000,'MaxIter',1000),...
			      network,external,internal,s_ext,L_int,T,N_int)) ;
 S           = zeros(n_S,1);
 S(internal) = L_int*abs(S_ind)+T;
 S(external) = s_ext;
 
 J            = network_velocities(S,network);
 Sdot         = N_int * J;

% --------------------------------------------------------------- 
 
function sdot_sq = calc_sdot_sq(s_ind,network,external,internal,s_ext,L,T,N_int)

  s            = zeros(size(network.metabolites));
  s(external)  = s_ext;
  s(internal)  = L * abs(s_ind) + T;
  v            = network_velocities(s,network);
  sdot_sq      = sum( (N_int * v).^2) + 100*(min(L * abs(s_ind) + T)<0);