Skip to content

Commit 59e442a

Browse files
edkerkedkerk
committed
feat: 9.1.0 curate r_0013 and biomass
1 parent 4794a50 commit 59e442a

File tree

4 files changed

+109
-105
lines changed

4 files changed

+109
-105
lines changed

‎code/modelCuration/v9_1_0.m‎

Lines changed: 77 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -58,21 +58,21 @@
5858
%% ========================================================================
5959
% Look for all proton symport/antiport reactions and make sure that they
6060
% only enter the cell.
61-
HcytIdx = getIndexes(model,'s_0794','mets'); % H+[c]
62-
HextIdx = getIndexes(model,'s_0796','mets'); % H+[e]
63-
64-
symporterIDs = find(model.S(HcytIdx,:) & model.S(HextIdx,:));
65-
for i = 1:length(symporterIDs)
66-
if strcmp(model.rxns(symporterIDs(i)), 'r_1258')
67-
% Ignore the sodium transporter, without it, the model does not work
68-
continue
69-
end
70-
if model.S(HextIdx,symporterIDs(i))<0
71-
model.lb(symporterIDs(i))=0;
72-
else
73-
model.ub(symporterIDs(i))=0;
74-
end
75-
end
61+
% HcytIdx = getIndexes(model,'s_0794','mets'); % H+[c]
62+
% HextIdx = getIndexes(model,'s_0796','mets'); % H+[e]
63+
%
64+
% symporterIDs = find(model.S(HcytIdx,:) & model.S(HextIdx,:));
65+
% for i = 1:length(symporterIDs)
66+
% if strcmp(model.rxns(symporterIDs(i)), 'r_1258')
67+
% % Ignore the sodium transporter, without it, the model does not work
68+
% continue
69+
% end
70+
% if model.S(HextIdx,symporterIDs(i))<0 % If defined H+[e] => H+[c]
71+
% model.lb(symporterIDs(i))=0;
72+
% else % If defined H+[c] => H+[e]
73+
% model.ub(symporterIDs(i))=0;
74+
% end
75+
% end
7676

7777
%% ========================================================================
7878
% This section balances reactions and ensures that a correct molecular
@@ -97,32 +97,11 @@
9797
model.metCharges(strcmp(model.mets,'s_3775'))=1;
9898

9999
% Balance the charge of all biomass component pseudo reactions by adding the required amount of H+
100-
Hidx = find(strcmp(model.mets,'s_0794'));
101-
% Protein
102-
rxnIdx = strcmp(model.rxns,'r_4047');
103-
metsStoch = model.S(:,rxnIdx);
104-
metsStoch(Hidx) = 0;
105-
model.S(Hidx,rxnIdx) = -sum(metsStoch.*model.metCharges,'omitnan');
106-
% RNA
107-
rxnIdx = strcmp(model.rxns,'r_4049');
108-
metsStoch = model.S(:,rxnIdx);
109-
metsStoch(Hidx) = 0;
110-
model.S(Hidx,rxnIdx) = -sum(metsStoch.*model.metCharges,'omitnan');
111-
% DNA
112-
rxnIdx = strcmp(model.rxns,'r_4050');
113-
metsStoch = model.S(:,rxnIdx);
114-
metsStoch(Hidx) = 0;
115-
model.S(Hidx,rxnIdx) = -sum(metsStoch.*model.metCharges,'omitnan');
116-
% Cofactor
117-
rxnIdx = strcmp(model.rxns,'r_4598');
118-
metsStoch = model.S(:,rxnIdx);
119-
metsStoch(Hidx) = 0;
120-
model.S(Hidx,rxnIdx) = -sum(metsStoch.*model.metCharges,'omitnan');
121-
% Ion
122-
rxnIdx = strcmp(model.rxns,'r_4599');
123-
metsStoch = model.S(:,rxnIdx);
124-
metsStoch(Hidx) = 0;
125-
model.S(Hidx,rxnIdx) = -sum(metsStoch.*model.metCharges,'omitnan');
100+
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
101+
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
102+
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
103+
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
104+
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
126105

127106
% Special case for SLIME rxns
128107
model.metCharges(find(contains(model.metNames,'chain')+contains(model.metNames,'backbone'))) = 0;
@@ -165,6 +144,12 @@
165144
% This section focuses on individual reactions that have the wrong
166145
% reversibility/direction/cofactor or should be completley removed
167146

147+
% Make GCY1 irreversible. Has a positive DeltaGo' (+20.9) and is part of a
148+
% transhydrogenase cycle (NADH -> NADPH) at the cost of one ATP. High
149+
% cytosolic NADPH/NADP ratio makes it thermodynamically infeasible that it
150+
% runs in reverse direction.
151+
model = setParam(model,'ub','r_0487',0);
152+
168153
% The mitochondrial ATP synthase is able to run in reverse, which occurs
169154
% in anaerobic conditions
170155
model = setParam(model,'lb','r_0226',-1000);
@@ -197,6 +182,40 @@
197182
model = setParam(model,'lb',{'r_4723','r_4724','r_4725'},0);
198183
model = setParam(model,'lb','r_4460',0);
199184

185+
% While r_0013 was elementary balanced, it was not charged balanced. The
186+
% reaction mechanism was incorrect. Corrected to mimic a combination of
187+
% MetaCyc rxns: R83-RXN and R147-RXN; or KEGG rxns: R07364 and R07395.
188+
model = changeRxns(model,'r_0013','5-(methylsulfanyl)-2,3-dioxopentyl phosphate[c] + H2O[c] + oxygen[c] => 4-methylthio-2-oxobutanoate[c] + formate[c] + 2 H+[c] + phosphate[c]',3);
189+
190+
% Represent ACP with formula "RHS"
191+
model.metFormulas(getIndexes(model,'s_1845','mets')) = {'RHS'};
192+
193+
% [~,metFormulae] = computeMetFormulae(model,'metMwRange','s_0338','fillMets','none')
194+
% model.metFormulas(getIndexes(model,'s_0329','mets')) = {'C17H28NO17P'};
195+
% model.metFormulas(getIndexes(model,'s_0330','mets')) = {'C33H58NO18P'};
196+
% model.metFormulas(getIndexes(model,'s_0331','mets')) = {'C19H30NO18P'};
197+
% model.metFormulas(getIndexes(model,'s_0334','mets')) = {'C39H68NO23P'};
198+
% model.metFormulas(getIndexes(model,'s_0337','mets')) = {'C44H78N2O27P2'};
199+
% model.metFormulas(getIndexes(model,'s_0338','mets')) = {'C41H38N2O41P3'};
200+
% model.metFormulas(getIndexes(model,'s_0339','mets')) = {'C50H88N2O32P2'};
201+
202+
% Copy annotations between same metabolites in separate compartments
203+
model.metFormulas(getIndexes(model,'s_4211','mets')) = model.metFormulas(getIndexes(model,'s_2885','mets'));
204+
model.metCharges(getIndexes(model,'s_4211','mets')) = model.metCharges(getIndexes(model,'s_2885','mets'));
205+
model.metFormulas(getIndexes(model,'s_4209','mets')) = model.metFormulas(getIndexes(model,'s_3826','mets'));
206+
model.metCharges(getIndexes(model,'s_4209','mets')) = model.metCharges(getIndexes(model,'s_3826','mets'));
207+
model.metMiriams(getIndexes(model,'s_4209','mets')) = model.metMiriams(getIndexes(model,'s_3826','mets'));
208+
209+
% r_4323 is an less precise half-reaction of r_4324 and will be removed
210+
model = removeReactions(model,'r_4323',true,true,true);
211+
212+
% r_4325 represents scaffolding during [Fe-S]-cluster synthesis, not a
213+
% metabolic process, and will therefore be removed
214+
model = removeReactions(model,'r_4325',true,true,true);
215+
216+
% r_0229
217+
%model = changeRxns(model,'r_0229','dethiobiotin[c] + polysulphur[c] <=> biotin[c] + 2 H+[c]',2)
218+
%dethiobiotin[c] + hydrogen sulfide[c] + 2 S-adenosyl-L-methionine[c] + 2 H+[c] <=> biotin[c] + 2 L-methionine[c] + 2 5'-Deoxyadenosine
200219
%% ========================================================================
201220
% Condition-specific gene expression. These can be enabled with scripts
202221
% Glycine cleavage only active when glycine is used as nitrogen source
@@ -214,6 +233,23 @@
214233
model = setParam(model,'eq',{'r_0252'},0);
215234
model.rxnNotes(ismember(model.rxns,{'r_0252'})) = {'Only active if growth medium contains carnitine'};
216235

236+
%% Rescale protein fraction so that biomass sums up to 1 g/gDCW
237+
% Protein is the largest fraction, so increasing
238+
[X,P] = sumBioMass(model, false);
239+
fprintf('Current biomass adds up to %.4f g/g. Protein fraction is scaled from %.4f to %.4f g/g to reach 1 g/g total biomass.\n', X, P, (1-X)+P)
240+
model = scaleBioMass(model,'biomass',1,'protein');
241+
242+
%% Degree of reduction of biomass
243+
% To align the degree of reduction of S. cerevisiae biomass to the
244+
% published value of 4.2 /Cmol (Lange and Heijnen, 2001, 10.1002/bit.10054)
245+
246+
DR = 75; % 3mmol (g CDW)−1s
247+
metIdx = getIndexes(model,{'s_1212','s_1207','s_0794'},'mets'); % NADPH[c], NADP[c], H+[c]
248+
bioIdx = getIndexes(model,'r_4041','rxns');
249+
250+
currCoeff = full(model.S(metIdx,bioIdx)); % Gather the current coefficients
251+
model.S(metIdx,bioIdx) = currCoeff + [-DR; +DR; -DR];
252+
217253
%% ========================================================================
218254

219255
%% DO NOT CHANGE OR REMOVE THE CODE BELOW THIS LINE.

‎code/modelTests/anaerobiosis.m‎

Lines changed: 29 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -1,67 +1,53 @@
11
clear; close all
22
% Load the model and apply the corrections for *all* models
3-
model = readYAMLmodel('../../model/yeast-GEM.yml');
3+
%model = loadYeastModel;
44
cd('../modelCuration/')
55
v9_1_0;
66

7-
% res=essentialGenes(model);
8-
% res
9-
107
%% Run growth tests
118
R2 = growth(model);
129

1310
%% Convert to anaerobic
1411
cd('../otherChanges/')
15-
model = anaerobicModel(model);
12+
modelAn = anaerobicModel(model);
1613

1714
cd('../modelTests/');
1815
%% flux predictions
19-
R2=anaerobic_flux_predictions(model);
16+
R2=anaerobic_flux_predictions(modelAn);
2017

18+
%% Plot
19+
plotAnaerobic(modelAn)
20+
%% Set glucose uptake rate and solve pFBA
21+
modelAn = setParam(modelAn,'eq','r_1714',-23);
22+
res=solveLP(modelAn,1);
2123

24+
v_AStr = res.x(getIndexes(modelAn,'r_1115','rxns'))
25+
v_ATPase = res.x(getIndexes(modelAn,'r_0227','rxns'))
2226

2327

24-
%% Set glucose uptake rate and solve pFBA
25-
model = setParam(model,'eq','r_1714',-23);
26-
res=solveLP(model,1);
27-
FLUX=res.x;
28-
29-
30-
%% Retrieve data for the main products
31-
v_glc=FLUX(getIndexes(model,'r_1714','rxns'),:);
32-
v_eth=FLUX(getIndexes(model,'r_1761','rxns'),:);
33-
v_CO2=FLUX(getIndexes(model,'r_1672','rxns'),:);
34-
v_gly=FLUX(getIndexes(model,'r_1808','rxns'),:);
35-
v_growth=FLUX(getIndexes(model,'r_4041','rxns'),:);
36-
37-
38-
%% Show relative accuracy of main extracellular products
39-
figure;
40-
%glycerol ethanol Co2
41-
%4.5 ± 0.4 31 ± 2 38 ± 10
42-
data=[4.5 31 38 0.36];
43-
sim=[v_gly v_eth v_CO2 v_growth];
44-
errorVal=[0.4 2 10 0.02];
45-
b1=bar(data./data,'FaceAlpha',0.5);hold on;b2=bar(sim./data,'FaceAlpha',0.5);
46-
hold on
47-
er = errorbar([1 2 3 4],data./data,errorVal./data,errorVal./data);
48-
er.Color = [0 0 0];
49-
er.LineStyle = 'none';
50-
legend({'data','simulation'});
51-
ylabel('Relative value');
52-
xticklabels({'Glycerol','Ethanol','CO2','Biomass'})
53-
5428
%% Pack flux results into table
55-
temp_model=model;
56-
temp_model.metNames=strcat( model.metNames, repmat('[',length(model.mets),1),model.compNames(temp_model.metComps),repmat(']',length(model.mets),1) );
57-
rxns_reacs=printRxnFormula(temp_model,'rxnAbbrList',model.rxns,'metNameFlag',1,'printFlag',0);
58-
tab=table(model.rxns,model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,model.grRules);
59-
tab = sortrows(tab,"Var4","descend");
29+
temp_model=modelAn;
30+
temp_model.metNames=strcat( modelAn.metNames, repmat('[',length(modelAn.mets),1),modelAn.compNames(temp_model.metComps),repmat(']',length(modelAn.mets),1) );
31+
rxns_reacs=printRxnFormula(temp_model,'rxnAbbrList',modelAn.rxns,'metNameFlag',1,'printFlag',0);
32+
tab=table(modelAn.rxns,modelAn.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,modelAn.grRules);
33+
%tab = sortrows(tab,"Var4","descend");
34+
35+
[massImbalance, imBalancedMass, imBalancedCharge, imBalancedRxnBool, elements, missingFormulaeBool, balancedMetBool] = checkMassChargeBalance(modelAn);
36+
37+
tabImbalance = [tab,table(imBalancedRxnBool,imBalancedCharge,imBalancedMass)];
38+
tabImbalance(~imBalancedRxnBool,:) = [];
39+
tabImbalance(contains(tabImbalance.Var2,'exchange'),:) = [];
40+
tabImbalance(contains(tabImbalance.Var2,'SLIME'),:) = [];
41+
%tab=table(model.rxns,model.rxnNames,rxns_reacs,abs(FLUX./v_glc),FLUX,model.grRules);
42+
43+
44+
45+
6046

6147

6248
%% Calculate formula and degree of reduciton of biomass
63-
[mwRange,metFormulae,elements,metEle]=computeMetFormulae(model,'metMwRange','s_0450','fillMets','none','printLevel',0);
64-
Biomass_index = find(strcmp(model.metNames,'biomass'));
49+
[mwRange,metFormulae,elements,metEle]=computeMetFormulae(modelAn,'metMwRange','s_0450','fillMets','none','printLevel',0);
50+
Biomass_index = find(strcmp(modelAn.metNames,'biomass'));
6551
%Degree of reduction per element. Order of the elements 'C', 'H', 'O', 'N'
6652
DR_per_ele = [4, 1, -2, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
6753
DR_per_Cmol=sum(metEle(Biomass_index,:).*DR_per_ele)/metEle(Biomass_index,1)

‎code/modelTests/growth.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -86,7 +86,7 @@
8686
model_origin = scaleBioMass(model_origin,'protein',0.289,'',false);
8787
model_origin = scaleBioMass(model_origin,'lipid',0.048,'',false);
8888
model_origin = scaleBioMass(model_origin,'RNA',0.077,'carbohydrate',false);
89-
89+
model_origin = setParam(model_origin,'ub','r_0472',1000); %Glutamate synthase repressed in excess nitrogen
9090
end
9191

9292
if mode1 == 2

‎code/otherChanges/anaerobicModel.m‎

Lines changed: 2 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,6 @@
1616
%
1717
% Usage: model = anaerobicModel(model)
1818

19-
20-
2119
%% Set environmental conditions
2220
% Remove heme a from the cofactor pseudoreaction (part of biomass)
2321
hemeIdx = getIndexes(model,'s_3714','mets');
@@ -39,17 +37,6 @@
3937
model.lb(strcmp(model.rxns,'r_1967')) = -1000; %nicotinate
4038
model.lb(strcmp(model.rxns,'r_1548')) = -1000; %(R)-pantothenate
4139

42-
%% Degree of reduction of biomass
43-
% To align the degree of reduction of S. cerevisiae biomass to the
44-
% published value of 4.2 /Cmol (Lange and Heijnen, 2001, 10.1002/bit.10054)
45-
46-
DR = 3; % 3mmol (g CDW)−1s
47-
metIdx = getIndexes(model,{'s_1212','s_1207','s_0794'},'mets'); % NADPH[c], NADP[c], H+[c]
48-
bioIdx = getIndexes(model,'r_4041','rxns');
49-
50-
currCoeff = full(model.S(metIdx,bioIdx)); % Gather the current coefficients
51-
model.S(metIdx,bioIdx) = currCoeff + [-DR; +DR; -DR];
52-
5340
%% Curations that are required to reach correct metabolic phenotypes during
5441
% anaerobic batch growth on minimal glucose media
5542

@@ -59,13 +46,7 @@
5946
% et al (2005) 10.1074 /jbc.M410573200) and not detected in proteome
6047
% (Sjöberg et al (2023) 10.1016/j.ymben.2024.01.007).
6148
model = setParam(model,'eq','r_0714',0);
62-
63-
% Make GCY1 irreversible. Has a positive DeltaGo' (+20.9) and is part of a
64-
% transhydrogenase cycle (NADH -> NADPH) at the cost of one ATP. High
65-
% cytosolic NADPH/NADP ratio makes it thermodynamically infeasible that it
66-
% runs in reverse direction.
67-
model = setParam(model,'ub','r_0487',0);
68-
49+
6950
% Block IDP2. It is strongly repressed in transcriptome (Tai et al (2005)
7051
% 10.1074 /jbc.M410573200) and not detected in proteome (Sjöberg et al
7152
% (2023) 10.1016/j.ymben.2024.01.007).
@@ -77,6 +58,7 @@
7758

7859
FADH2_prod=0.08;
7960
metIdx = getIndexes(model,{'s_0689','s_0687','s_0794'},'mets'); % FADH2[c], FAD[c], H+[c]
61+
bioIdx = getIndexes(model,'r_4041','rxns');
8062

8163
currCoeff = full(model.S(metIdx,bioIdx)); % Gather the current coefficients
8264
model.S(metIdx,bioIdx) = currCoeff + [FADH2_prod; -FADH2_prod; -2*FADH2_prod];

0 commit comments

Comments
 (0)