search

SysBioChalmers / yeast-GEM

The consensus GEM for Saccharomyces cerevisiae
http://sysbiochalmers.github.io/yeast-GEM/
Creative Commons Attribution 4.0 International
95 stars 46 forks source link

Proposed changes to the anaerobic model #352

Open davidhcsic opened 1 year ago

davidhcsic commented 1 year ago

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 commented 1 year ago

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:

davidhcsic commented 1 year ago

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: