Proposed changes to the anaerobic model

[davidhcsic]

After inspecting the results obtained with the anaerobic version model and comoparing these with transcriptomics and fluxomics, I propose some changes to improve the accuracy of anaerobic—bioreactor data from Nissen et al. 1997.
nissen.xlsx

I tried constraining all measured extracellular or just glucose along with GAM tunning. The less constrained version is also working consistently good. PPP fluxes are well recovered at lower dilution rates (we have no data at higher rates) and TCA is always consistent. In some cases, the model predicts diphosphate (from phosphate) production which could need fixing.
flux_analysis.xlsx

Fluxomic data
D=0.1h-1. Minimal media. Jouhten, Paula, et al. "Oxygen dependence of metabolic fluxes and energy generation of Saccharomyces cerevisiae CEN. PK113-1A." BMC systems biology 2 (2008): 1-19.
https://pubmed.ncbi.nlm.nih.gov/18613954/

Batch mu=0.23, contains one anaerobic experiment with glucose (comparable with D=0.2). YNB without aminoacids.
Wasylenko, Thomas M., and Gregory Stephanopoulos. "Metabolomic and 13C‐metabolic flux analysis of a xylose‐consuming Saccharomyces cerevisiae strain expressing xylose isomerase." Biotechnology and bioengineering 112.3 (2015): 470-483.
https://pubmed.ncbi.nlm.nih.gov/25311863/

Transcriptomics were taken from:
Tai, S. L., Boer, V. M., Daran-Lapujade, P., Walsh, M. C., de Winde, J. H., Daran, J. M., & Pronk, J. T. (2005).
Two-dimensional transcriptome analysis in chemostat cultures: combinatorial effects of oxygen availability and
macronutrient limitation in Saccharomyces cerevisiae. Journal of Biological Chemistry, 280(1), 437-447.
https://pubmed.ncbi.nlm.nih.gov/15496405/
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE1723
GSE1723_family.txt

function fluxtable=test_proposed_changes_to_anaerobic_model()
clear;
close all;

%% Set this variable to 1 to apply proposed changes in the model
change_model=0;
%% Chose dilution rate from 0.1 (1) to 0.4 (4)
Di=1;
%% Set to 1 to apply constraints on glycerol, ethanol, pyruvate and acetate.
fully_constrained=0;

%% Check GAM implementation
if change_model && fully_constrained
    
    %% Choose GAM that improves PPP for 0.1 and 0.2
    %% 0.3 chose 50 for consistency
    %% 0.4 highest possible
    GAM=[50 50 50 48];
    
elseif change_model && ~fully_constrained
    
    %% Chose GAM that improves ethanol
    GAM=[50 52 50 49];
    
elseif ~change_model && fully_constrained
    
    %% Choose GAM that improves PPP for 0.1 and 0.2
    %% 0.3 and 0.4 highest possible
    GAM=[48 46 42 40];
    
elseif ~change_model && ~fully_constrained
    
    %% Chose GAM that improves ethanol
    GAM=[50 49 48.25 46.5];    
    
end

[data text_pars]=xlsread('nissen.xlsx','data');
variables=text_pars(1:end,1);
Ncarbons=data(:,1);
data(:,1)=[];

for i=1:length(variables)
    eval([variables{i} '=data(i,:)*Ncarbons(i);' ]);
end

load yeast-GEM;

%% Nissen et al. give experimental error for some variables. For the rest we assume as tol=0.01.
tol=0.01;

%% Glucose
model.lb(findRxnIDs(model,'r_1714'))=-glucose_v(Di)*(1+0.006);
model.ub(findRxnIDs(model,'r_1714'))=-glucose_v(Di)*(1-0.006);

%% Growth
model.lb(findRxnIDs(model,'r_4041'))=D(Di)*(1-tol);
model.ub(findRxnIDs(model,'r_4041'))=D(Di)*(1+tol);

if fully_constrained

    %% Ethanol
    model.lb(findRxnIDs(model,'r_1761'))=glucose_v(Di).*(ethanol_y(Di))*(1-tol);
    model.ub(findRxnIDs(model,'r_1761'))=glucose_v(Di).*(ethanol_y(Di))*(1+tol);
    
    %% Carbon dioxide
    model.lb(findRxnIDs(model,'r_1672'))=glucose_v(Di).*CO2_y(Di)*(1-tol);
    model.ub(findRxnIDs(model,'r_1672'))=glucose_v(Di).*CO2_y(Di)*(1+tol);
    
    %% succinate UB
    model.lb(findRxnIDs(model,'r_2056'))=glucose_v(Di).*succinate_y(Di)*(1-tol);
    model.ub(findRxnIDs(model,'r_2056'))=glucose_v(Di).*succinate_y(Di)*(1+tol);
    
    %% Acetic
    model.lb(findRxnIDs(model,'r_1634'))=glucose_v(Di)*acetic_y(Di)*(1-tol);
    model.ub(findRxnIDs(model,'r_1634'))=glucose_v(Di)*acetic_y(Di)*(1+tol);
    
    %% Pyruvate
    model.lb(findRxnIDs(model,'r_2033'))=glucose_v(Di)*pyruvate_y(Di)*(1-tol);
    model.ub(findRxnIDs(model,'r_2033'))=glucose_v(Di)*pyruvate_y(Di)*(1+tol);
    
    %% glycerol exchange
    model.lb(findRxnIDs(model,'r_1808'))=glucose_v(Di)*(glycerol_y(Di))*(1-0.06);
    model.ub(findRxnIDs(model,'r_1808'))=glucose_v(Di)*(glycerol_y(Di))*(1+0.06);
    
end


%% Ammonium
% 1128	'r_1654'	'ammonium exchange'	'ammonium[extracellular]  <=> '	0.144967213
model.lb(findRxnIDs(model,'r_1654'))=-1000;
model.ub(findRxnIDs(model,'r_1654'))=0;

%% Step 1 Original model: Taken from the original function anaerobicModel
model = changeObjective(model, {'r_1654'},1);
% %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;
%% 2-oxoglutarate + L-glutamine -> 2 L-glutamate (alternative pathway)
model.ub(findRxnIDs(model,'r_0472'))=0;

if change_model    
    
    %% Step 2 1796	r_2045_REV'	'serine transport (reversible)'	'L-serine [mitochondrion]  -> L-serine [cytoplasm] '	0.020342928	''
    % Hem25 has been proposed as a transporter of glycine into the mitochondria, where it is used for heme biosynthesis [55].
    % Mitochondrial carriers that transport serine or proline have not been found [98].
    % Mitochondrial Carriers and Substrates Transport Network: A Lesson from Saccharomyces cerevisiae
    % https://www.yeastgenome.org/locus/S000002916
    model.lb(findRxnIDs(model,'r_2045'))=0;
    model.ub(findRxnIDs(model,'r_2045'))=0;
    
    
    %% Step 3  r_0659No1 Isocitrate IDP2 / YLR174W
    %% In step 3 we closed IDP2 in light of the transcription data (3 replicates). IPD2 and MDH2 are known genes for glyoxylate shunt and unlikely to have a
    % major role in anaerobiosis. This rationale can proably applied to most gluconeogenic genes (e.g. ICL, MDH2, etc) and other mestabolic genes 
    % regulated by CAT8 and or ADR1.
    % 10154_at IDP2		
    % Aerobic glucose limited		
    % 661.7	586	759.3
    % Anaerobic carbon limited		
    % 14.8	10.1	10.7
    % FC 44.70945946	58.01980198	70.96261682
    % Log2(FC)	5.482508199	5.858473467	6.148987306
    % Tai, S. L., Boer, V. M., Daran-Lapujade, P., Walsh, M. C., de Winde, J. H., Daran, J. M., & Pronk, J. T. (2005).
    % Two-dimensional transcriptome analysis in chemostat cultures: combinatorial effects of oxygen availability and
    % macronutrient limitation in Saccharomyces cerevisiae. Journal of Biological Chemistry, 280(1), 437-447.
    model.lb(findRxnIDs(model,'r_0659'))=0;
    model.ub(findRxnIDs(model,'r_0659'))=0;
    
    
    %% Block ATP synthase, not likley to operate in anaerobic conditions
    model.lb(findRxnIDs(model,'r_0226'))=0;
    model.ub(findRxnIDs(model,'r_0226'))=0;
    
    
    %% Add reverse ATP synthase. Please check stoichiometry.
    % Acin‐Perez, Rebeca, et al. "Inhibition of ATP synthase reverse activity restores energy homeostasis in mitochondrial pathologies."
    % The EMBO Journal 42.10 (2023): e111699.
    % 's_0397 + 3 s_0794 + s_1326  -> s_0437 + 2 s_0799 + s_0807 '
    % 'ADP[mitochondrion] + 3 H+[cytoplasm] + phosphate[mitochondrion]  -> ATP[mitochondrion] + 2 H+[mitochondrion] + H2O[mitochondrion] '
    % model = addReaction(model,'ATPsynt_rev','metaboliteList',{'s_0397','s_0794','s_1326','s_0437','s_0799','s_0807'},'stoichCoeffList',[1 1 1 -1 -1 -1], 'reversible',false);
    model = addReaction(model,'ATPsynt_rev','metaboliteList',{'s_0397','s_0794','s_1326','s_0437','s_0799','s_0807'},'stoichCoeffList',[1 3 1 -1 -2 -1], 'reversible',false);
    
    
    %% STEP 5 remove glycine production (we might be able to relax this after including other changes).
    %% I commented this. It was a necessary step in the procedure to conclude that we could not block 
    %% glutamate synthase (GOGAT,0472). 
    %% 1254	'r_1810'	'L-glycine exchange'	'L-glycine[extracellular]  -> '	0.047670503	1.063579824	'
    %model.ub(findRxnIDs(model,'r_1810'))=0;
    
    
    %% Step 6, we removed the GOGAT constraint.
    %% 2-oxoglutarate + L-glutamine -> 2 L-glutamate (alternative pathway)
    model.lb(findRxnIDs(model,'r_0472'))=0;
    model.ub(findRxnIDs(model,'r_0472'))=1000;
    
    
    %% Step7: Mitochondrial Malate dehydrogenase should be allowed to operate in reverse. Forward was blocked.
    % 571	'r_0713'	'malate dehydrogenase'	'(S)-malate[mitochondrion] + NAD[mitochondrion]  -> H+[mitochondrion] + NADH[mitochondrion] + oxaloacetate[mitochondrion] '	0.027368804	0.61062536	'YKL085W'
    %https://equilibrator.weizmann.ac.il/reaction?query=NAD++%2B+%28S%29-Malate+%3C%3D%3E+NADH+%2B+Oxaloacetate&p_h=7.5&p_mg=3.0&ionic_strength=0.25&reactantsCoeff=-1.0&reactantsId=11&reactantsName=NAD+&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=1.0&reactantsId=13&reactantsName=NADH&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=1.0&reactantsId=48&reactantsName=Oxaloacetate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=-1.0&reactantsId=92&reactantsName=%28S%29-Malate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&submit=Reverse
    %NADH(aq) + Oxaloacetate(aq) ⇌ NAD (aq) + (S)-Malate(aq)
    %deltaG0=-26.5;
    %deltaG1=-26.5;
    model.lb(findRxnIDs(model,'r_0713'))=-1000;
    model.ub(findRxnIDs(model,'r_0713'))=0;
    
    
    %% Step 8, eliminate carnitine shuttle based in transcriptomics (3 replicates)
    %   225	'r_0254'	'carnitine O-acetyltransferase'	'coenzyme A[mitochondrion] + O-acetylcarnitine[mitochondrion]  -> (R)-carnitine[mitochondrion] + acetyl-CoA[mitochondrion] '	0.035829599	0.799394632	YML042W'
    %   870	'r_1120'	'carnithine-acetylcarnithine carrier'	'(R)-carnitine[mitochondrion] + O-acetylcarnitine[cytoplasm]  -> (R)-carnitine[cytoplasm] + O-acetylcarnitine[mitochondrion] '	0.035829599	0.799394632	YOR100C'
    % 	9668_at	YML042W	
    % 	Aerobic glucose limited		
    % 	578.4	452.6	631.7
    % 	Anaerobic carbon limited		
    % 	72.4	57.2	34.1
    %   FC	7.988950276	7.912587413	18.52492669
    %   Log2(FC)	2.99800595	2.984149532	4.211395928 			
    % 	8437_at	YOR100C	
    % 	Aerobic glucose limited		
    % 	280	214.1	243.4
    % 	Anaerobic carbon limited		
    % 	15.8	6.6	8.6
    %   FC	17.72151899	32.43939394	28.30232558
    %   Log2(FC)	4.147430364	5.019674961	4.822848698
    % Tai, S. L., Boer, V. M., Daran-Lapujade, P., Walsh, M. C., de Winde, J. H., Daran, J. M., & Pronk, J. T. (2005).
    % Two-dimensional transcriptome analysis in chemostat cultures: combinatorial effects of oxygen availability and
    % macronutrient limitation in Saccharomyces cerevisiae. Journal of Biological Chemistry, 280(1), 437-447.
    %% Other arguments for abolishing carnitine shuttle can also be found in:
    % "The transport of Pyr and Oaa across the mitochondrial membrane were included in the model but the
    % transport of AcCoA, the final step of the cytosolic PDH bypass, was omitted since exogenous carnitine would be
    % required for the carnitine shuttle to be active [70-72], and carnitine was not provided in the medium. In addition
    % carnitine acetyltransferase activity has not been detected in S. cerevisiae grown in anaerobic chemostats at 0.1 h-1
    % [35]" Jouhten, Paula, et al. "Oxygen dependence of metabolic fluxes and energy generation of Saccharomyces
    % cerevisiae CEN. PK113-1A." BMC systems biology 2 (2008): 1-19.
    %% AND
    % "since measurements of the activity of carnitine acetyltransferase carried out in this study showed no activity
    % under the present growth conditions (see Results) only the reaction catalysed by acetyl-CoA synthetase has been
    % included in the model (reaction 21)." Nissen, Torben L., et al. "Flux distributions in anaerobic, glucose-limited 
    % continuous cultures of Saccharomyces cerevisiae." Microbiology 143.1 (1997): 203-218.
    %" Export of acetyl CoA is also accomplished through its conjugation to carnitine by the carnitine acetyl transferase Cat2p,
    % which is localized in both peroxisomes and mitochondria. This pathway is only possible when yeast cells are grown in rich
    %media that contain carnitine, which otherwise cannot be synthesized by
    %the yeast cell." Compagno, C., Dashko, S., & Piškur, J. (2014). Introduction to carbon metabolism in yeast.
    % In Molecular mechanisms in yeast carbon metabolism (pp. 1-19). Berlin, Heidelberg: Springer Berlin Heidelberg.
    model.lb(findRxnIDs(model,'r_0254'))=0;
    model.ub(findRxnIDs(model,'r_0254'))=0;
    model.lb(findRxnIDs(model,'r_1120'))=0;
    model.ub(findRxnIDs(model,'r_1120'))=0;
    
    
    %% Add reversible exchange of malate by phosphate. This reversibility exists on IMM904
    %% http://bigg.ucsd.edu/universal/reactions/MALtm
    % Identification of a Novel Gene Encoding the Yeast Mitochondrial Dicarboxylate Transport Protein via Overexpression,
    % Purification, and Characterization of Its Protein Product
    % "First, the reconstituted DTP maintains a strict requirement for intraliposomal counteranion and thus catalyzes an obligatory
    % exchange reaction (Table I). Malonate, malate, succinate, and phosphate support this exchange, whereas α-ketoglutarate, citrate,
    % isocitrate, phosphoenolpyruvate, pyruvate, and ADP do not. These characteristics are identical to those observed with the native
    % transporter following functional reconstitution of the purified protein (20, 21, 23). Second, the overexpressed DTP displays the
    % same external substrate specificity profile as the internal profile described above, a property that is also characteristic of
    % the native carrier in isolated mitochondria (1, 9) and the purified transporter (19, 20, 22, 23).
    %% DSH: My interpretation of this is that the transporter is reversible.
    % 959	'r_1226'	'malate transport'	'(S)-malate[cytoplasm] + phosphate[mitochondrion]  -> (S)-malate[mitochondrion] + phosphate[cytoplasm] '	0	0	'YLR348C'
    model.lb(findRxnIDs(model,'r_1226'))=-1000;
    
    
    %% Step 10 This reversibility is blocked in the cytoplasm.
    %% 343	r_0447'	formate-tetrahydrofolate ligase'	ATP[mitochondrion] + formate[mitochondrion] + THF[mitochondrion]  <=> 10-formyl-THF[mitochondrion] + ADP[mitochondrion] + phosphate[mitochondrion] '	0.023386519	-0.521778214	'YBR084W'
    %deltaG1=2.1
    %deltaG0=2.1
    %% Seems somewhat reversible. I think 
    % Still think it can be reversed for coherence with the cytplasmatic enzyme.
    %https://equilibrator.weizmann.ac.il/reaction?query=ADP+%2B+Orthophosphate+%2B+10-Formyltetrahydrofolate+%3C%3D%3E+ATP+%2B+Formate+%2B+Tetrahydrofolate&p_h=7.5&p_mg=3.0&ionic_strength=0.25&reactantsCoeff=1.0&reactantsId=6&reactantsName=ATP&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=-1.0&reactantsId=10&reactantsName=ADP&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=-1.0&reactantsId=12&reactantsName=Orthophosphate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=1.0&reactantsId=41&reactantsName=Formate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=1.0&reactantsId=74383&reactantsName=Tetrahydrofolate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&reactantsCoeff=-1.0&reactantsId=224&reactantsName=10-Formyltetrahydrofolate&reactantsAbundance=1.000&reactantsAbundanceUnit=millimolar&reactantsPhase=aqueous&submit=Reverse
    %'r_0447'	'formate-tetrahydrofolate ligase'	ATP[mitochondrion] + formate[mitochondrion] + THF[mitochondrion]  <=> 10-formyl-THF[mitochondrion] + ADP[mitochondrion] + phosphate[mitochondrion] '	0.023386519	-0.521778214	'YBR084W'
    %'r_0447'	'formate-tetrahydrofolate ligase'	ATP[mitochondrion] + formate[mitochondrion] + THF[mitochondrion]  <=> 10-formyl-THF[mitochondrion] + ADP[mitochondrion] + phosphate[mitochondrion] '	0.023386519	-0.521778214	'YBR084W'
    model.lb(findRxnIDs(model,'r_0447'))=0;
    
    
    %% Please check diphosphate production. 
    %% 3887	'r_4527'	'diphosphate exchange'	'diphosphate[extracellular]  -> '	0.149172412	3.328190324
      
    model = changeObjective(model, {'r_1654'},1);

    
end

%% Change GAM, add anaerobuc lipid consumption and set biomass changes for protein, RNA and carbohydrates.
model=set_anaerobic_and_biomass(model,protein(Di),RNA(Di),GAM(Di),mannan(Di),trehalose(Di),glycogen(Di),other_cs(Di));
flux=solveLP(model,1);
flux=flux.x;

%% Using the fluxomics data-sets as a reference, the TCA is consistent at all growth rates. PPP is also looking good.

%% Sometimes there is diphosphate production.
% 3820	'r_4460'	'Pyrophosphate transport in via proton symport'	2 H+[extracellular] + diphosphate[extracellular]  <=> diphosphate[cytoplasm] + 2 H+[cytoplasm] '
% 3887	'r_4527'	'diphosphate exchange'	'diphosphate[extracellular]  -> '

model.metNames=strcat(model.metNames,repmat('[',length(model.mets),1),model.compNames(model.metComps),repmat(']',length(model.mets),1));
v_glc=flux(findRxnIDs(model,'r_1714'));
flux_r=abs(flux./v_glc);
fluxtable=table((1:length(model.rxns))',model.rxns,model.rxnNames,printRxnFormula(model,'metNameFlag',true,'rxnAbbrList',model.rxns,'printFlag',0),flux_r,flux,model.grRules);
fluxtable = sortrows(fluxtable,'Var5','descend');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function model = set_anaerobic_and_biomass(model,P,R,GAM,mannan,trehalose,glycogen,other_cs)


%% Set anaerobic factors and blocked reactions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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'};
[~,met_index] = ismember(mets,model.mets);
model.S(met_index,strcmp(model.rxns,'r_4598')) = 0;

%3rd change: Changes media to anaerobic (no O2 uptake and allows sterol
%            and fatty acid exchanges)
model.lb(strcmp(model.rxns,'r_1992')) = 0;        %O2
model.lb(strcmp(model.rxns,'r_1757')) = -1000;    %ergosterol
model.lb(strcmp(model.rxns,'r_1915')) = -1000;    %lanosterol
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
model.lb(strcmp(model.rxns,'r_2189')) = -1000;    %oleate

%% Get current contents and calculate conversion factors for proteins and carbs:
[Pbase,Cbase,Rbase]=calculateContent(model);
Pfactor = P/Pbase;
Rfactor = R/Rbase;
fullGAM = GAM;

%% Change biomass composition:
protPos = strcmp(model.metNames,'protein');
carbPos = strcmp(model.metNames,'carbohydrate');
rnaPos = strcmp(model.metNames, 'RNA');

bioRxn  = model.S(protPos,:) == -1;
protRxn = model.S(protPos,:) == 1;
carbRxn = model.S(carbPos,:) == 1;
rnaRxn =  model.S(rnaPos,:) == 1;

bioPos = strcmp(model.rxnNames,'biomass pseudoreaction');
for i = 1:length(model.mets)
    S_ix  = model.S(i,bioPos);
    isGAM = sum(strcmp({'ATP','ADP','H2O','H+','phosphate'},model.metNames{i})) == 1;
    if S_ix ~= 0 && isGAM
        model.S(i,bioPos) = sign(S_ix)*fullGAM;
    end
end

for i = 1:length(model.mets)
    
    Sbio  = model.S(i,bioRxn);
    Sprot = model.S(i,protRxn);
    Scarb = model.S(i,carbRxn);
    Srna = model.S(i,rnaRxn);
    
    if Sprot ~= 0 && Sprot ~= 1
        %Variable aa content in biomass eq:
        model.S(i,protRxn) = round(Sprot*Pfactor,4);
        
        
        %% Adjusted carbohydrates in a case by case manner
        % Name	ID	%	MW	Si
        % (1->3)-beta-D-glucan	s_0001	0.135	180.16	0.749
        % (1->6)-beta-D-glucan	s_0004	0.045	180.16	0.25
        % chitin	s_0509	0	221.21	0
        % glycogen	s_0773	0.065	180.16	0.361
        % mannan	s_1107	0.128	180.16	0.711
        % trehalose	s_1520	0.047	342.296	0.138
        
        %% other_cs is assumed to be the sum  of (1->3)-beta-D-glucan and (1->6)-beta-D-glucan
        % and is scaled accordinglly. Ignored chitin
    elseif Scarb ~= 0 && Scarb ~= 1
        if (strcmp(model.mets{i},'s_0001'))
            model.S(i,carbRxn) = model.S(i,carbRxn)*other_cs/(0.135+0.045);
        elseif (strcmp(model.mets{i},'s_0004'))
            model.S(i,carbRxn) = model.S(i,carbRxn)*other_cs/(0.135+0.045);
        elseif(strcmp(model.mets{i},'s_0509'))
            %% leave uspared
        elseif(strcmp(model.mets{i},'s_0773'))
            model.S(i,carbRxn) = -glycogen/180.16*1000;
        elseif(strcmp(model.mets{i},'1107'))
            model.S(i,carbRxn) = -mannan/180.16*1000;
        elseif(strcmp(model.mets{i},'s_1520'))
            model.S(i,carbRxn) = -trehalose/342.296*1000;
        end
        
    elseif Srna ~= 0 && Srna ~= 1
        model.S(i,rnaRxn) = round(Srna*Rfactor,4);
    end
    
    
end

end

%% This was taken from an earlier version of the model and I have manually
% removed the 18 g/mol...
function [Ptot,Ctot,Rtot] = calculateContent(model)


% MW aminoacids [g/mol]:
aas = {'s_0404'	89.09-18       % A     Alanine         ala
    's_0542'	121.16-18      % C     Cysteine        cys
    's_0432'	133.11-18      % D     Aspartic acid   asp
    's_0748'	147.13-18      % E     Glutamic acid   glu
    's_1314'	165.19-18      % F     Phenylalanine   phe
    's_0757'	75.07-18       % G     Glycine         gly
    's_0832'	155.15-18      % H     Histidine       his
    's_0847'	131.17-18      % I     Isoleucine      ile
    's_1099'	146.19-18      % K     Lysine          lys
    's_1077'	131.17-18      % L     Leucine         leu
    's_1148'	149.21-18      % M     Methionine      met
    's_0430'	132.12-18      % N     Asparagine      asn
    's_1379'	115.13-18      % P     Proline         pro
    's_0747'	146.14-18      % Q     Glutamine       gln
    's_0428'	174.2-18       % R     Arginine        arg
    's_1428'	105.09-18      % S     Serine          ser
    's_1491'	119.12-18      % T     Threonine       thr
    's_1561'	117.15-18      % V     Valine          val
    's_1527'	204.23-18      % W     Tryptophan      trp
    's_1533'	181.19-18};    % Y     Tyrosine        tyr

% MW carbohidrates [g/mol]:
carbs = {'s_0001'	180.16-18       % (1->3)-beta-D-glucan
    's_0004'	180.16-18       % (1->6)-beta-D-glucan
    's_0509'	221.21-18       % chitin
    's_0773'	180.16-18       % glycogen
    's_1107'	180.16-18       % mannan
    's_1520'	342.296-18 };	% trehalose


% CMP [cytoplasm]	s_0526 323.1965
% GMP [cytoplasm]	s_0782 363.22
% UMP [cytoplasm]	s_1545 324.1813

rnas={'s_0423' 347.2212-18   %AMP
    's_0526' 323.1965-18   %CMP
    's_0782' 363.22-18     %GMP
    's_1545' 324.1813-18 };  %UMP}

%Initialize protein and carb content:
Ptot = 0;
Ctot = 0;
Rtot=0;

%Count protein/carb content in the corresponding pseudo-rxn:
protPos = strcmp(model.metNames,'protein');
carbPos = strcmp(model.metNames,'carbohydrate');
rnaPos = strcmp(model.metNames,'RNA');
protRxn = model.S(protPos,:) == 1;
carbRxn = model.S(carbPos,:) == 1;
rnaRxn = model.S(rnaPos,:) == 1;

for i = 1:length(model.mets)
    
    posP = strcmp(aas(:,1),model.mets{i});
    posC = strcmp(carbs(:,1),model.mets{i});
    posR = strcmp(rnas(:,1),model.mets{i});
    
    if sum(posP) == 1
        Sprot = abs(model.S(i,protRxn));            % mmol/gDW
        Ptot  = Ptot + Sprot*aas{posP,2}/1000;      % mmol/gDW * g/mmol = g/gDW
    elseif sum(posC) == 1
        Scarb = abs(model.S(i,carbRxn));            % mmol/gDW
        Ctot  = Ctot + Scarb*carbs{posC,2}/1000;    % mmol/gDW * g/mmol = g/gDW
    elseif sum(posR) == 1
        Srna = abs(model.S(i,rnaRxn));            % mmol/gDW
        Rtot  = Rtot + Srna*rnas{posR,2}/1000;    % mmol/gDW * g/mmol = g/gDW
    end
end

end

[edkerk]

Thanks, we'll have to go through this to see if the changes affect the anaerobic simulations as well.

To simplify this, could you please clean up the code a little by editing your post and:

• Remove any lines that are not necessary to run (for instance, all the gene_name= commands were probably just used by you when investigating the changes, but are not at all required to run. But don't remove any comments that are useful to understand your code.
• Provide the nissen.xlsx file that is loaded in the beginning.
• Notice that you first allow for reverse flux through r_0714, but later you run set_anaerobic_and_biomass which will again block this reaction.
• In our call you mentioned that you had arguments why certain curations were done, for instance due to low gene expression of certain enzymes. Can you provide such reasoning in the comments, including reference and data?
• With the final model, what happens if you relax the constraints on CO2, succinate, pyruvate, and only constrain glucose and oxygen to Nissen et al., do you get a reasonable flux distribution?
• I don't understand what is the purpose of adding another reaction as objective function?
• What is meant with some protonation issue relating phosphate?

[davidhcsic]

Thanks, we'll have to go through this to see if the changes affect the anaerobic simulations as well.

Thank you. I went through your comments and tried to polish things a bit.

To simplify this, could you please clean up the code a little by editing your post and:

• Remove any lines that are not necessary to run (for instance, all the gene_name= commands were probably just used by you when investigating the changes, but are not at all required to run. But don't remove any comments that are useful to understand your code.
• Provide the nissen.xlsx file that is loaded in the beginning.

Done.

• Notice that you first allow for reverse flux through r_0714, but later you run set_anaerobic_and_biomass which will again block this reaction.

Ups. I had missed that. Updated.

• In our call you mentioned that you had arguments why certain curations were done, for instance due to low gene expression of certain enzymes. Can you provide such reasoning in the comments, including reference and data?

Improved this. Let me know if more detail is needed.

• With the final model, what happens if you relax the constraints on CO2, succinate, pyruvate, and only constrain glucose and oxygen to Nissen et al., do you get a reasonable flux distribution?

Still looks quite reasonable. I have uploaded an excel with the solutions. FC 1 stands for fully constrained and MC1 for model changes.

• I don't understand what is the purpose of adding another reaction as objective function?

The goal was to preserve the idea that glutamate dehydrogensase should be (?) the core route for ammonia assmilation. I removed the secondary objective and I the benefit of using these changes still holds.

• What is meant with some protonation issue relating phosphate?

There is some diphosphate production from phosphate. My first guess is that it is realated with a proton somewhere.