From b7e5b77f867a8c23a28e536c269a51d9d266182d Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 21 Aug 2024 21:07:24 +0200 Subject: [PATCH] feat: initial model curations as defined in May --- code/otherChanges/anaerobicModel.m | 113 ++++++++++-- code/otherChanges/test_Sjoberg_20240524.m | 211 ++++++++++++++++++++++ 2 files changed, 308 insertions(+), 16 deletions(-) create mode 100644 code/otherChanges/test_Sjoberg_20240524.m diff --git a/code/otherChanges/anaerobicModel.m b/code/otherChanges/anaerobicModel.m index 3ee6d8b4..79d5f448 100644 --- a/code/otherChanges/anaerobicModel.m +++ b/code/otherChanges/anaerobicModel.m @@ -4,21 +4,27 @@ % % Inputs: model (struct) aerobic model % Output: model (struct) anaerobic model -% +% % Usage: model = anaerobicModel(model) % -%1st change: Refit GAM and NGAM to exp. data, change biomass composition -GAM = 30.49; %Data from Nissen et al. 1997 + +GAM = 55.2; + + P = 0.461; %Data from Nissen et al. 1997 -NGAM = 0; %Refit done in Jouthen et al. 2012 +NGAM = 1; model = changeGAM(model,GAM,NGAM); model = scaleBioMass(model,'protein',P,'carbohydrate',false); -%2nd change: Removes the requirement of heme a, NAD(PH), coenzyme A in the biomass equation -% (not used under anaerobic conditions) -mets = {'s_3714','s_1198','s_1203','s_1207','s_1212','s_0529'}; +% %2nd change: Removes the requirement of heme a, NAD(PH), coenzyme A in the biomass equation +% % (not used under anaerobic conditions) + +% I disagree. Please do not remove NAD(P)H and CoA, include the exhange +% reactions for the relevant vitamins instead. See below /Gustav +mets = {'s_3714'}; %,'s_1198','s_1203','s_1207','s_1212','s_0529'}; + [~,met_index] = ismember(mets,model.mets); model.S(met_index,strcmp(model.rxns,'r_4598')) = 0; @@ -30,16 +36,91 @@ model.lb(strcmp(model.rxns,'r_1994')) = -1000; %palmitoleate model.lb(strcmp(model.rxns,'r_2106')) = -1000; %zymosterol model.lb(strcmp(model.rxns,'r_2134')) = -1000; %14-demethyllanosterol -model.lb(strcmp(model.rxns,'r_2137')) = -1000; %ergosta-5,7,22,24(28)-tetraen-3beta-ol +%remove this due to NADH recycling to ergosterol +model.lb(strcmp(model.rxns,'r_2137')) = 0; %ergosta-5,7,22,24(28)-tetraen-3beta-ol model.lb(strcmp(model.rxns,'r_2189')) = -1000; %oleate -%4th change: Blocked pathways for proper glycerol production -%Block oxaloacetate-malate shuttle (not present in anaerobic conditions) -model.lb(strcmp(model.rxns,'r_0713')) = 0; %Mithocondria -model.lb(strcmp(model.rxns,'r_0714')) = 0; %Cytoplasm -%Block glycerol dehydroginase (only acts in microaerobic conditions) -model.ub(strcmp(model.rxns,'r_0487')) = 0; -%Block 2-oxoglutarate + L-glutamine -> 2 L-glutamate (alternative pathway) -model.ub(strcmp(model.rxns,'r_0472')) = 0; + + +% Added exchange of vitamins enabling NAD(P)H and CoA syntheis in anaerobic +% conditions /Gustav +model.lb(strcmp(model.rxns,'r_1967')) = -1000; %nicotinate +model.lb(strcmp(model.rxns,'r_1548')) = -1000; %(R)-pantothenate + + +%% Changes to the model that give correct phenotype for anaerobic batch growht on minimal glucose media + +% Inhibit MDH2 during excess glucose = anaerobic conditions. +model = setParam(model,'lb','r_0714',0); +model = setParam(model,'ub','r_0714',0); + +%Glutamate synthae repressed in excess nitrogen +model.ub(strcmp(model.rxns,'r_0472'))=0; + +% GCY1 has a positive DeltaG and is part of a transhydrogenase cycle NADH -> NADPH +model.ub(strcmp(model.rxns,'r_0487')) = 0; + +%glycine cleavage respressed by presence of glucose +model.ub(strcmp(model.rxns,'r_0501'))=0; %glycine cleavage, mitochondrion +model.lb(strcmp(model.rxns,'r_0501'))=0; +model.ub(strcmp(model.rxns,'r_0507'))=0; %glycine cleavage complex (lipoylprotein), mitochondrion +model.lb(strcmp(model.rxns,'r_0507'))=0; +model.ub(strcmp(model.rxns,'r_0509'))=0; %glycine cleavage complex (lipoamide), mitochondrion +model.lb(strcmp(model.rxns,'r_0509'))=0; + +% MAE1 and IDP are likely not major mitochondrial NADPH sources. +% model.ub(strcmp(model.rxns,'r_0719')) = 0; % malic enzyme (MAE1), mitochondrion +% model.ub(strcmp(model.rxns,'r_2131')) = 0; % isocitrate dehydrogenase (IDP1), mitochondrion +% %model.ub(strcmp(model.rxns,'r_0659')) = 0; % isocitrate dehydrogenase (IDP2), mitochondrion + +%========================================================================= +%Speculative, this cleans up alternate sources of mitochondrial pyruvate +%(not from the transporter). No effect on glycerol production +%model.ub(strcmp(model.rxns,'r_0718')) = 0; % malic enzyme (MAE1), mitochondrion, NADH reaction, acts as major mitochondrial pyruvate source +%model.ub(strcmp(model.rxns,'r_4701')) = 0; % IRC7, Cysteine desulphydrase, enables growth on cysteine as nitrogen source + + + + +% heme a[c] heme a bigg.metabolite/hemeA;chebi/CHEBI:24479;kegg.compound/C15670;metanetx.chemical/MNXM53309;sbo/SBO:0000247 C49H55FeN4O6 c s_3714 -3% cofactor[c] cofactor sbo/SBO:0000649 c s_4205 #NUM! +% r_4598 cofactor pseudoreaction 0.00019 coenzyme A[c] + 1e-05 FAD[c] + 0.00265 NAD[c] + 0.00015 NADH[c] + 0.00057 NADP(+)[c] + 0.0027 NADPH[c] + 0.00099 riboflavin[c] + 1.2e-06 TDP[c] + 6.34e-05 THF[c] + 1e-06 heme a[c] => cofactor[c] +% model.S(find(strcmp(model.mets,'s_3714')),strcmp(model.rxns,'r_4598'))=0; + +% A reaction converting NADH to NAD + at 3 mmol (g CDW)−1 was coupled to the growth reaction to give the correct ratio of glycerol production to +% glucose consumption, and to align the degree of reduction (as defined by (Heijnen, 1994)) of the modeled biomass to a published value for +% degree of reduction of S. cerevisiae biomass at 4.2 C-mol−1 (Lange and Heijnen, 2001). +% NAD[c] NAD bigg.metabolite/nad;chebi/CHEBI:57540;kegg.compound/C00003;metanetx.chemical/MNXM8;sbo/SBO:0000247 C21H26N7O14P2 c s_1198 -1 +% NADH[c] NADH bigg.metabolite/nadh;chebi/CHEBI:57945;kegg.compound/C00004;metanetx.chemical/MNXM10;sbo/SBO:0000247 C21H27N7O14P2 c s_1203 -2 +% H+[c] H+ bigg.metabolite/h;chebi/CHEBI:24636;kegg.compound/C00080;metanetx.chemical/MNXM1;sbo/SBO:0000247 H c s_0794 1 + +%% 3mmol (g CDW)−1s +RD=3; +% RD=0; + +% Updated the adjustment of degree of reduciton /Gustav + +%Try NADPH for DR balance instead. Also, try to convert NADH to NADPH at an adjustable ratio +% is not required. But we can consider if NADH or NADPH should be used to +% balance the degree of reduction of the biomass +NADH_NADPH=0; + +%% NADPH +model.S(find(strcmp(model.mets,'s_1212')),strcmp(model.rxns,'r_4041'))=model.S(find(strcmp(model.mets,'s_1212')),strcmp(model.rxns,'r_4041'))-RD-NADH_NADPH; +%% NADP + model.S(find(strcmp(model.mets,'s_1207')),strcmp(model.rxns,'r_4041'))=model.S(find(strcmp(model.mets,'s_1207')),strcmp(model.rxns,'r_4041'))+RD+NADH_NADPH; +%% NADH +model.S(find(strcmp(model.mets,'s_1203')),strcmp(model.rxns,'r_4041'))=model.S(find(strcmp(model.mets,'s_1203')),strcmp(model.rxns,'r_4041'))+NADH_NADPH; +%% NAD + model.S(find(strcmp(model.mets,'s_1198')),strcmp(model.rxns,'r_4041'))=model.S(find(strcmp(model.mets,'s_1198')),strcmp(model.rxns,'r_4041'))-NADH_NADPH; +%% H+ + model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4041'))=model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4041'))-RD; + + +%% Adjust Mw of biomass to 1000 + +% Biomass_MW=computeMetFormulae(model,'metMwRange','s_0450','fillMets','none','printLevel',0); +% model.S(:,strcmp(model.rxns,'r_4041')) = model.S(:,strcmp(model.rxns,'r_4041'))*1000/mean(Biomass_MW); +% model.S(find(strcmp(model.mets,'s_0450')),strcmp(model.rxns,'r_4041')) = 1; +% end diff --git a/code/otherChanges/test_Sjoberg_20240524.m b/code/otherChanges/test_Sjoberg_20240524.m new file mode 100644 index 00000000..3666d309 --- /dev/null +++ b/code/otherChanges/test_Sjoberg_20240524.m @@ -0,0 +1,211 @@ +clear; +close all; +load('yeast-GEM.mat') + +% this script makes corrections that should be applied to *all* models, not +% just the anaerobic. Adjustments to the anaerobic model are found in the +% anaerobicModel function + + +%========================================================================= + +% Create a pathway for anaplerosis. Check alternatives? this one might work +% on other diacids also? If considered important for "flux reasons", dig +% into this transporter literature +% 'r_1264' 'succinate transport' 'phosphate[mitochondrion] + succinate[cytoplasm] -> phosphate[cytoplasm] + succinate[mitochondrion] ' +%model.lb(findRxnIDs(model,'r_1264'))=-1000; %DIC1. I am not sure its actually a phosphate/succinate antiport. It not clear from the article. It could also be specific for both phosphate and succinate + +%========================================================================= + + +% look for all proton symport/antiport reactions and make sure that they +% only enter the cell. + +%'s_0794' H+ cytosol +%'s_0796' H+ extracellular +symporterIDs = intersect(find(model.S(find(strcmp(model.mets,'s_0796')),:)),find(model.S(find(strcmp(model.mets,'s_0794')),:)),'stable'); +for i = 1:length(symporterIDs) + if model.S(find(strcmp(model.mets,'s_0796')),symporterIDs(i))<0 + model.lb(symporterIDs(i))=0; + else + model.ub(symporterIDs(i))=0; + end + + if find(strcmp(model.rxns,'r_1258'))==symporterIDs(i) % ignore the sodium transporter, without it, the model doesn't work + model.lb(symporterIDs(i))=-1000; + model.ub(symporterIDs(i))=1000; + end +end + +%========================================================================= + +% This section balances reactions and ensures that a correct molecular +% weight can be calculated for the biomass + +% Set the charge of all biomass components to 0, should be applied to *all* models, not just anaerobic. +model.metCharges(strcmp(model.mets,'s_3717'))=0; % Protein +model.metCharges(strcmp(model.mets,'s_3718'))=0; % Carbohydrate +model.metCharges(strcmp(model.mets,'s_3719'))=0; % RNA +model.metCharges(strcmp(model.mets,'s_3720'))=0; % DNA +model.metCharges(strcmp(model.mets,'s_3746'))=0; % Lipid backbone +model.metCharges(strcmp(model.mets,'s_3747'))=0; % Lipid chain +model.metCharges(strcmp(model.mets,'s_4205'))=0; % Cofactor +model.metCharges(strcmp(model.mets,'s_4206'))=0; % Ion + +% Make the charge of K and Na 1+, should be applied to *all* models, not just anaerobic. +model.metCharges(strcmp(model.mets,'s_1373'))=1; +model.metCharges(strcmp(model.mets,'s_1374'))=1; +model.metCharges(strcmp(model.mets,'s_3776'))=1; +model.metCharges(strcmp(model.mets,'s_1437'))=1; +model.metCharges(strcmp(model.mets,'s_1438'))=1; +model.metCharges(strcmp(model.mets,'s_3775'))=1; + +% Balance the charge of all biomass component pseudo reactions by adding the required amount of H+ +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4047')) = -sum(model.S(:,strcmp(model.rxns,'r_4047')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4049')) = -sum(model.S(:,strcmp(model.rxns,'r_4049')).*model.metCharges,'omitnan'); % RNA +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4050')) = -sum(model.S(:,strcmp(model.rxns,'r_4050')).*model.metCharges,'omitnan'); % DNA +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4598')) = -sum(model.S(:,strcmp(model.rxns,'r_4598')).*model.metCharges,'omitnan'); % Cofactor +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4599')) = -sum(model.S(:,strcmp(model.rxns,'r_4599')).*model.metCharges,'omitnan'); % Ion + + +% Special case for SLIME rxns +model.metCharges(find(contains(model.metNames,'chain')+contains(model.metNames,'backbone'))) = 0; + +% Now, based on the charge balance, find all the reactions that are +% imbalanced, add or remove hydrogen as necessary + + +% Balance the charge of all imbalanced SLIME reactions by adding the required amount of H+, +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_3975')) = -sum(model.S(:,strcmp(model.rxns,'r_3975')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_3976')) = -sum(model.S(:,strcmp(model.rxns,'r_3976')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_3977')) = -sum(model.S(:,strcmp(model.rxns,'r_3977')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_3978')) = -sum(model.S(:,strcmp(model.rxns,'r_3978')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4076')) = -sum(model.S(:,strcmp(model.rxns,'r_4076')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4077')) = -sum(model.S(:,strcmp(model.rxns,'r_4077')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4078')) = -sum(model.S(:,strcmp(model.rxns,'r_4078')).*model.metCharges,'omitnan'); % Protein +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4079')) = -sum(model.S(:,strcmp(model.rxns,'r_4079')).*model.metCharges,'omitnan'); % Protein + + +% Make the charge of HS (hydrogen sulfide) -1 +model.metCharges(strcmp(model.mets,'s_0841'))=-1; +model.metCharges(strcmp(model.mets,'s_3906'))=-1; +model.metCharges(strcmp(model.mets,'s_4263'))=-1; + +% Balance the reaction r_4702, 'L-cysteine:2-oxoglutarate aminotransferase' +% by adding a proton as reactant +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4702'))=-1; + +% Balance the reaction r_4703, 'L-cysteine:2-oxoglutarate aminotransferase' +% by adding a proton as product +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_4703'))=1; + +% Balance the reactions 'r_0774' and 'r_0775', 'NAPRtase' by removing H+ consumption +% and adding a H2O as a reactant +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_0774'))=0; +model.S(find(strcmp(model.mets,'s_0803')),strcmp(model.rxns,'r_0774'))=-1; +model.S(find(strcmp(model.mets,'s_0794')),strcmp(model.rxns,'r_0775'))=0; +model.S(find(strcmp(model.mets,'s_0803')),strcmp(model.rxns,'r_0775'))=-1; + +%========================================================================= +% This section focuses on individual reactions that have the wrong +% reversibility/direction/cofactor or should be completley removed + +% rename r_0227, it is the plasma membrane ATPase, not a cytosolic ATPase +model.rxnNames(strcmp(model.rxns,'r_0227')) = {'ATPase, plasma membrane'}; + +% make sure both formate-THF ligases are reversible. +model.lb(strcmp(model.rxns,'r_0446')) = -1000; + +% ADE3 and MIS1, methylenetetrahydrofolate dehydrogenase (NADP+) [EC 1.5.1.5] +% make irrevrersible +model.ub(strcmp(model.rxns,'r_0732')) = 0; +model.ub(strcmp(model.rxns,'r_0733')) = 0; + +% There is no evidence for this PFK1 side reaction in yeast. Consider +% removing the reaction completley +model.ub(strcmp(model.rxns,'r_0887')) = 0; + +% TYR1 incorrectly annotated as using NAD, should be NADP +model.S(find(strcmp(model.mets,'s_1212')),strcmp(model.rxns,'r_0939'))=0; %NADPH +model.S(find(strcmp(model.mets,'s_1207')),strcmp(model.rxns,'r_0939'))=0; %NADP +model.S(find(strcmp(model.mets,'s_1203')),strcmp(model.rxns,'r_0939'))=1; %NADH +model.S(find(strcmp(model.mets,'s_1198')),strcmp(model.rxns,'r_0939'))=-1;%NAD + +% Make esterification reactions irreversible. positive deltaG +model.ub(strcmp(model.rxns,'r_4713')) = 0; %diethyl succinate +model.ub(strcmp(model.rxns,'r_4714')) = 0; %monoethyl succinate + +% Make polyphosphate hydrolase and diphosphate transport over cell membrane +% both irreversible +model = setParam(model,'lb',{'r_4723','r_4724','r_4725'},0); +model = setParam(model,'lb','r_4460',0); + +%========================================================================= + + +% Convert to anaerobic +model = anaerobicModel(model); + +% Glucose uptake rate +model = setParam(model,'eq','r_1714',-23); + +% Solve +res=solveLP(model,1); +FLUX=res.x; + +v_glc=FLUX(getIndexes(model,'r_1714','rxns'),:); +v_eth=FLUX(getIndexes(model,'r_1761','rxns'),:); +v_CO2=FLUX(getIndexes(model,'r_1672','rxns'),:); +v_gly=FLUX(getIndexes(model,'r_1808','rxns'),:); +v_growth=FLUX(getIndexes(model,'r_4041','rxns'),:); +%% +%glycerol ethanol Co2 +%4.5 ± 0.4 31 ± 2 38 ± 10 +data=[4.5 31 38 0.36]; +sim=[v_gly v_eth v_CO2 v_growth]; +errorVal=[0.4 2 10 0.02]; +b1=bar(data./data,'FaceAlpha',0.5);hold on;b2=bar(sim./data,'FaceAlpha',0.5); + +hold on + +er = errorbar([1 2 3 4],data./data,errorVal./data,errorVal./data); +er.Color = [0 0 0]; +er.LineStyle = 'none'; + + +legend({'data','simulation'}); +ylabel('Relative value'); +xticklabels({'Glycerol','Ethanol','CO2','Biomass'}) + +temp_model=model; +temp_model.metNames=strcat( model.metNames, repmat('[',length(model.mets),1),model.compNames(temp_model.metComps),repmat(']',length(model.mets),1) ); + + +% Pack everything into a nice table +rxns_reacs=printRxnFormula(temp_model,'rxnAbbrList',model.rxns,'metNameFlag',1,'printFlag',0); +%[massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(model); +%tab=table(model.rxns,model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,model.grRules,imBalancedRxnBool,imBalancedCharge,imBalancedMass); +tab=table(model.rxns,model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,model.grRules); +tab = sortrows(tab,"Var4","descend"); + +printFluxes(model,FLUX,true) + + + + + +%% Calculate formula and degree of reduciton of biomass + +[mwRange,metFormulae,elements,metEle]=computeMetFormulae(model,'metMwRange','s_0450','fillMets','none','printLevel',0); +mwRange + +% Print out the results of the +Biomass_index = find(strcmp(model.metNames,'biomass')); + +%Degree of reduction per element. Order of the elements 'C', 'H', 'O', 'N' +DR_per_ele = [4, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]; +DR_per_Cmol=sum(metEle(Biomass_index,:).*DR_per_ele)/metEle(Biomass_index,1) + +Biomass_formula_Cmol=metEle(Biomass_index,:)/metEle(Biomass_index,1); + +MW_per_Cmol_min = mwRange/metEle(Biomass_index,1) \ No newline at end of file