cobratoolbox

cobratoolbox copied to clipboard

createMultipleSpeciesModel should, by default, if given only a single model for input, output a model that performs equivalently under FBA as the original model. For instance, given the model iAF1260, if I run optimizeCbModel on it, I get an objective value of 0.7367, though if I run createMultipleSpeciesModel first I get an objective value of 395.7353 - this doesn't seem ideal, though if there is a reason for this that we want to keep, it might be good to document the reasons for the difference, and possible how to recover the original value using createMultipleSpeciesModel. Note that currently this relates to #1227, though only tangentially I believe.

>> optimizeCbModel(iAF1260)

ans = 

  struct with fields:

         full: [2382×1 double]
          obj: 0.7367
        rcost: [2382×1 double]
         dual: [1668×1 double]
        slack: [1668×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0700
        basis: [1×1 struct]
            x: [2382×1 double]
            f: 0.7367
            y: [1668×1 double]
            w: [2382×1 double]
            v: [2382×1 double]

>> msModel =  createMultipleSpeciesModel({iAF1260}, {'modelIn'});
>> msModel.csense = repmat(['E'], 1, length(msModel.mets));
>> optimizeCbModel(msModel)

ans = 

  struct with fields:

         full: [2681×1 double]
          obj: 395.7353
        rcost: [2681×1 double]
         dual: [1967×1 double]
        slack: [1967×1 double]
       solver: 'gurobi'
    algorithm: 'default'
         stat: 1
     origStat: 'OPTIMAL'
         time: 0.0630
        basis: [1×1 struct]
            x: [2681×1 double]
            f: 395.7353
            y: [1967×1 double]
            w: [2681×1 double]
            v: [2681×1 double]

I hereby confirm that I have:

• [X] Tried to solve the issue on my own

Caveat: While we have looked at this a bit, I'm now at the stage of needing to look into it in more detail, as nothing obvious seemed to show up when last I checked.

• [X] Retried to run my code with the latest version of The COBRA Toolbox

Using the latest from master

• [X] Checked that a similar issue has not already been opened

This is due to createMultiSpeciesModel building a new extracellular environment, replacing the existing one. As such, there are no longer any uptake constraints, as it is in any instance likely, that the total uptakes and individual uptake constraints will have to be adjusted to the actual community.

@tpfau Thanks, I see what you mean - in fact each LB is lower than what I remember as the default: -999999 compared to -1000 (I'm quite rusty in general - I haven't been active with this software for nearly 4 years).

What you say makes sense, and if there is any more than one model, it becomes much more complicated.

Still, I think it might be useful for testing and validation purposes to have a way to do this. I can try to add this into the toolbox with some basic FBA/FVA tests if you think it would be worthwhile (we are thinking of doing it anyway for our own research).

You are welcome to add an option to the CreateMultiSpeciesModel functionality to allow existing flux bounds to be carried over between the model.

@tpfau @bbarker I think that is a good idea to allow existing bounds on exchange reactions to be carried over. It would need to be specified which of the joined models (if there is more than one) the exchange bounds should be copied from. Another important consideration: the exchange reactions in the joined model are created and named based on the extracellular metabolites in each original model. However, in some cases, the name of the exchange reaction in the original model is not "EX_metabolite ID(e)", but something different.

@almut-heinken
Personally, I think that any bounds should be only carried to the model they originate from and not to the overall model.
i.e. if model MA has an uptake max of 0.4 for compound X and model MB has one of 0.7, than the respective bounds of the model <-> extracellular should be adapted (i.e. the MAIEX_compound[u]tr reactions).

@tpfau That would be possible. However, we currently assume that the paired species can freely decide how to distribute the available nutrient resources and exchange metabolites with each other. That would no longer be true then.

@almut-heinken From my understanding, the exchange fluxes that are set in models are normally measured uptakes, or rather medium depletion measured when an individual microbe is cultured. As such this (at least to me) kind of indicates that this is the "maximal" uptake of the compound for the specific organism. If it would be total available nutrient, then of course you are right, but I haven't seen that type of constraint regularly, since most studies give data in mmol/gDW/h, which does not make any sense if its a absolute amount. And if its a bound on the actual uptake reaction (i.e. M1_glc[e] <-> M1_glc[c]), then it would anyways already be carried over. For maximal total available medium amounts I think that whoever uses the model will have to ensure that those are correct manually, since they are unlikely to be derivable from an individual model.