how to generate model for anaerobic conditions

[jiuia]

Hello everyone! Now I'm confused because I don't understand how to use Y8GEM to calculate FBA under anaerobic conditions. Yes, I ran this code “/otherChanges/anaerobicModel” and got the anaerobic model(Anaerobic_y8.mat). The problem now is that when I run this code("Anaerobic_y8.mat"), the program keeps reporting errors.

like that <
model1=load('Anaerobic_y8.mat');
model1 = changeRxnBounds(model1,'r_1714',6.240241765,'u');
model1 = changeObjective(model1,'r_2111',1);
sol = optimizeCbModel(model1,'max');
printFluxes(model1,sol.x)

errors:
errors1 "Reference to a non-existent field 'rxns'.
Error setParam (line 49):index=find(strcmp(rxnList{i},model.rxns),1);

errors2 changeRxnBounds (line 33)
rxnID = findRxnIDs(model,rxnNameList);

errors3 (line 2)
model1 = changeRxnBounds(model1,'r_1714',6.240241765,'u');

Does anyone know how to solve this problem please?

[edkerk]

If you run the following MATLAB code for the yeast-GEM repository root, this should give a functional model:

initCobraToolbox(false)
cd ./code
model = loadYeastModel(); % loading
cd otherChanges\
modelAnaerobic = anaerobicModel(model);
optimizeCbModel(modelAnaerobic)
model1 = changeRxnBounds(modelAnaerobic,'r_1714',6.240241765,'u');
model1 = changeObjective(model1,'r_2111',1);
sol=optimizeCbModel(model1,'max');
printFluxes(model1,sol.x)

However, note that you're using a mixture of RAVEN and COBRA functions (printFluxes is from RAVEN, changeRxnBounds, changeObjective and optimizeCbModel are from COBRA). This is a potential source of conflict (although it works out fine here). Generally, if you load a model with COBRA, you should use COBRA functions, and vice versa for RAVEN. nb. loadYeastModel uses COBRA functions, so you're best off using COBRA functions after using loadYeastModel. If you want to use RAVEN functions instead, you should run model=ravenCobraWrapper(model) to convert the model to RAVEN format.

[edkerk]

Please report the full commands that you run, to be able to replicate what you're doing. At the moment I'm not sure what is in Anaerobic_Y8GEM.mat.

But check what model3 contains after you loaded the mat file, it might contain a nested model structure. For instance model3.model. In that case, you should run model3 = model3.model;. After that, model3.S should be an existing field.

[jiuia]

sorry I missed your message for a long time！

Here is the full code I ran.
`model3 = load('Anaerobic_Y8GEM.mat'); % loading Anaerobic model
model2=ravenCobraWrapper(model3) %Transformation of the model from COBRA to RAVEN
modelAnaerobic = model2;
optimizeCbModel(modelAnaerobic);

model=setParam(modelAnaerobic,'ub',{ 'r_1714' 'r_1762' 'r_1171' 'r_1697'},[6.240241765 8.680804536 1.129925807 11.08604813]); % glucose ethanol glycerol co2
model=setParam(modelAnaerobic,'obj',{'r_2111'},1);
sol=solveLP(model);
printFluxes(model,sol.x,false);`

[jiuia]

Could you please show me a code that can run under anaerobic conditions? The code works in the environment of RAVEN in matlab. please,I really don't see what's wrong! @edkerk

[feiranl]

Hi, we don't have any model in the repo named 'Anaerobic_Y8GEM.mat', thus, we cannot know where is the error. But I do see an issue in your code. if you want to constrain the glucose uptake, it should be negative value and constrains the lower bound rather than the upper bound. Please check the code chunk below. Either of them can give you the functional constrained anaerobic model. The second constrains the glucose uptake, ethanol glycerol co2 production rate upper bound, but the first code chunk only constraints the glucose uptake.

initCobraToolbox(false)
cd ./code
model = loadYeastModel(); % loading
cd otherChanges
model_ana = anaerobicModel(model);
optimizeCbModel(model_ana)
model_ana = changeRxnBounds(model_ana,'r_1714',-6.240241765,'l'); % constrains glucose uptake 
model_ana = changeObjective(model_ana,'r_2111',1);
sol=optimizeCbModel(model_ana,'max');
printFluxes(model_ana,sol.x)

initCobraToolbox(false)
cd model
load('yeast-GEM.mat')
cd ../code/otherChanges
model_ana = anaerobicModel(model);
model_ana=setParam(model_ana,'ub',{'r_1762' 'r_1171','r_1697'},[ 8.680804536 1.129925807 11.08604813]);  %  ethanol glycerol co2
model_ana=setParam(model_ana,'lb',{ 'r_1714'},-6.240241765); % glucose
sol = optimizeCbModel(model_ana);
sol.f

[jiuia]

@feiranl Thank you so much for your advice, it's been particularly helpful and has helped me to see where I'm going wrong!
For your two suggestions, the first form still reports an error when run! But the second form runs with a successful result. However, I still have a question for you in there. I have calculated the metabolic flow separately using data from three different specific growth rates in anaerobic constant culture conditions, e.g. µ=0.1, 0.2, 0.3. The results of the model are calculated as follows, e.g. growth= 0.1389, 0.2790, 0.3325. We can see that the experimental values and the model simulations are very different. Is this normal, please? Or do you think I've done something wrong there!
Looking forward to your reply!

Here is the code I ran.
initCobraToolbox(false) cd model load('yeast-GEM.mat') cd ../code/otherChanges model_ana = anaerobicModel(model); model_ana=setParam(model_ana,'ub',{'r_1762' 'r_1171','r_1697'},[ 8.680804536 1.129925807 11.08604813]); % ethanol glycerol co2 model_ana=setParam(model_ana,'lb',{ 'r_1714'},-6.240241765); % glucose model=setParam(model_ana,'obj',{'r_2111'},1); sol = optimizeCbModel(model_ana,'max'); printFluxes(model_ana,sol.x)

[feiranl]

Hi, @jiuia
We noticed that the growth on carbon limited anaerobic condition is a little bit higher compared with experimental report. But here is the suggestion to tune the GAM of the model (currently in the model is 30.49 for anaerobic condition, but in order to have a lower predicted growth rate, you can increase the GAM ). See the below code for tuning the GAM and NGAM.

GAM = 50;
NGAM = 0;

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

% change NGAM
pos = strcmp(model.rxnNames,'non-growth associated maintenance reaction');%NGAM
model.lb(pos) = NGAM;
model.ub(pos) = NGAM;

• Growth simulation of Yeast8.5
  

• Experimental data for growth on glucose anaerobic condition used in the simulation, you can check wether the data here is deviating a lot from your data.

glucose uptake rate[mmol/gDW/h]	growth rate[1/h]
2.3000000	0.03
2.8599999	0.047
5.8800001	0.1
11.310000	0.19
16.750000	0.28
22.180000	0.369

[jiuia]

@feiranl sorry, I replied to you a few days ago, but i don't know what happened. It didn't succeed! this is the appendix i find.
appendix.pdf

[feiranl]

Hi @jiuia, The appendix you suggested is from GECKO paper, in which the GAM and NGAM is fitted for the enzyme constrained model rather than the basic GEM here we are talking about.

The previous model prediction is higher than the experimental observation, which means that either using the previous GAM or the biomass composition can not represent the cell condition well. You can improve that by either fitting the GAM or adjusting the biomass composition. I suggested 50 for GAM because I found using that value, the mean error of the prediction is lower compared with the experimental observations.

[jiuia]

@feiranl Ok, thank you so much! Your ideal is very useful to me, and i will try to analyze my data.

[jiuia]

@feiranl Hello, I have a question for you, when you run this model under anaerobic conditions, is there a flux change of material in the EMP pathway? And, the flux changes of these substances are not in my run results!
Like this :

EXCHANGE FLUXES:
r_1549	((R,R)-2,3-butanediol exchange):	0.69804
r_1654	(ammonium exchange):	-1.2788
r_1672	(carbon dioxide exchange):	20.2284
r_1714	(D-glucose exchange):	-12.4821
r_1757	(ergosterol exchange):	-0.0075266
r_1761	(ethanol exchange):	17.495
r_1808	(glycerol exchange):	3.0734
r_1814	(glycolaldehyde exchange):	1.3024e-05
r_1832	(H+ exchange):	-0.054087
r_1861	(iron(2+) exchange):	-6.1287e-06
r_1994	(palmitoleate exchange):	-0.018804
r_2005	(phosphate exchange):	-1.016
r_2020	(potassium exchange):	-0.00073182
r_2049	(sodium exchange):	-0.00080036
r_2056	(succinate exchange):	0.021626
r_2060	(sulphate exchange):	-0.01733
r_2100	(water exchange):	2.9018
r_2111	(growth):	0.2016
r_2189	(oleate exchange):	-0.0062141
r_4527	(diphosphate exchange):	0.48147
r_4593	(chloride exchange):	-0.00026007
r_4594	(Cu2(+) exchange):	-0.00013286
r_4595	(Mn(2+) exchange):	-0.00055038
r_4596	(Zn(2+) exchange):	-0.0001508
r_4597	(Mg(2+) exchange):	-0.0002505
r_4600	(Ca(2+) exchange):	-4.3748e-05

[feiranl]

@jiuia yes, the EMP pathway carries fluxes. The printed fluxes is only for exchange reactions. you can check the fluxes for all reactions in sol.x . It should be a matrix with the same length of reactions. Below are codes for checking some reactions in EMP pathway.

rxns = {'r_0534','r_0467','r_0886','r_0892','r_0450','r_1054','r_0892'}
printRxnFormula(model,'rxnAbbrList',rxns,'metNameFlag',true);
[~,idx] = ismember(rxns,model.rxns);
sol.x(idx)

Printed result:

r_0534	ATP [cytoplasm] + D-glucose [cytoplasm] 	->	ADP [cytoplasm] + D-glucose 6-phosphate [cytoplasm] + H+ [cytoplasm] 
r_0467	D-glucose 6-phosphate [cytoplasm] 	<=>	D-fructose 6-phosphate [cytoplasm] 
r_0886	ATP [cytoplasm] + D-fructose 6-phosphate [cytoplasm] 	->	ADP [cytoplasm] + D-fructose 1,6-bisphosphate [cytoplasm] + H+ [cytoplasm] 
r_0892	1,3-bisphospho-D-glycerate [cytoplasm] + ADP [cytoplasm] 	<=>	3-phosphonato-D-glycerate(3-) [cytoplasm] + ATP [cytoplasm] 
r_0450	D-fructose 1,6-bisphosphate [cytoplasm] 	<=>	dihydroxyacetone phosphate [cytoplasm] + glyceraldehyde 3-phosphate [cytoplasm] 
r_1054	dihydroxyacetone phosphate [cytoplasm] 	<=>	glyceraldehyde 3-phosphate [cytoplasm] 
r_0892	1,3-bisphospho-D-glycerate [cytoplasm] + ADP [cytoplasm] 	<=>	3-phosphonato-D-glycerate(3-) [cytoplasm] + ATP [cytoplasm] 

ans =

    1.0000
    0.8704
    0.5390
    1.6717
    0.5390
    0.7388
    1.6717