model_embedding

[network,network_CoHid,network_combined, res] = model_embedding(kinetic_models, network, network_CoHid, me_options)

 Inputs
   kinetic_models
   network
   network_CoHid
   me_options
     me_options.v_network
     me_options.id   structure describing mappings between the models. 
                                     Default:
       me_options.id.metabolites_network        = 'metabolites';
       me_options.id.reactions_network          = 'actions';
       me_options.id.metabolites_kinetic_models = {'metabolites', ..};
       me_options.id.reactions_kinetic_models   = {'actions', ..};
     me_options.fba_constraints   (refer to the network model; are extended to augmented network below)

Remarks:
 kinetic models in list 'kinetic models' are ordered by priority: highest priority first
 if the network is automatically extended, the FBA constraints have to be given for the 
 final extended network

 Outputs
   network
   network_CoHid
   network_combined
   res

0001 function [network_combined, res,network_aug,network_aug_CoHid,network_aug_CoSplit] = model_embedding(kinetic_models, network, network_CoHid, me_options)
0002 
0003 %[network,network_CoHid,network_combined, res] = model_embedding(kinetic_models, network, network_CoHid, me_options)
0004 %
0005 % Inputs
0006 %   kinetic_models
0007 %   network
0008 %   network_CoHid
0009 %   me_options
0010 %     me_options.v_network
0011 %     me_options.id   structure describing mappings between the models.
0012 %                                     Default:
0013 %       me_options.id.metabolites_network        = 'metabolites';
0014 %       me_options.id.reactions_network          = 'actions';
0015 %       me_options.id.metabolites_kinetic_models = {'metabolites', ..};
0016 %       me_options.id.reactions_kinetic_models   = {'actions', ..};
0017 %     me_options.fba_constraints   (refer to the network model; are extended to augmented network below)
0018 %
0019 %Remarks:
0020 % kinetic models in list 'kinetic models' are ordered by priority: highest priority first
0021 % if the network is automatically extended, the FBA constraints have to be given for the
0022 % final extended network
0023 %
0024 % Outputs
0025 %   network
0026 %   network_CoHid
0027 %   network_combined
0028 %   res
0029 
0030 me_options_default = me_default_options(length(kinetic_models));
0031 me_options         = join_struct(me_options_default, me_options);
0032 if isempty(me_options.id),
0033   me_options.id = me_options_default.id;
0034 end
0035 
0036 % -------------------------------------------------------------------------------
0037 % Preparations
0038 
0039 % All compartments in kinetic models must also appear in network model
0040 % and have the same sizes
0041 
0042 display('  Checking the use of compartments');
0043 
0044 if ~isfield(network,'compartments'),
0045   network.compartments      = {};
0046   network.compartment_sizes = [];
0047 end
0048 
0049 for it = 1:length(kinetic_models),
0050   if ~isfield(kinetic_models{it},'compartments'),
0051     kinetic_models{it}.compartments = {};
0052     kinetic_models{it}.compartment_sizes = [];
0053   end
0054   missing_compartments = setdiff(kinetic_models{it}.compartments,network.compartments);
0055   if length(missing_compartments),
0056     error(sprintf('Some compartments from model %d are missing in network model',it));
0057   end
0058   ll = label_names(kinetic_models{it}.compartments,network.compartments);
0059   if sum(kinetic_models{it}.compartment_sizes(find(ll)) ~= network.compartment_sizes(ll(find(ll)))),
0060     kinetic_models{it}.compartment_sizes(find(ll))
0061     network.compartment_sizes(ll(find(ll)))
0062     error('Compartment sizes do not match between models');
0063   end
0064 end
0065 
0066 
0067 % -------------------------------------------------------------------------------
0068 % Phase 1: Map the kinetic models onto the network model
0069 
0070 % Match elements and, if necessary,
0071 % add metabolites and reactions to the network model
0072 
0073 [mapping_metabolites, mapping_reactions, covered_metabolites, covered_reactions, shared_metabolites, shared_reactions, network_aug, network_aug_CoHid,network_aug_CoSplit] = embedding_element_mapping(kinetic_models,network,me_options);
0074 
0075 ind_ext = label_names(me_options.set_external,network_aug.metabolites);
0076 if length(ind_ext),
0077 network_aug.external(label_names(me_options.set_external,network_aug.metabolites) ) = 1;
0078 network_aug.external(label_names(me_options.set_internal,network_aug.metabolites) ) = 0;
0079 end
0080 
0081 % -------------------------------------------------------------------------------
0082 % Phase 2: Determine consistent fluxes
0083 
0084 % o determine concentrations and fluxes within the kinetic models
0085 % o determine equilibrium concentrations within kinetic models
0086 % o kinetic models are considered one after the other (corresponding to their
0087 %   priority order: each model can predefine concentrations for the subsequent models!
0088 
0089 [collect_v_all,collect_v_kinetic,c_stat,v_stat,collect_c,collect_v,collect_mu,mu] = model_embedding_consistent_fluxes(network_aug,kinetic_models,me_options,mapping_metabolites,mapping_reactions);
0090 
0091 for it = 1:length(kinetic_models),
0092   kinetic_models{it}.s_init = c_stat{it};
0093 end
0094 
0095 %netgraph_concentrations(kinetic_models{1},log(c_stat{1}+1),[],1,struct('arrowvalues',v_stat{1},'arrowstyle','fluxes'))
0096 
0097 % ---------------------------------------
0098 % Project fluxes on network model, where fluxes are then constrained to fluxes on the kinetic model parts
0099 
0100 % Scale network fluxes to make them comparable to kinetic fluxes
0101 
0102 ind_v_kinetic = find(isfinite(collect_v_kinetic));
0103 ratio         = median(collect_v_kinetic(ind_v_kinetic) ./ collect_v_all(ind_v_kinetic) );
0104 
0105 v_mean = ratio * collect_v_all;
0106 v_std  = guess_flux_std(v_mean);
0107 v_sign = nan * v_mean;
0108 v_fix  = nan * v_mean;
0109 
0110 if length(me_options.fba_constraints),
0111   nr = length(me_options.fba_constraints.v_sign);
0112   if isstruct(me_options.fba_constraints),
0113     v_sign(1:nr) = me_options.fba_constraints.v_sign;
0114     v_fix(1:nr)  = me_options.fba_constraints.v_fix;
0115   end
0116 end
0117 
0118 % Impose fluxes from kinetic model(s) as fixed constraint
0119 v_fix(isfinite(collect_v_kinetic)) = collect_v_kinetic(isfinite(collect_v_kinetic));
0120 
0121 if ~isfield(network_aug,'reaction_names'),
0122   network_aug.reaction_names = network_aug.actions;
0123 end
0124 
0125 try 
0126   v_proj = project_fluxes(network_aug.N, find(network_aug.external), v_mean, v_std, v_sign, struct('method','euclidean'),v_fix);
0127 catch
0128   warning('Steady state fluxes of kinetic models cannot be embedded into the larger network; I am now trying to find an approximate solution');
0129   v_mean(find(isfinite(v_fix))) = v_fix(find(isfinite(v_fix)));
0130   v_std  = guess_flux_std(v_mean);
0131   v_proj = project_fluxes(network_aug.N, find(network_aug.external), v_mean, v_std, v_sign, struct('method','euclidean'));
0132 
0133 end
0134 
0135 % ---------------------------------------
0136 % Determine enyzme adjustments for the kinetic models
0137 % for each kinetic model: v[stationary within network] = enzyme_adjustment * v[original]
0138 
0139 for it = 1:length(kinetic_models);
0140   if sum(v_stat{it}==0),
0141     error('Vanishing fluxes in submodel');
0142   end
0143   enzyme_adjustment{it} = v_proj(mapping_reactions{it}) ./ v_stat{it};
0144   if min(abs(enzyme_adjustment{it}))<10^-5,
0145     error('Difference in fluxes too extreme for enzyme adjustment');
0146   end  
0147 end
0148 
0149 
0150 % -----------------------------------------------------------------------
0151 % Collect concentrations and chemical potentials from kinetic models
0152 
0153 c_fix  = nan * ones(size(network_aug.metabolites));
0154 
0155 for it = length(kinetic_models):-1:1;
0156   c_fix(mapping_metabolites{it}) = c_stat{it};
0157 end
0158 
0159 % Determine predefined dmu within each kinetic model;
0160 % copy them into a vector dmu_fix for the entire network
0161 
0162 dmu_fix = nan * ones(size(network_aug.actions));
0163 for it = 1:length(kinetic_models);
0164   mmm = mu{it};
0165   mmm(isnan(mmm)) = 0;
0166   my_dmu = kinetic_models{it}.N'* mmm;
0167   my_dmu(find(abs(kinetic_models{it}.N') * isnan(mu{it}))) = nan;
0168   dmu_fix(mapping_reactions{it}) = my_dmu;
0169 end
0170 
0171 % Check if the flux distribution is thermo-feasible, and determine a feasible dmu vector
0172 
0173 [eba_feasible, dmu] = eba_feasible_lp(v_proj, network_aug.N, dmu_fix,[],[],8); 
0174 
0175 if ~eba_feasible, 
0176   warning('Flux distribution is infeasible for given chemical potentials.');
0177   display('Workaround: Inventing chemical potential differences');
0178   %% ALTERNATIVE: INVENT dmu .. THIS HAS TO BE FIXED;
0179   %% PROBLEMS: original models may be infeasible;
0180   %% numeric calculation of stationary and equilibrium concentrations
0181   %% is not accurate enough to obtain dmu that agree with rate signs
0182   [eba_feasible, dmu] = eba_feasible_lp(v_proj, network_aug.N,[],100); 
0183   if ~eba_feasible,
0184     warning('Infeasible flux distribution');
0185     display('Next workaround: Inventing chemical potential differences');
0186     dmu = -sign(v_proj);
0187   end
0188 end
0189 
0190 
0191 % -------------------------------------------------------------------------------
0192 % Phase 3: Insert standard rate laws and build the combined model
0193 
0194 % Problem to be solved: predefine the c and mu values
0195 % of communicating metabolites in elasticity sampling
0196 
0197 % run elasticity sampling
0198 
0199 [es_options, es_constraints] = es_default_options(network_aug);
0200 
0201 es_constraints.v_fix         = v_proj; 
0202 es_constraints.log_c_fix     = log(c_fix); 
0203 es_constraints.dmu_fix       = dmu;
0204 es_constraints.dmu_limit     = 10;      
0205 es_constraints.dmu_limit_min = 10^-5;
0206 
0207 result = es_reference_state(network_aug, es_options, es_constraints);
0208 
0209 if me_options.verbose,
0210   display('  Check: Predefined and resulting concentrations');
0211   [c_fix,  result.c]
0212   display('  Check: Predefined and resulting chem. pot differences');
0213   [dmu,   -result.A]
0214 end
0215 
0216 
0217 % -------------------------------------------------------------------------------
0218 % Outputs
0219 
0220 network_combined                              = network_aug;
0221 network_combined.kinetics                     = struct;
0222 network_combined.kinetics.type                = 'embedded_kinetic_models';
0223 network_combined.kinetics.kinetic_models      = kinetic_models;
0224 network_combined.kinetics.kinetics_network    = result.kinetics;
0225 network_combined.kinetics.mapping_metabolites = mapping_metabolites;
0226 network_combined.kinetics.mapping_reactions   = mapping_reactions;
0227 network_combined.kinetics.enzyme_adjustment   = enzyme_adjustment;
0228 network_combined.s_init                       = result.c;
0229 
0230 res.v_mean              = v_mean;
0231 res.c_fix               = c_fix;
0232 res.dmu_fix             = dmu_fix;
0233 res.mu                  = result.mu;
0234 res.A                   = result.A;
0235 res.c_stat              = c_stat;
0236 res.v_stat              = v_stat;
0237 res.v_proj              = v_proj;
0238 res.c_combined          = network_combined.kinetics.kinetics_network.c;
0239 res.v_combined          = network_velocities(res.c_combined,network_combined);
0240 res.collect_c           = collect_c;
0241 res.collect_v           = collect_v;
0242 res.collect_mu          = collect_mu;
0243 res.covered_reactions   = covered_reactions;
0244 res.covered_metabolites = covered_metabolites;
0245 res.shared_metabolites  = shared_metabolites;
0246 res.shared_reactions    = shared_reactions;
0247 res.mapping_metabolites = mapping_metabolites;
0248 res.mapping_reactions   = mapping_reactions;
0249