model_embedding_consistent_fluxes

 [collect_v_all,collect_v_kinetic,c_stat,v_stat,collect_c,collect_v,collect_mu,mu_stat] = ...
    model_embedding_consistent_fluxes(network,kinetic_models,me_options,mapping_metabolites,mapping_reactions)

 o determine concentrations and fluxes within the kinetic models
 o determine equilibrium concentrations within kinetic models
 o kinetic models are considered one after the other (corresponding to their
   priority order: each model can predefine concentrations for the subsequent models!

0001 function [collect_v_all,collect_v_kinetic,c_stat,v_stat,collect_c,collect_v,collect_mu,mu_stat] = model_embedding_consistent_fluxes(network,kinetic_models,me_options,mapping_metabolites,mapping_reactions)
0002 
0003 % [collect_v_all,collect_v_kinetic,c_stat,v_stat,collect_c,collect_v,collect_mu,mu_stat] = ...
0004 %    model_embedding_consistent_fluxes(network,kinetic_models,me_options,mapping_metabolites,mapping_reactions)
0005 %
0006 % o determine concentrations and fluxes within the kinetic models
0007 % o determine equilibrium concentrations within kinetic models
0008 % o kinetic models are considered one after the other (corresponding to their
0009 %   priority order: each model can predefine concentrations for the subsequent models!
0010 
0011 % collect_v_all
0012 % collect_v_kinetic
0013 %
0014 % values in individual kinetic models
0015 % (list of vectors, one for each kinetic model; vector sizes like in individual kinetic models)
0016 %  c_stat
0017 %  v_stat
0018 %  mu_stat
0019 %
0020 % values in global network; computed by going through kinetic models iteratively;
0021 % each time, some concentrations may be imposed by previous kinetic models.
0022 %  collect_c
0023 %  collect_v
0024 %  collect_mu
0025 
0026 % values resulting from each model's own concentrations
0027 collect_c_own= nan * ones(length(network.metabolites),length(kinetic_models));
0028 collect_v_own= nan * ones(length(network.actions),    length(kinetic_models));
0029 
0030 % values possibly changed by imposed concentrations
0031 collect_c    = nan * ones(length(network.metabolites),length(kinetic_models));
0032 collect_v    = nan * ones(length(network.actions),    length(kinetic_models));
0033 collect_c_eq = nan * ones(length(network.metabolites),length(kinetic_models));
0034 collect_mu   = nan * ones(length(network.metabolites),length(kinetic_models));
0035 
0036 % collects all values fixed so far
0037 collect_c_all    = nan * ones(length(network.metabolites),1);
0038 collect_c_eq_all = nan * ones(length(network.metabolites),1);
0039 collect_v_all    = nan * ones(length(network.actions),1);
0040 collect_v_kinetic= nan * ones(length(network.actions),1);
0041 
0042 % insert fluxes for network (if given)
0043 if isfield(me_options,'v_network'),
0044   collect_v_all(1:length(me_options.v_network)) = me_options.v_network;
0045 end
0046 
0047 for it = 1:length(kinetic_models),
0048 
0049   %% set "_own" variables:
0050   %% concentrations and rates as given by kinetic models individually
0051   nn             = kinetic_models{it};
0052   c_stat_own{it} = nn.s_init;
0053   v_stat_own{it} = network_velocities(c_stat_own{it},nn);
0054   collect_c_own(mapping_metabolites{it},it) = c_stat_own{it};
0055   collect_v_own(mapping_reactions{it},it)   = v_stat_own{it};
0056 
0057   %% set "_stat" variables:
0058   %% ... but some concentrations get fixed by preceding models!
0059   ind_predefined_c = find(isfinite(collect_c_all(mapping_metabolites{it})));
0060   predefined_c     = collect_c_all(mapping_metabolites{it}(ind_predefined_c));
0061   predefined_c_eq  = collect_c_eq_all(mapping_metabolites{it}(ind_predefined_c));
0062   nn.external(ind_predefined_c) = 1;
0063   c_stat{it} = c_stat_own{it};
0064   c_stat{it}(ind_predefined_c) = predefined_c;
0065   v_stat{it} = network_velocities(c_stat{it}, kinetic_models{it});
0066 
0067   if me_options.enforce_stationarity,
0068     c_stat{it}  = network_steady_state(nn,c_stat{it},10000);
0069     v_stat{it}  = network_velocities(c_stat{it}, kinetic_models{it});
0070     v_stat{it}  = es_make_fluxes_stationary(kinetic_models{it},  v_stat{it});
0071     v_stat{it}  = es_make_fluxes_stationary(kinetic_models{it},  v_stat{it});
0072   end
0073 
0074   %% compute equilibrium concentrations (only possible for reversible rate laws!!!
0075   nnn          = nn;
0076   nnn.external = 0 * nnn.external;
0077   c_eq{it}     = c_stat{it};
0078   c_eq{it}(ind_predefined_c) = predefined_c_eq;
0079   c_eq{it}    = real(network_steady_state(nnn, c_eq{it},10));%000);
0080 
0081   %% for safety reasons (in case of irreversible rate laws)
0082   c_eq{it} = real(c_eq{it});
0083   c_eq{it}(c_eq{it}<10^-15) = 10^-15;
0084   
0085   %% compare rate signs and driving forces
0086   %[ sign(v_stat{it}) .* sign(-kinetic_models{it}.N'*log(c_stat{it}./c_eq{it}))]
0087   mu_stat{it}   = RT * log(c_stat{it}./c_eq{it});
0088 
0089   collect_c(mapping_metabolites{it},it)    = c_stat{it};
0090   collect_v(mapping_reactions{it},it)      = v_stat{it};
0091   collect_c_eq(mapping_metabolites{it},it) = c_eq{it};
0092   collect_mu(mapping_metabolites{it},it)   = mu_stat{it};
0093   
0094   ind_undet                   = find(isnan(collect_c_all));
0095   collect_c_all(ind_undet)    = collect_c(ind_undet,it);
0096   collect_c_eq_all(ind_undet) = collect_c_eq(ind_undet,it);
0097   ind_undet                   = find(isnan(collect_v_all));
0098   collect_v_all(ind_undet)    = collect_v(ind_undet,it);
0099   ind_det                     = find(isfinite(collect_v(:,it)));
0100   collect_v_kinetic(ind_det)  = collect_v(ind_det,it);
0101 
0102 end
0103 
0104 % ---------------------------------------
0105 % if there are overlaps between kinetic models and if their
0106 % concentrations or equilibrium concentrations differ -> error!
0107 
0108 for it = 1:length(kinetic_models);
0109   if find([c_eq{it}==0]+[isfinite(c_eq{it})==0]), warning('Irreversible reaction encountered'); end
0110 end
0111 
0112 for it = 1:length(network.metabolites),
0113   dum = isfinite(collect_c_eq(it,:));
0114   if length(unique(collect_c_eq(it,dum))) > 1, 
0115     warning(sprintf('Inconsistent equilibrium concentration of metabolite %s',network.metabolites{it})); 
0116   end
0117 end