A problem with the Fig.2b in the article of yeast8

[vivarock233-byte]

I just wanted to repeat the outcome of Fig2b. in the article of yeast8 https://doi.org/10.1038/s41467-019-11581-3

For this figure, I tried to simulate the FCC（ flux control coefficient ）using the model of ecYeast9. But I just couldn't to repeat the same outcome. And even the growth in the only carbon source of glucitol or ribose is impossible. Could you give me a code about how to simulate this outcome, thanks~

The model I used is yeast 9.0.2 constructed with GECKO 3.2.4 after stage 3 or step 52 in https://github.com/SysBioChalmers/GECKO/blob/main/tutorials/full_ecModel/protocol.m

[edkerk]

A few things to note:

• The figure in the paper is from yeast8, you try to use yeast9, so that the outcome is not the same might be anticipated.
• The figure in the paper is from an ec-model generated by GECKO 1, that would have different kcat values from the GECKO 3 version that you are attempting now.

So the same code that would work for GECKO1-yeast8 might not work for GECKO3-yeast9.

Can you please provide the code that you are trying to run, identify the problem with this, and mention what result is expected. Otherwise it is very hard to troubleshoot your question.

[vivarock233-byte]

The first part is to look out what nutrition it take

adapterLocation = fullfile(findGECKOroot,'ecModelGEM9','ecModelGEM9Adapter.m');
ModelAdapter = ModelAdapterManager.setDefault(adapterLocation);
ModelAdapter = ModelAdapterManager.getDefault();
params = ModelAdapter.getParameters();
model=loadEcModel('ecYeastGEM_stage3.yml');

sol = solveLP(model, 1);
exc_idx = find(contains(model.rxns, 'r_1714') | contains(model.rxnNames, 'exchange')); 
for i = 1:length(exc_idx)
    idx = exc_idx(i);
    flux = sol.x(idx);
    if model.lb(idx) < 0 && abs(flux) > 1e-6
        rxnName = model.rxnNames{idx};       
            fprintf('%s (Flux: %.4f)\n', rxnName, flux);
    end
end

and then get the output

ammonium exchange (Flux: -2.7580)
D-glucose exchange (Flux: -23.0195)
H+ exchange (Flux: -0.1226)
iron(2+) exchange (Flux: -0.0000)
oxygen exchange (Flux: -1.2004)
phosphate exchange (Flux: -0.1204)
potassium exchange (Flux: -0.0016)
sodium exchange (Flux: -0.0017)
sulphate exchange (Flux: -0.0371)
water exchange (Flux: 9.6263)
chloride exchange (Flux: -0.0006)
Cu2(+) exchange (Flux: -0.0003)
Mn(2+) exchange (Flux: -0.0012)
Zn(2+) exchange (Flux: -0.0003)
Mg(2+) exchange (Flux: -0.0005)
Ca(2+) exchange (Flux: -0.0001)
prot_pool_exchange (Flux: -125.0000)

then I tried to simulate the outcome using the code

carbon_sources = {
    'Glucose',  'r_1714';
    'Acetate',  'r_1634';
    'Ethanol',  'r_1761';
    'Glycerol', 'r_1808';
    'Glucitol', 'r_1713';
    'Galactose', 'r_1710';
    'Ribose',    'r_2020';
};

target_enzymes = {
    'Tdh1 (P00360)',  'P00360';
    'Pgi1 (P12709)',  'P12709'; 
    'Pfk1 (P16861)',  'P16861'; 
    'Fba1 (P14540)',  'P14540';
    'Cdc19 (P00549)', 'P00549'; % Pyruvate kinase
    % -ATP Synthase
    'Oli1 (P00853)',  'P00853'; % ATP synthase subunit 9 
    'Atp1 (P07251)',  'P07251'; % F1 alpha subunit 
    'Atp2 (P00830)',  'P00830'; % F1 beta subunit
    % Respiration
    'Cox1 (P00401)',  'P00401'; % Cytochrome c oxidase subunit 1
    'Cit1 (P00814)',  'P00814'; % Citrate synthase 
};

num_sources = size(carbon_sources, 1);
num_enzymes = size(target_enzymes, 1);
fcc_matrix = zeros(num_enzymes, num_sources);

for c = 1:num_sources
    source_name = carbon_sources{c, 1};
    source_rxn  = carbon_sources{c, 2};
    
    fprintf('\nsimulation of the carbon source : %s ...\n', source_name);
    temp_model = model; 
    temp_model = setParam(temp_model, 'lb', 'r_1714', 0); 
    
    % 2. open the current carbon source 
    temp_model = setParam(temp_model, 'lb', source_rxn, -1000);
    sol_ref = solveLP(temp_model, 1);

    if isempty(sol_ref.x) || sol_ref.f < 1e-6
        fprintf('  %f it cannot growth \n',sol_ref.x);
        continue;
    end
    
    mu_ref = full(sol_ref.f); 
    fprintf('  basic rate : %.4f h^-1\n', mu_ref);
    
    % calculate the FCC
    for e = 1:num_enzymes
        gene_name = target_enzymes{e, 1};
        prot_id   = target_enzymes{e, 2};
        met_name  = ['prot_' prot_id]; %
        met_idx = find(strcmp(temp_model.mets, met_name));
        
        if isempty(met_idx)
            met_idx = find(contains(temp_model.mets, prot_id));
        end
        
        if isempty(met_idx)
            continue; 
        end
        
        rxns_consuming = find(temp_model.S(met_idx, :) < 0);
        if isempty(rxns_consuming) || sum(abs(sol_ref.x(rxns_consuming))) < 1e-9
            continue; 
        end
        

        perturbation = 10; 
        
        model_pert = temp_model;
        model_pert.S(met_idx, rxns_consuming) = temp_model.S(met_idx, rxns_consuming) * (1/perturbation);
        
        sol_pert = solveLP(model_pert, 1);
        
        if ~isempty(sol_pert.x)
            mu_new = full(sol_pert.f);
            
            % FCC 公式: (% change in flux) / (% change in kcat)
            % change in kcat = 1% = 0.01
            fcc_val = ((mu_new - mu_ref) / mu_ref) / (perturbation - 1);
            
            % to avoid the noise 
            if fcc_val < 1e-6, fcc_val = 0; end
            fcc_matrix(e, c) = fcc_val;
            
            if fcc_val > 0.01
                fprintf('    %s (FCC = %.4f)\n', gene_name, fcc_val);
            end
        end
    end
end

and the outcome is below , all the FCC I calculated are all the 0.

simulation of the carbon source : Glucose ...
  basic rate : 0.4310 h^-1
simulation of the carbon source : Acetate ...
  basic rate : 0.1216 h^-1
simulation of the carbon source : Ethanol ...
  basic rate : 0.1736 h^-1
simulation of the carbon source : Glycerol ...
  basic rate : 0.2054 h^-1
simulation of the carbon source : Glucitol ...
   it cannot grow 
simulation of the carbon source : Galactose ...
  basic rate : 0.2388 h^-1
simulation of the carbon source : Ribose ...
   it cannot grow 

I think maybe I get something wrong in the culture medium or the method I used to calculate the FCC . Thank you for your help !

[edkerk]

Just a quick question, have you tested if the non-ec-model (yeast9) can grow on glucitol and ribose?

[vivarock233-byte]

Ok, I find out yeast9 cannot grow on glucitol and ribose.
But so why dose it can grow in yeast8 , no in yeast9. Can I pretend to think that the fig2.b is not correct after the curation ?
In my opinion, the fig2.b represent some gene can affect the growth to some extent . If this gene can function , no matter how the model change, it still can work. It's all my confusion. Thanks !

[edkerk]

Interesting insight, there has been curation between yeast8 vs. yeast9 that "broke" glucitol and ribose catabolism. There have been quite some curations between these two models, so the best would be to narrow down between which two releases this happened. So I suggested you systematically go through the yeast-GEM released versions between 8.0.0 and the latest 9.0.2. You might find this function useful to load previous model version.