embedding_yeast_ccm_make_model

0001 function [model_name, kinetic_models,network,network_CoHid,me_options] = embedding_yeast_ccm_make_model(position_file,network_file, biomass_file, matlab_file)
0002 
0003 eval(default('matlab_file','[]'));
0004 
0005   model_name = 'embedding_yeast_ccm'; % ##
0006 
0007 % Insert yeast ccm model (with kinetics, and without biomass production)
0008 % into yeast ccm + biomass network
0009 
0010 % ---------------------------------------------------------------
0011 % load yeast ccm glycolysis model
0012 % see script yeast_ccm_kinetic_biomass_make_model
0013 
0014 model_name                 = 'yeast_ccm_kinetic';
0015 kk                         = load(network_file);
0016 
0017 kinetic_model_gly          = kk.network;
0018 kinetic_model_gly.s_init   = 0.1*ones(size(kinetic_model_gly.metabolites));
0019 
0020 ind_mets  = label_names({'D-glucose_out','alpha,alpha-trehalose_out','glycerol_out','ATP','pyruvate_out','ethanol_out'},kinetic_model_gly.metabolite_names);
0021 ext_mets = label_names({'ATP','ADP','AMP','NAD+','NADH','UDP','UTP','diphosphate','phosphate'},kinetic_model_gly.metabolite_names);
0022 kinetic_model_gly.s_init(ind_mets)   = [1 2 1 0.5 0.01 0.001]';
0023 kinetic_model_gly.external(ext_mets) = 1;
0024 
0025 kinetic_model_gly_CoHid = kk.network_CoHid;
0026 
0027 
0028 % ---------------------------------------------------------------
0029 % load prepared yeast_ccm plus biomass network
0030 % see script yeast_ccm_kinetic_biomass_make_model.m
0031 
0032 
0033 % ---------------------------------------------------------------
0034 % Kinetic model biomass part
0035 
0036 kk = load(biomass_file);
0037 
0038 kinetic_model_bm = kk.network_biomass;
0039 
0040 ind_rm = label_names({'prod_Ala', 'prod_Arg', 'prod_Asn', 'prod_Asp', 'prod_Cys', 'prod_Glm', 'prod_Gln', 'prod_Gly', 'prod_His', 'prod_Ile', 'prod_Leu', 'prod_Lys', 'prod_Met', 'prod_Orn', 'prod_Phe', 'prod_Pro', 'prod_Ser', 'prod_Thr', 'prod_Trp', 'prod_Tyr', 'prod_Val', 'prod_CTP', 'prod_GTP', 'prod_UTP', 'prod_RIBU5P'},kinetic_model_bm.actions);
0041 kinetic_model_bm          = network_remove_action(kinetic_model_bm,ind_rm);
0042 kinetic_model_bm.external = ones(size(kinetic_model_bm.metabolites));
0043 kinetic_model_bm.external(label_names({'CARBHYD:cyt','RNA:cyt','DNA:cyt','PROT:cyt','LIPID:cyt'},kinetic_model_bm.metabolites)) = 0;
0044 
0045 kinetic_model_bm.kinetics.type              = 'numeric';
0046 kinetic_model_bm.kinetics.velocity_function = @embedding_yeast_ccm_bm_velocities;
0047  %% NOTE THAT INDICES ARE HARD_CODED WITHIN VELOCITY FUNCTION!!!
0048 
0049 kinetic_model_bm.kinetics.parameters = struct('mu',1,'CARBHYD',1,'DNA',1,'LIPID',1,'PROT',1,'RNA',1);
0050 kinetic_model_bm.s_init              = ones(size(kinetic_model_bm.metabolites));
0051 
0052 
0053 % ---------------------------------------------------------------
0054 % Network model
0055 
0056 network          = kk.network_combined;
0057 network_CoHid    = kk.network_combined_CoHid;
0058 network          = netgraph_read_positions(network,position_file,[0,0],1,0);
0059 network_CoHid    = netgraph_read_positions(network_CoHid,position_file,[0,0],1,0);
0060 
0061 
0062 % ----------------------------------------------------
0063 
0064 % netgraph_edit_positions(network_CoHid,position_file,0);
0065 
0066 
0067 % figure(1); netgraph_draw(kinetic_model_gly_CoHid);
0068 
0069 
0070 
0071 % ----------------------------------------------------
0072 % extract only reactions to be modelled kinetically
0073 
0074 % ----------------------------------------------------
0075 % fit concentrations to match approximately the predefined fluxes
0076 
0077 % ind_v_gly = label_names(kinetic_model_gly.metabolites,network.metabolites);
0078 % v_gly = v(ind_v_gly);
0079 % c_gly_init = kinetic_model_gly.kinetics.c;
0080 %
0081 % [c_gly, u_ratio_gly] = network_fit_concentrations_to_rates(network,v_gly);
0082 
0083 % ----------------------------------------------------
0084 % plot fluxes (=> see hepatonet model)
0085 
0086 % ---------------------------------------------------------------
0087 % determine flux distribution on network by FBA
0088 
0089 fba_constraints = fba_default_options(network);
0090 
0091 % no import of some metabolites
0092 no_consumption = {'ethanol_out', 'glycerol_out', 'pyruvate_out', 'alpha,alpha-trehalose_out'};
0093 ext_met        = network.metabolite_names(find(network.external));
0094 ind_no_cons    = label_names(no_consumption,ext_met);
0095 fba_constraints.ext_sign(ind_no_cons) = 1;
0096 
0097 ind_macromolecule_dilution = label_names({'prod_CARBOHYD:dil','prod_DNA:dil','prod_PROT:dil','prod_LIPD:dil','prod_RNA:dil'},network.actions);
0098 
0099 % target: production of macromolecules (equal weight for all five sorts)
0100 fba_constraints.zv(ind_macromolecule_dilution) = 1;
0101 
0102 % no consumption of macromolecules
0103 fba_constraints.v_sign(ind_macromolecule_dilution) = 1;
0104 
0105 reaction_fixed_signs = {'prod_Ala', 'prod_Arg', 'prod_Asn', 'prod_Asp', 'prod_Cys', 'prod_Glm', 'prod_Gln', 'prod_Gly', 'prod_His', 'prod_Ile', 'prod_Leu', 'prod_Lys', 'prod_Met', 'prod_Orn', 'prod_Phe', 'prod_Pro', 'prod_Ser', 'prod_Thr', 'prod_Trp', 'prod_Tyr', 'prod_Val', 'prod_CTP', 'prod_GTP', 'prod_UTP', 'prod_CARBHYD:cyt', 'prod_DNA:cyt', 'prod_LIPID:cyt', 'prod_PROT:cyt', 'prod_RNA:cyt', 'prod_CARBOHYD:dil', 'prod_DNA:dil', 'prod_PROT:dil', 'prod_LIPD:dil', 'prod_RNA:dil', 'prod_RIBU5P'};
0106 fixed_signs     = ones(size( reaction_fixed_signs));
0107 
0108 fba_constraints.v_sign(label_names(reaction_fixed_signs,network.reaction_names)) = fixed_signs;
0109 
0110 v       = fba(network,fba_constraints);
0111 
0112 % figure(2); netgraph_concentrations(network_CoHid,[],v,1,struct('actstyle','none','arrowstyle','fluxes','arrowsize',0.1));
0113 
0114 % figure(3); netgraph_concentrations(network_CoHid,[],10*v,1,struct('actstyle','none'))
0115 
0116 % --------------------------------------------------------------------------------
0117 % ids for mapping the models
0118 % which name fields should be used for mapping the elements?
0119 
0120 id = [];
0121 
0122 % ---------------------------------------------------------------
0123 % save
0124 
0125 kinetic_models = {kinetic_model_gly, kinetic_model_bm};
0126 
0127 me_options.set_internal    = {'RIBU5P'};
0128 me_options.id              = id;
0129 me_options.fba_constraints = fba_constraints;
0130 me_options.enforce_stationarity = 1;
0131 me_options.arrowsize       = 0.02;
0132 me_options.squaresize      = 0.02;
0133 me_options.text_offset     = [0.02,0.01];
0134 me_options.arrowvaluesmax  = 0.05;
0135 me_options.legendFontsize  = 4;
0136 me_options.fontsize  = 6;
0137 me_options.position_file   = position_file;
0138 me_options.cofactors       = {'NADPH','NADP','NH4','s_0764_b', 's_1434_b', 's_0470_b', 's_0766_b', 's_0446', 's_0400', 's_0434', 's_1207', 's_1082', 's_1087', 's_1430', 's_1411', 'CTP', 'O2', 'GTP','SO4'};
0139 me_options.v_network       = v;
0140 
0141 me_options.c_init          = [];
0142 me_options.Tmax            = 20; 
0143 
0144 if length(matlab_file),
0145   save(matlab_file,'kinetic_models','network','network_CoHid','me_options');
0146 end