createCommModel

[jfoldi81]

I ran this code:
createCommModel(comm,options)
comm is a cell array of models that I have run previously so I do not believe that this issue lies with them.
options is a cell array of the biomass reaction names for each model as strings - {["R_biomass345"]} {["R_biomass536"]} {["R_biomass525"]}. I only did spBm because I wanted to try the simplest version of the function at first without adding the optional options.

This was the error:
Error: File: getCobraComParams.m Line: 42 Column: 29
Using the dot operator to index into the output of function 'paramDefault' is not supported.

Error in createCommModel (line 95)
rxnField,metField,geneField,compField,sepUtEx,addExRxns,verbFlag] = getCobraComParams(...

[shjchan]

Hi @jfoldi81 ,

Sorry for the late response. options is a matlab structure (see the documentation for the available fields). For example, try the following code:

model = readCbModel('iJO1366.mat');
options=struct();
options.spBm={'BIOMASS_Ec_iJO1366_core_53p95M'; 'BIOMASS_Ec_iJO1366_core_53p95M'};
[modelCom,infoCom,indCom, msg] = createCommModel({model;model},options);

[jfoldi81]

No problem at all for the late reply @shjchan ! The beginning of the semester is always busy. I ran the code you provided and ran into the same error as before. Please let me know if there's any other info regarding my setup that may be helpful in resolving the error.

[shjchan]

@jfoldi81
Can you screen shot or copy all the error messages? 'cause I have it run successfully. And please also paste the output of running generateSystemConfigReport.

[shjchan]

@jfoldi81
I suspect that might be due to stricter code checking in later Matlab version? I am using R2017b. But anyway, those lines are not needed. If you just comment them out (lines 37 - 44), it should run. You can git pull to update the code if you are using git.

[jfoldi81]

@shjchan
I commented out those lines and it works now, thank you so much! My other question was regarding the 'media' being used in the function. I want to test the community on different media conditions (modeled as exchange reactions with flux values just like in the normal COBRA toolbox) but could not tell if this was something to add in the 'options' or if there was a different place for me to add the media being used.
For single model simulations, I was using .txt files that contained exchange reaction IDs and flux values and wanted to know how to insert these into the multi-species model. It seems it may be through either the addMetCom or addRxnCom functions but I could not tell from the documentation if that was the case or if it was elsewhere (neither of those two had flux values as an input which is why it was unclear to me if they would be used for media).

[shjchan]

You don't need to add reactions in general to change the media.
printUptakeBound will show the uptake bounds (which are negative numbers) for single organism models.
printUptakeBoundCom will show uptake bounds for community models.
printFluxVector is useful in showing the exchange fluxes in a flux distribution (e.g., simulated by FBA).

And you can change the lower bounds for exchange reactions to change the media.
See 'Example 1: Calculating growth rates' in this tutorial to see how to use the function changeRxnBounds to change the media condition. But each model might have a different naming convention for the exchange reactions.
If you are working with community models building using createCommModel here, the good thing is that in the infoCom and indCom, there is all the information needed to locate exchange reactions. Ask if you have questions about that. Might not be intuitive to understand all of them.

[shjchan]

For the test code I just asked you to run for the E. coli model, that won't actually generate a useful community model because some parameters were not set correctly. Try this instead:

model = readCbModel('iJO1366.mat');
options=struct();
options.spBm={'BIOMASS_Ec_iJO1366_core_53p95M'; 'BIOMASS_Ec_iJO1366_core_53p95M'};
% give an abbreviation to the members, which will be used in naming the model elements
options.spAbbr = {'Ec1';'Ec2'};
% single community exchange reactions for both uptake and export
options.sepUtEx = false;
% IMPORTANT: correct identifier to tell what are considered as extracellular metabolites from individual models
options.metExId = '_e'
[modelCom,infoCom,indCom, msg] = createCommModel({model;model},options);

You shell see something like the attached picture:

The diagram visualizes how different spaces are connected through these defined compartments and exchange reactions.

[jfoldi81]

@shjchan
I used the modified code you wrote and was able to get the same results as you posted. I had a couple of follow up questions.

First, I'm getting this error for one of my models that does not have an exchange reaction for NO3. Is there a way to bypass this error/change the metabolites which it tests? I haven't been using NO3 as a compound when running the model otherwise so I have not run into this issue before.

Compile metabolites and reactions...
Input model 2 org2:
	Exchange reaction added for no3[e], EX_no3(e): no3[e] <=>
Error using cell/ismember (line 34)
Input A of class cell and input B of class cell must be cell arrays of character vectors, unless one is a
character vector.

Error in findRxnIDs (line 18)
    [tmp,rxnID] = ismember(rxnList,model.rxns);

Error in createCommModel (line 280)
            spBm(j) = findRxnIDs(modelCell{j}, spBm0(j));

Error in steadyCOMpractice (line 23)
[modelCom,infoCom,indCom, msg] = createCommModel({sbmls(1);sbmls(2)},options);

Second, I tried running a community with the same model twice (different than the one with the NO3 error) in a similar vein to the double e. coli model, but am running into this error:

Compile metabolites and reactions...
Add community metabolites and exchange reactions...
Merge gene-reaction rules...
Array indices must be positive integers or logical values.

Error in infoCom2indCom (line 72)
        indCom.spBm = modelCom.rxns(infoCom.spBm);

Error in createCommModel (line 776)
infoCom = infoCom2indCom(modelCom,indCom,true,spAbbr,spName);

Error in steadyCOMpractice (line 23)
[modelCom,infoCom,indCom, msg] = createCommModel({sbmls(1);sbmls(1)},options);

Last, after using your code for the double e. coli model, I tried running SteadyCom by using the code from your tutorial and swapping EcCom for modelCom.

options = struct();
options.GRguess = 0.5; % initial guess for max. growth rate
options.GRtol = 1e-6; % tolerance for final growth rate
options.algorithm = 3; % use the default algorithm (simple guessing for bounds, followed by matlab fzero)
[sol, result] = SteadyCom(modelCom, options);

Running this immediately gave a Model infeasible at maintenance. Time elapsed: 1 / 1 sec message and I was not sure if I was missing an additional parameter to get SteadyCOM to work or if for this sample community this was the expected result since it is built on two identical models.

If you need any more information or files for these issues please let me know and I'm happy to provide any more information that would be helpful to you!

[shjchan]

@jfoldi81 ,

It will be good for you if you learn to debug the code. Usually what you can do is to add a break point in the matlab editor at the line where an error occurs. Then run the code again and Matlab will stop at that line and then you can look at the variables at that line to see why there is an error. See here.

And it is important to understand the error message. here in the first error:

Error using cell/ismember (line 34)
Input A of class cell and input B of class cell must be cell arrays of character vectors, unless one is a
character vector.

Error in findRxnIDs (line 18)
    [tmp,rxnID] = ismember(rxnList,model.rxns);

Error in createCommModel (line 280)
            spBm(j) = findRxnIDs(modelCell{j}, spBm0(j));

Input A of class cell and input B of class cell must be cell arrays of character vectors, unless one is a
character vector.
So this suggests the input variable spBm0 is problematic. Then you should trace back the code how spBm0 is defined. I guess probably it is an empty vector [] maybe because you didn't set that correctly in the input options.spBm. I am not sure exactly why this happens maybe a little more sanity check is needed to prevent this. It has nothing to do with the message about NO3.

For the second error, array indices must be positive integers or logical values. This suggests that again the index spBm in indCom is not correct. Make sure you input a correct biomass reaction for each model.

For the last thing, whenever you encounter

model infeasible

no matter you are using FBA, FVA or any functions, make sure you have set up the uptake media correctly (see the previous comment about using printUptakeBoundCom). The model returned by createCommModel has all the uptake shut down, therefore the infeasibility if you run steadyCom directly.

[shjchan]

@jfoldi81

I actually have a function for finding minimal media findMinimalMedia. Think it may be useful for you. I just added that to this repo.

[jfoldi81]

@shjchan
Thank you so much for the debugging tips! I'm still getting used to MATLAB and I was unaware of the breakpoint functionality which is unbelievably helpful. As a result, I figured out the first two issues while also getting a lot of MATLAB debugging practice which will serve me well in the future.
And thank you for adding findMinimalMedia! My current approach can get a bit time-intensive so I will be sure to give yours a try. I checked its documentation and that of solveCobraMILP to get a general idea of how it works but may have a few questions when I try implementing it. Thanks again!

[shjchan]

Something you need to be aware of: to solve MILP for larger models (say

100 reactions), practically you will need to install gurobi or cplex (both

free for students). Using the glpk that comes with the cobra toolbox installation almost won't do you any good.

…

_______________________________________________ Siu Hung Joshua Chan Assistant Professor Department of Chemical and Biological Engineering Colorado State University Phone: 970-491-6078 Office: Scott 254 Mailing address: 1370 Campus Delivery, Fort Collins, CO 80523 Schedule a meeting using this link <https://freebusy.io/shjchan@rams.colostate.edu>.

On Sat, Sep 5, 2020 at 12:23 AM jfoldi81 ***@***.***> wrote: @shjchan <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fshjchan&data=02%7C01%7Cjoshua.chan%40colostate.edu%7C70b55bd1cff7410971ea08d85164249d%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637348837849774627&sdata=Dqgc6zs5OVt6DjfyFz4MRt1BAhM6d671yMR0V2Hfv%2FE%3D&reserved=0> Thank you so much for the debugging tips! I'm still getting used to MATLAB and I was unaware of the breakpoint functionality which is unbelievably helpful. As a result, I figured out the first two issues while also getting a lot of MATLAB debugging practice which will serve me well in the future. And thank you for adding findMinimalMedia! My current approach can get a bit time-intensive so I will be sure to give yours a try. I checked its documentation and that of solveCobraMILP to get a general idea of how it works but may have a few questions when I try implementing it. Thanks again! — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fchan-csu%2FinhouseCobraFunctions%2Fissues%2F1%23issuecomment-687559671&data=02%7C01%7Cjoshua.chan%40colostate.edu%7C70b55bd1cff7410971ea08d85164249d%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637348837849774627&sdata=xCCVa%2B7EWevX3QQ%2BxbheMrYsYuFpwAWaItcjtauvW3E%3D&reserved=0>, or unsubscribe <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2Fnotifications%2Funsubscribe-auth%2FABQC5ZV6IGGOEIPDQIDXVJ3SEHKMLANCNFSM4QKC7FPA&data=02%7C01%7Cjoshua.chan%40colostate.edu%7C70b55bd1cff7410971ea08d85164249d%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637348837849784613&sdata=pna0mwzeQMghNAQzmrrOA05JXQytkU2P3ThaY7TvSbE%3D&reserved=0> .