Tutorial 4 — Fix an erroneous model¶

The power of GEMs comes from their size — but that size makes errors almost inevitable, whether you build a model yourself or use someone else's. This exercise is a systematic round of quality control on a deliberately broken version of the small yeast model, smallYeastBad.yml.

Model validation is iterative: some errors only surface once others are fixed. The RAVEN philosophy is that it is more valuable to try to make a model do something it should not (e.g. create mass from nothing) than to test only for what it should do.

tutorial4.m leaves several gaps for you to fill; the completed version is in tutorial4_solutions.m. Many of these checks are also wrapped by the gapReport function, but doing them step by step is far more instructive.

Note

Many of these changes are easier to do in the Excel sheet. They are done here in code just to avoid having several model files. It is assumed that you have already completed Tutorials 2–3 and are somewhat familiar with linear programming.

Step by step¶

1. Can the model make something from nothing?¶

Import the broken model, close all uptake, and maximise the sum of the producing exchange reactions. Any non-zero solution is a problem.

model = readYAMLmodel('smallYeastBad.yml');
model = setParam(model, 'eq', {'glcIN', 'o2IN'}, [0 0]);
model = setParam(model, 'obj', {'acOUT' 'biomassOUT' 'co2OUT' 'ethOUT' 'glyOUT'}, [1 1 1 1 1]);
sol = solveLP(model);
printFluxes(model, sol.x, true); %Nothing produced, good

2. Provoke hidden errors¶

An error may be masked because it costs energy or redox power. Add temporary "free" cofactor reactions (ATP ⇌ ADP + Pᵢ, and similar for NADH/NADPH) to provoke the model, then minimise the sum of fluxes for an interpretable solution:

rxns.rxns = {'FREE_ATP'; 'FREE_NADH'; 'FREE_NADPH'};
rxns.equations = {'ATP <=> ADP + phosphate'; 'NAD(+) <=> NADH'; 'NADP(+) <=> NADPH'};
model = addRxns(model, rxns, 2, 'c');
sol = solveLP(model, 1);
printFluxes(model, sol.x, false, [], [], '%rxnID (%rxnName):%flux\n\t%eqn\n');

Answer to Question 2 — production of ethanol from nothing

Lots of ethanol is produced. ADH1 should only produce one unit of ethanol. Change its equation with model = changeRxns(model, 'ADH1', 'acetaldehyde[c] + NADH[c] => ethanol[c] + NAD(+)[c]', 3);

3. Allow every metabolite to be excreted¶

Add a second column to model.b so RAVEN treats it as lower/upper bounds on the equality constraints — letting anything be excreted. This exposes errors that need a partner metabolite to be dumped.

model.b = [model.b inf(numel(model.b), 1)];
sol = solveLP(model, 1);
printFluxes(model, sol.x, false, 10^-5, [], '%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');

Answer to Question 3 — two unbalanced reactions

By looking at the reactions which were unbalanced and that were in the flux list, one can see that two reactions should each result in only one unit of product:

FBP should give one unit of F6P: change it to beta-D-fructofuranose 1,6-bisphosphate[c] => beta-D-fructofuranose 6-phosphate[c] + phosphate[c].
PFK should give one unit of F16P: change it to ATP[c] + beta-D-fructofuranose 6-phosphate[c] => ADP[c] + beta-D-fructofuranose 1,6-bisphosphate[c].

Apply each with changeRxns(model, …, 3).

4. `canConsume`¶

Set all uptakes and production to 0, drop the excretion column again, and check which metabolites can be consumed without any production:

model = setParam(model, 'eq', getExchangeRxns(model), 0);
model.b = model.b(:, 1);
I = canConsume(model);
disp(model.mets(I)); %These 12 metabolites can be consumed without any production

canConsume reports the metabolites the model can consume even when no production is allowed. Allow all uptake again, then force the uptake of one of them — here CO2 — and study the fluxes (a negative output means input):

model.b = [ones(numel(model.b), 1) * -1000 model.b];
model = setParam(model, 'eq', {'co2OUT'}, -1);
sol = solveLP(model);
printFluxes(model, sol.x, false, 10^-5, [], '%rxnID (%rxnName):\n\t%eqn\n\t%flux\n'); %Now it works

Answer to Question 4 — PDC is missing a product

PDC converts pyruvate (3 carbons) to acetaldehyde (2 carbons) without any other products — CO2 is missing. This would be simpler to change in the Excel file (or using changeRxns), but as an exercise one can add the cytosolic CO2 coefficient directly to the stoichiometric matrix:

Irxn = ismember(model.rxns, 'PDC');
Imet = ismember(model.mets, 'CO2_c');
model.S(Imet, Irxn) = 1; %The coefficient is 1.0
constructEquations(model, Irxn)

The solution is now infeasible, meaning it is no longer possible to force uptake of CO2 without any output.

5. Gap-filling against a reference model¶

tutorial4.m also imports a reference model to draw reactions from, allows all excretion, and reports gaps before filling them:

model = readYAMLmodel('smallYeastBad.yml');
refModel = readYAMLmodel('smallYeast.yml');
model.b = [model.b inf(numel(model.b), 1)];
sol = solveLP(model, 1);
printFluxes(model, sol.x, false, 10^-5, [], '%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');
I = canConsume(model);
disp(model.mets(I));
gapReport(model, {refModel});

[newConnected, cannotConnect, addedRxns, newModel] = fillGaps(model, {refModel}, false);
disp(addedRxns);
disp(newConnected);
disp(cannotConnect);

gapReport automates this whole quality-control workflow, while fillGaps adds reactions from the reference model (here the Tutorial 3 small yeast model) to connect the remaining gaps.

Full script¶

% tutorial4
%   This script contains the list of functions necessary for running
%   Tutorial 4, see Tutorial 4 in "RAVEN tutorials.docx" for more details.
%   Several key stages may be missing, try to fill these gaps before
%   checking the solutions in tutorial4_solutions. It is assumed that the
%   user is somewhat familiar with linear programming.
%
%   It is assumed that the user has already completed Tutorials 2-3

%Import the model to gap-fill and the reference model to draw reactions from
model=readYAMLmodel('smallYeastBad.yml');
refModel=readYAMLmodel('smallYeast.yml');
model.b=[model.b inf(numel(model.b),1)];
sol=solveLP(model,1);
printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');
I=canConsume(model);
disp(model.mets(I));
gapReport(model,{refModel});

[newConnected, cannotConnect, addedRxns, newModel]=fillGaps(model,{refModel},false);
disp(addedRxns);
disp(newConnected);
disp(cannotConnect);