Description of embedding_element

0001 function [mapping_metabolites, mapping_reactions, covered_metabolites, covered_reactions, shared_metabolites, shared_reactions, network, network_CoHid, network_CoSplit] = embedding_element_mapping(kinetic_models,network,me_options)
0002   
0003 % -----------------------------------------------------------------------------------------
0004 % build the mappings for each kinetic model:
0005 %     mapping_metabolites{it}
0006 %     mapping_reactions{it}
0007 %  if elements do not yet exist in the network, create them
0008 %
0009 % Reactions to be mapped must have the same directions!!!
0010 %  (This is checked at the end of this function)
0011 
0012 nm_orig = length(network.metabolites);
0013 
0014 for it = 1:length(kinetic_models),
0015   
0016   %% mapping for metabolites
0017   mapping_metabolites{it} = label_names(kinetic_models{it}.(me_options.id.metabolites_kinetic_models{it}),network.(me_options.id.metabolites_network));
0018   ind_nonmapped           = find(mapping_metabolites{it}==0);
0019   nonmapped_metabolites   = kinetic_models{it}.metabolites(ind_nonmapped);
0020   nonmapped_metabolite_id = kinetic_models{it}.(me_options.id.metabolites_kinetic_models{it})(ind_nonmapped);
0021   
0022   if ind_nonmapped, 
0023     [network,indices] = network_add_metabolites(network,nonmapped_metabolites,kinetic_models{it}.external(ind_nonmapped));
0024     network.(me_options.id.metabolites_network)(indices) = column(nonmapped_metabolite_id);
0025     network.(me_options.id.metabolites_network) = column(network.(me_options.id.metabolites_network));
0026     %network.(me_options.id.metabolites_network)(end-length(ind_nonmapped)+1:end) = nonmapped_metabolite_id;
0027     display('  Adding metabolites (from kinetic model) to network model:'); 
0028     fprintf(mytable(nonmapped_metabolites,1))
0029     %%display('  Adding them to the network. Metabolites in augmented network:');
0030     %%network.metabolites
0031     if ~strcmp(me_options.id.metabolites_network,'metabolites'),
0032       network.(me_options.id.metabolites_network) = [network.(me_options.id.metabolites_network); kinetic_models{it}.(me_options.id.metabolites_kinetic_models{it})(ind_nonmapped)];
0033     end
0034     mapping_metabolites{it} = label_names(kinetic_models{it}.(me_options.id.metabolites_kinetic_models{it}),network.(me_options.id.metabolites_network));  
0035   end
0036   
0037   %% mapping for reactions
0038 
0039   mapping_reactions{it} = label_names(kinetic_models{it}.(me_options.id.reactions_kinetic_models{it}),network.(me_options.id.reactions_network));
0040   ind_nonmapped         = find(mapping_reactions{it}==0);
0041   nonmapped_reactions   = kinetic_models{it}.actions(ind_nonmapped);
0042 
0043   if ind_nonmapped,
0044 
0045     nn = network;
0046 
0047     for itt = 1:length(ind_nonmapped),
0048       n = kinetic_models{it}.N(:,ind_nonmapped(itt));
0049       reactants_kinetic_id = kinetic_models{it}.(me_options.id.metabolites_kinetic_models{it})(find(n));
0050       ll = label_names(reactants_kinetic_id,nn.(me_options.id.metabolites_network));
0051       reactants = nn.metabolites(ll);
0052       stoich_coeffs = n(find(n));
0053       reversible = kinetic_models{it}.reversible(ind_nonmapped(itt));
0054       name       = kinetic_models{it}.actions{ind_nonmapped(itt)};
0055       external   = nn.external(label_names(reactants_kinetic_id,nn.(me_options.id.metabolites_network)));
0056       nn = network_add_reaction(nn,reactants,stoich_coeffs,name,reversible,external);
0057       if ~strcmp(me_options.id.reactions_network,'actions'),
0058         nn.(me_options.id.reactions_network) = [ nn.(me_options.id.reactions_network); kinetic_models{it}.(me_options.id.reactions_kinetic_models{it})(ind_nonmapped(itt))];
0059       end
0060     end
0061 
0062     display('  Adding  reactions (from kinetic model) to network model:'); 
0063     fprintf(mytable(nonmapped_reactions,1))
0064     %%display('  Adding them to the network. Reactions in augmented network:');
0065     %%nn.actions
0066 
0067     display('  Resetting field "kinetics" to cs rate law');
0068     nn.kinetics = set_kinetics(nn,'cs');
0069 
0070     display('  Resetting field "graphics_par"');
0071     nn = netgraph_make_graph(nn);
0072     nn = netgraph_read_positions(nn,me_options.position_file,[0,0],1,0);
0073 
0074     mapping_reactions{it} = label_names(kinetic_models{it}.(me_options.id.reactions_kinetic_models{it}),nn.(me_options.id.reactions_network));
0075 
0076     network       = nn;
0077   end
0078 
0079 end
0080 
0081 network_CoHid = netgraph_simple_graph(network,me_options.cofactors);
0082 network_CoHid = netgraph_read_positions(network_CoHid,me_options.position_file,[0,0],1,0);
0083 network_CoHid.graphics_par.me_options.omitreactions = me_options.omitreactions;
0084 if isfield(network_CoHid,'metabolite_names'),
0085   network_CoHid.graphics_par.metnames = escape_underscores(network_CoHid.metabolite_names);
0086 end
0087 
0088 network_CoSplit = netgraph_split_metabolites(network,me_options.cofactors);
0089 network_CoSplit = netgraph_read_positions(network_CoSplit,me_options.position_file,[0,0],1,0);
0090 network_CoSplit.graphics_par.me_options.omitreactions = me_options.omitreactions;
0091 if isfield(network_CoSplit,'metabolite_names'),
0092   network_CoSplit.graphics_par.metnames = escape_underscores(network_CoSplit.metabolite_names);
0093 end
0094 
0095 
0096 % ---------------------------------------------------------------------
0097 % make vectors showing which network elements are covered by which kinetic models
0098 
0099 covered_metabolites = zeros(size(network.metabolites)); 
0100 shared_metabolites  = zeros(size(network.metabolites)); 
0101 covered_reactions   = zeros(size(network.actions)); 
0102 shared_reactions    = zeros(size(network.actions)); 
0103 
0104 for  it = length(kinetic_models):-1:1,
0105   %% bit vectors for newly covered
0106   d1 = zeros(size(covered_metabolites));
0107   d2 = zeros(size(covered_reactions));
0108   d1(mapping_metabolites{it}) = 1;
0109   d2(mapping_reactions{it})   = 1;
0110   o1 = d1 .* covered_metabolites;
0111   o2 = d2 .* covered_reactions;
0112 
0113   covered_metabolites(find(d1)) = it;
0114   covered_metabolites(find(o1)) = it;
0115   shared_metabolites(find(o1))  = 1;
0116 
0117   covered_reactions(find(d2))   = it;
0118   covered_reactions(find(o2))   = it;
0119   shared_reactions(find(o2))    = 1;
0120 
0121   covered_metabolites(mapping_metabolites{it}) = it;
0122   covered_reactions(mapping_reactions{it})     = it;
0123   
0124 end
0125 
0126 if find(shared_reactions), 
0127   display('  Reactions appear in several kinetic models');
0128   network.actions(find(shared_reactions))
0129 end
0130 
0131 if find(shared_metabolites), 
0132   display('  Metabolites appear in several kinetic models');
0133   network.metabolites(find(shared_metabolites))
0134 end
0135 
0136 
0137 % ---------------------------------------------------------------------
0138 % check if reaction stoichiometries match
0139 
0140 for it = 1:length(kinetic_models),
0141   dN = double(network.N(mapping_metabolites{it}, mapping_reactions{it}) - kinetic_models{it}.N);
0142   if norm(full(dN)) > 0,
0143     ind_conflict = find(sum(abs(dN)));
0144     warning('Conflict between stoichiometries in kinetic models and network');
0145     display(sprintf('Kinetic model number %d',it));
0146     display(sprintf('Reaction %d ',ind_conflict));
0147     formulae_kin = network_print_formulae(kinetic_models{it});
0148     formulae_kin(ind_conflict)
0149     display(sprintf('Updating the global network with stoichiometries from model %d',it));
0150     network.N(mapping_metabolites{it}, mapping_reactions{it}(ind_conflict)) = kinetic_models{it}.N(:,ind_conflict);
0151   else 
0152     display('  Stoichiometries are consistent')
0153   end
0154 end
0155 
0156 
0157 % ---------------------------------------------------------------------
0158 % in principle, metabolites are set internal or external in the combined model
0159 % the same way they were in the network; but all metabolites that are internal
0160 % in one of the kinetic models are also set internal in the combined model
0161 % and all metabolites that did not exist in the network are set like in the kinetic models
0162 
0163 display('  Setting internal metabolites (from kinetic models) internal in network');   
0164 for it = 1:length(kinetic_models),
0165   network.external(mapping_metabolites{it}(find(kinetic_models{it}.external==0))) = 0;
0166   ind_new = find(mapping_metabolites{it}>nm_orig);
0167   network.external(mapping_metabolites{it}(ind_new)) = kinetic_models{it}.external(ind_new);
0168 end
embedding_element_mapping

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
