steadycom-based community metabolic network modeling-detailed code
0 Preface
Description of the content of this series: This series is an analysis and introduction to the types of models in biological research. Microbial Communities Modeling Each Other Series part of the . Introduce the metabolic network model in stages, from the construction of the entry-level environment to the sharing of the whole process of guiding the wet experiment process. Explore the interaction methods in the community based on biological models, and explore the evolutionary relationship between different organisms.
The directory article address of the series: https://blog.csdn.net/qq_33980829/article/details/118958463
1 Introduction to community metabolic network model
- Synthetic communities are usually assembled by microorganisms with multiple metabolic types, and there are multiple types of interaction patterns between microorganisms. The most typical ones include metabolic cross-mutual nutrition, substrate competition, and metabolic inhibition (antibiotics, antibacterial peptides, etc.). (To be added)
- Metabolic networks of microbial genome-based communities to study optimal growth rates under community co-cultivation conditions. (Genome sources can include: metagenomic prediction, NCBI download, homology modeling).
- Reference: Predicted Metabolic Function of the Gut Microbiota of Drosophila melanogaster10.1128/mSystems.01369-20
2. Modeling process
2.1 Obtaining the genome/metabolic network model of the target strain in the community
- When the target strain does not have a completed metabolic network (including the relationship between biomass and metabolites), the automatic modeling website is used for modeling.
- Predicting the linear relationship between biological growth rate and metabolite accumulation based on metabolic experiments (determining metabolites at a specific growth rate)
- Save metabolic network in standard sbml3 format or xls, mat
2.2 Read the metabolic network into the software and modify it manually
- Using cobar to read metabolic network in matlab
- Modify the metabolic label of the biomass in the read network to bio as the prefix e.g. biomass
- Check the metabolite names in the exocrine community in the metabolic network data, and modify it according to the standardized format
- Check the integrity of each strain's metabolism
2.3 Constructing the metabolic network into a community metabolic network to predict the metabolic rate of different strains in the community
- Establishment of community metabolism model using cobar tool
- Constraints for designing communities Maximum metabolic growth rate, basal rate
- Using steadycom to analyze the maximum growth rate under the metabolic balance of the community,
2.4 Analysis of the changes of metabolites in the community under different metabolic rates
- Analyze the amount of change in species at different growth rates
- The effect of cross-cotrophic strength in model communities on community stability
2.5 Add fitness function to constrain the community
3. Correction of official examples
3.1 Description
- The identification name of the biomass in the officially given case is not compatible with the actual matching process (the code identifies bio and it is Ec_bio in the improved metabolic model)
- There are incompatible parts of version changes, and some code needs to be commented
3.2 Modify the commented code
% function [sol, res] = ecoli_basicTest
iAF1260 = load('iAF1260.mat')
ecoli=iAF1260.model;
%model_1 = readCbModel('ecoli_1.2_cobra.xml')
% convert the compartment format from e.g., '_c' to '[c]'
ecoli.mets = regexprep(ecoli.mets, '_([^_]+)$', '\[$1\]');
% make all empty cells in cell arrays to be empty string
fieldToBeCellStr = {'metFormulas'; 'genes'; 'grRules'; 'metNames'; 'rxnNames'; 'subSystems'};
for j = 1:numel(fieldToBeCellStr)
if isfield(ecoli, fieldToBeCellStr{j})
ecoli.(fieldToBeCellStr{j})(cellfun(@isempty, ecoli.(fieldToBeCellStr{j}))) = {''};
end
end
ecoli = addReaction(ecoli,{'METt3pp',''},'met__L[c] + h[c] => met__L[p] + h[p]');
argH = {'ARGSL'}; % essential for arginine biosynthesis
lysA = {'DAPDC'}; % essential for lysine biosynthesis
metA = {'HSST'}; % essential for methionine biosynthesis
ilvE = {'PPNDH'}; % essential for phenylalanine biosynthesis
argO = {'ARGt3pp'}; % Evidence for an arginine exporter encoded by yggA (argO) that is regulated by the LysR-type transcriptional regulator ArgP in Escherichia coli.
lysO = {'LYSt3pp'}; % Distinct paths for basic amino acid export in Escherichia coli: YbjE (LysO) mediates export of L-lysine
yjeH = {'METt3pp'}; % YjeH is a novel L-methionine and branched chain amino acids exporter in Escherichia coli
yddG = {'PHEt2rpp'}; % YddG from Escherichia coli promotes export of aromatic amino acids.
% % auxotrophic for Lys and Met, not exporting Phe
Ec1 = ecoli;
Ec1 = changeRxnBounds(Ec1, [lysA; metA; yddG], 0, 'b');
% % auxotrophic for Arg and Phe, not exporting Met
Ec2 = ecoli;
Ec2 = changeRxnBounds(Ec2, [argH; yjeH; ilvE], 0, 'b');
% % Auxotrophic for Arg and Phe, not exporting Lys
Ec3 = ecoli;
Ec3 = changeRxnBounds(Ec3, [argH; lysO; ilvE], 0, 'b');
% % Auxotrophic for Lys and Met, not exporting Arg
Ec4 = ecoli;
Ec4 = changeRxnBounds(Ec4, [argO; lysA; metA], 0, 'b');
% extracellular metabolites (met[e])
metEx = strcmp(getCompartment(ecoli.mets),'e');
% the corresponding exchange reactions
rxnExAll = find(sum(ecoli.S ~= 0, 1) == 1);
[rxnEx, ~] = find(ecoli.S(metEx, rxnExAll)'); % need to be in the same order as metEx
rxnEx = rxnExAll(rxnEx);
% exchange rate
lbEx = ecoli.lb(rxnEx);
nameTagsModel = {'Ec1'; 'Ec2'; 'Ec3'; 'Ec4'};
EcCom = createMultipleSpeciesModel({Ec1; Ec2; Ec3; Ec4}, nameTagsModel);
%EcCom.csense = char('E' * ones(1,numel(EcCom.mets))); % correct the csense
clear Ec1 Ec2 Ec3 Ec4
[EcCom.infoCom, EcCom.indCom] = getMultiSpeciesModelId(EcCom, nameTagsModel);
disp(EcCom.infoCom);
rxnBiomass = strcat(nameTagsModel, 'Ec_biomass_iAF1260_core_59p81M'); % biomass reaction names
rxnBiomassId = findRxnIDs(EcCom, rxnBiomass); % ids
EcCom.infoCom.spBm = rxnBiomass; % .spBm for organism biomass reactions
EcCom.indCom.spBm = rxnBiomassId;
[yn, id] = ismember(strrep(ecoli.mets(metEx), '[e]', '[u]'), EcCom.infoCom.Mcom); % map the metabolite name
assert(all(yn)); % must be a 1-to-1 mapping
EcCom.lb(EcCom.indCom.EXcom(:,1)) = lbEx(id); % assign community uptake bounds
EcCom.ub(EcCom.indCom.EXcom(:,1)) = 1e5;
EcCom.lb(EcCom.indCom.EXsp) = repmat(lbEx(id), 1, 4); % assign organism-specific uptake bounds
%
%
% only allow to take up the amino acids that one is auxotrophic for
exRate = 1; % maximum uptake rate for cross feeding AAs
% Ec1
EcCom = changeRxnBounds(EcCom, {'Ec1IEX_arg__L[u]tr'; 'Ec1IEX_phe__L[u]tr'}, 0, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec1IEX_met__L[u]tr'; 'Ec1IEX_lys__L[u]tr'}, -exRate, 'l');
% Ec2
EcCom = changeRxnBounds(EcCom, {'Ec2IEX_arg__L[u]tr'; 'Ec2IEX_phe__L[u]tr'}, -exRate, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec2IEX_met__L[u]tr'; 'Ec2IEX_lys__L[u]tr'}, 0, 'l');
% Ec3
EcCom = changeRxnBounds(EcCom, {'Ec3IEX_arg__L[u]tr'; 'Ec3IEX_phe__L[u]tr'}, -exRate, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec3IEX_met__L[u]tr'; 'Ec3IEX_lys__L[u]tr'}, 0, 'l');
% Ec4
EcCom = changeRxnBounds(EcCom, {'Ec4IEX_arg__L[u]tr'; 'Ec4IEX_phe__L[u]tr'}, 0, 'l');
EcCom = changeRxnBounds(EcCom, {'Ec4IEX_met__L[u]tr'; 'Ec4IEX_lys__L[u]tr'}, -exRate, 'l');
% allow production of anything for each member
EcCom.ub(EcCom.indCom.EXsp(:)) = 1;
options = struct();
options.GRguess = 0.01; % initial guess for max. growth rate
options.GRtol = 1e-6; % tolerance for final growth rate
options.algorithm = 1; % use the default algorithm (simple guessing for bounds, followed by matlab fzero)
[sol, res] = SteadyCom(EcCom, options);
% FIXME, Note from Brandon Barker: the below code is not working and
% almost seems to be tacked on; haven't had a chance
% to ascertain its provenance as yet.
% %%
% growth = "growth";
% if model_1.rxns{find(model_1.c)} == ""
% model_1.rxns{find(model_1.c)} = 'growth';
% else
% growth = model_1.rxns(find(model_1.c));
% end
% model_2 = model_1;
% model_3 = model_1;
% % to avoid errors associated with PubChemID column in the createMultipleSpeciesModel
% %model_2.metPubChemID = num2cell(model_2.metPubChemID)
% %model_3.metPubChemID = num2cell(model_3.metPubChemID)
% multi_model = createMultipleSpeciesModel({model_1; model_3}, {'Org1'; 'Org2'});
% [multi_model.infoCom, multi_model.indCom] = getMultiSpeciesModelId(multi_model, {'Org1'; 'Org2'});
% multi_model.csense = char('E' * ones(1,numel(multi_model.mets))); % correct the csense
% disp( multi_model.infoCom)
% disp( multi_model.indCom)
% rxnBiomass = strcat({'Org1'; 'Org2'}, growth); % biomass reaction names
% rxnBiomassId = findRxnIDs(multi_model, rxnBiomass); % ids
% multi_model.infoCom.spBm = rxnBiomass; % .spBm for organism biomass reactions
% multi_model.indCom.spBm = rxnBiomassId;
options.optBMpercent = 100;
n = size(EcCom.S, 2); % number of reactions in the model
% options.rxnNameList is the list of reactions subject to FVA. Can be reaction names or indices.
% Use n + j for the biomass variable of the j-th organism. Alternatively, use {'X_j'}
% for biomass variable of the j-th organism or {'X_Ec1'} for Ec1 (the abbreviation in EcCom.infoCom.spAbbr)
options.rxnNameList = {'X_Ec1'; 'X_Ec2'; 'X_Ec3'; 'X_Ec4'};
options.optGRpercent = [89:0.2:99, 99.1:0.1:100]; % perform FVA at various percentages of the maximum growth rate, 89, 89.1, 89.2, ..., 100
[fvaComMin,fvaComMax] = SteadyComFVA(EcCom, options);
%%
% Similar to the output by |fluxVariability|, |fvaComMin| contains the minimum
% fluxes corresponding to the reactions in |options.rxnNameList|. |fvaComMax|
% contains the maximum fluxes. options.rxnNameList can be supplied as a (#rxns
% + #organism)-by-K matrix to analyze the variability of the K linear combinations
% of flux/biomass variables in the columns of the matrix. See |help SteadyComFVA|
% for more details.
%
% We would also like to compare the results against the direct use of FBA
% and FVA by calling |optimizeCbModel| and |fluxVariability|:
optGRpercentFBA = [89:2:99 99.1:0.1:100]; % less dense interval to save time because the results are always the same for < 99%
nGr = numel(optGRpercentFBA);
[fvaFBAMin, fvaFBAMax] = deal(zeros(numel(options.rxnNameList), nGr));
% change the objective function to the sum of all biomass reactions
EcCom.c(:) = 0;
EcCom.c(EcCom.indCom.spBm) = 1;
EcCom.csense = char('E' * ones(1, numel(EcCom.mets)));
s = optimizeCbModel(EcCom); % run FBA
grFBA = s.f;
for jGr = 1:nGr
fprintf('Growth rate %.4f :\n', grFBA * optGRpercentFBA(jGr)/100);
[fvaFBAMin(:, jGr), fvaFBAMax(:, jGr)] = fluxVariability(EcCom, optGRpercentFBA(jGr), 'max', EcCom.infoCom.spBm, 2);
end
%%
% Plot the results to visualize the difference (see also Figure 2 in ref.
% [1]):
grComV = result.GRmax * options.optGRpercent / 100; % vector of growth rates tested
lgLabel = {'{\itEc1 }';'{\itEc2 }';'{\itEc3 }';'{\itEc4 }'};
col = [235 135 255; 0 235 0; 255 0 0; 95 135 255 ]/255; % color
f = figure;
% SteadyCom
subplot(2, 1, 1);
hold on
x = [grComV(:); flipud(grComV(:))];
for j = 1:4
y = [fvaComMin(j, :), fliplr(fvaComMax(j, :))];
p(j, 1) = plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2);
p(j, 1).Color = col(j, :);
end
tl(1) = title('\underline{SteadyCom}', 'Interpreter', 'latex');
tl(1).Position = [0.7 1.01 0];
ax(1) = gca;
ax(1).XTick = 0.66:0.02:0.74;
ax(1).YTick = 0:0.2:1;
xlim([0.66 0.74])
ylim([0 1])
lg = legend(lgLabel);
lg.Box = 'off';
yl(1) = ylabel('Relative abundance');
xl(1) = xlabel('Community growth rate (h^{-1})');
% FBA
grFBAV = grFBA * optGRpercentFBA / 100;
x = [grFBAV(:); flipud(grFBAV(:))];
subplot(2, 1, 2);
hold on
% plot j=1:2 only because 3:4 overlap with 1:2
for j = 1:2
y = [fvaFBAMin(j, :), fliplr(fvaFBAMax(j, :))] ./ x';
% it is possible some values > 1 because the total biomass produced is
% only bounded below when calling fluxVariability. Would be strictly
% equal to 1 if sum(biomass) = optGRpercentFBA(jGr) * grFBA is constrained. Treat them as 1.
y(y>1) = 1;
p(j, 2)= plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2);
p(j, 2).Color = col(j, :);
end
tl(2) = title('\underline{Joint FBA}', 'Interpreter', 'latex');
tl(2).Position = [0.55 1.01 0];
ax(2) = gca;
ax(2).XTick = 0.52:0.02:0.58;
ax(2).YTick = 0:0.2:1;
xlim([0.52 0.58])
ylim([0 1])
xl(2) = xlabel('Community growth rate (h^{-1})');
yl(2) = ylabel('Relative abundance');
ax(1).Position = [0.1 0.6 0.5 0.32];
ax(2).Position = [0.1 0.1 0.5 0.32];
lg.Position = [0.65 0.65 0.1 0.27];
%%
% The direct use of FVA compared to FVA under the SteadyCom framework gives
% very little information on the organism's abundance. The ranges for almost all
% growth rates span from 0 to 1. In contrast, |SteadyComFVA| returns results with
% the expected co-existence of all four mutants. When the growth rates get closer
% to the maximum, the ranges shrink to unique values.
%
%
%% Analyze Pairwise Relationship Using SteadyComPOA
% Now we would like to see at a given growth rate, how the abundance of an organism
% influences the abundance of another organism. We check this by iteratively fixing
% the abundance of an organism at a level (independent variable) and optimizing
% for the maximum and minimum allowable abundance of another organism (dependent
% variable). This is what |SteadyComPOA| does.
%
% Set up the option structure and call |SteadyComPOA|. |Nstep| is an important
% parameter to designate how many intermediate steps are used or which values
% between the min and max values of the independent variable are used for optimizing
% the dependent variable. |savePOA| options must be supplied with a non-empty
% string or a default name will be used for saving the POA results. By default,
% the function analyzes all possible pairs in |options.rxnNameList|. To analyze
% only particular pairs, use |options.pairList|. See |help SteadyComPOA |for more
% details.
%%
options.savePOA = ['POA' filesep 'EcCom']; % directory and fila name for saving POA results
options.optGRpercent = [99 90 70 50]; % analyze at these percentages of max. growth rate
% Nstep is the number of intermediate steps that the independent variable will take different values
% or directly the vector of values, e.g. Nsetp = [0, 0.5, 1] implies fixing the independent variable at the minimum,
% 50% from the min to the max and the maximum value respectively to find the attainable range of the dependent variable.
% Here use small step sizes when getting close to either ends of the flux range
a = 0.001*(1000.^((0:14)/14));
options.Nstep = sort([a (1-a)]);
[POAtable, fluxRange, Stat, GRvector] = SteadyComPOA(EcCom, options);
%%
% POAtable is a _n_-by-_n_ cell if there are _n_ targets in |options.rxnNameList|.
% |POAtable{i, i}| is a _Nstep_-by-1-by-_Ngr_ matrix where _Nstep _is the number
% of intermediate steps detemined by |options.Nstep| and _Ngr _is the number of
% growth rates analyzed. |POAtable{i, i}(:, :, k)| is the values at which the
% _i_-th target is fixed for the community growing at the growth rate |GRvector(k)|.
% POAtable{i, j} is a _Nstep_-by-2-by-_Ngr_ matrix where |POAtable{i, j}(:, 1,
% k)| and |POAtable{i, j}(:, 2, k)| are respectively the min. and max. values
% of the _j_-th target when fixing the _i_-th target at the corresponding values
% in |POAtable{i, i}(:, :, k)|. |fluxRange |contains the min. and max. values
% for each target (found by calling |SteadyComFVA|). |Stat |is a _n_-by-_n-by-Ngr_
% structure array, each containing two fields: |*.cor|, the correlatiion coefficient
% between the max/min values of the dependent variable and the independent variable,
% and |*.r2|, the R-squred of linear regression. They are also outputed in the
% command window during computation. All the computed results are also saved in
% the folder 'POA' starting with the name 'EcCom', followed by 'GRxxxx' denoting
% the growth rate at which the analysis is performed.
%
% Plot the results (see also Figure 3 in ref. [1]):
nSp = 4;
spLab = {'{\it Ec1 }';'{\it Ec2 }';'{\it Ec3 }';'{\it Ec4 }'};
mark = {'A', 'B', 'D', 'C', 'E', 'F'};
nPlot = 0;
for j = 1:nSp
for k = 1:nSp
if k > j
nPlot = nPlot + 1;
ax(j, k) = subplot(nSp-1, nSp-1, (k - 2) * (nSp - 1) + j);
hold on
for p = 1:size(POAtable{1, 1}, 3)
x = [POAtable{j, j}(:, :, p);POAtable{j, j}(end:-1:1, :, p);...
POAtable{j, j}(1, 1, p)];
y = [POAtable{j, k}(:, 1, p);POAtable{j, k}(end:-1:1, 2, p);...
POAtable{j, k}(1, 1, p)];
plot(x(~isnan(y)), y(~isnan(y)), 'LineWidth', 2)
end
xlim([0.001 1])
ylim([0.001 1])
ax(j, k).XScale = 'log';
ax(j, k).YScale = 'log';
ax(j, k).XTick = [0.001 0.01 0.1 1];
ax(j, k).YTick = [0.001 0.01 0.1 1];
ax(j, k).YAxis.MinorTickValues=[];
ax(j, k).XAxis.MinorTickValues=[];
ax(j, k).TickLength = [0.03 0.01];
xlabel(spLab{j});
ylabel(spLab{k});
tx(j, k) = text(10^(-5), 10^(0.1), mark{nPlot}, 'FontSize', 12, 'FontWeight', 'bold');
end
end
end
lg = legend(strcat(strtrim(cellstr(num2str(options.optGRpercent(:)))), '%'));
lg.Position = [0.7246 0.6380 0.1700 0.2015];
lg.Box='off';
subplot(3, 3, 3, 'visible', 'off');
t = text(0.2, 0.8, {'% maximum';'growth rate'});
for j = 1:nSp
for k = 1:nSp
if k>j
ax(j, k).Position = [0.15 + (j - 1) * 0.3, 0.8 - (k - 2) * 0.3, 0.16, 0.17];
ax(j, k).Color = 'none';
end
end
end
%%
% There are two patterns observed. The two pairs showing negative correlations,
% namely Ec1 vs Ec4 (panel D) and Ec2 vs Ec3 (panel C) are indeed competing for
% the same amino acids with each other (Ec1 and Ec4 competing for Lys and Met;
% Ec2 and Ec4 competing for Arg and Phe). Each of the other pairs showing positive
% correlations are indeed the cross feeding pairs, e.g., Ec1 and Ec2 (panel A)
% cross feeding on Arg and Lys. See ref. [1] for more detailed discussion.
%
%
%% Parallelization and Timing
% SteadyCom in general can be finished within 20 iterations, i.e. solving 20
% LPs (usually faster if using Matlab |fzero|) for an accuracy of 1e-6 for the
% maximum community growth rate. The actual computation time depends on the size
% of the community metabolic network. The current |EcCom| model has 6971 metabolites
% and 9831 reactions. It took 18 seconds for a MacBook Pro with 2.5 GHz Intel
% Core i5, 4 GB memory running Matlab R2016b and Cplex 12.7.1.
%
% Since the FVA and POA analysis can be time-consuming for large models with
% a large number of reactions to be analyzed, SteadyComFVA and SteadyComPOA support
% parrallelization using the Matlab Distributed Computing Toolbox (|parfor| for
% SteadyComFVA and |spmd| for SteadyComPOA).
%
% Test SteadyComFVA with 2 threads:
%%
options.rxnNameList = EcCom.rxns(1:100); % test FVA for the first 50 reactions
options.optGRpercent = 99;
options.algorithm = 1;
options.threads = 1; % test single-thread computation first
options.verbFlag = 0; % no verbose output
tic;
[minF1, maxF1] = SteadyComFVA(EcCom, options);
t1 = toc;
if isempty(gcp('nocreate'))
parpool(2); % start a parallel pool
end
options.threads = 2; % two threads (0 to use all available workers)
tic;
[minF2, maxF2] = SteadyComFVA(EcCom, options); % test single-thread computation first
t2 = toc;
fprintf('Maximum difference between the two solutions: %.4e\n', max(max(abs(minF1 - minF2)), max(abs(maxF1 - maxF2))));
fprintf('\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n', t1, t2);
%%
% If there are many reactions to be analyzed, use |options.saveFVA| to give
% a relative path for saving the intermediate results. Even though the computation
% is interrupted, by calling |SteadyComFVA| with the same |options.saveFVA|, the
% program will detect previously saved results and continued from there.
%
% Test SteadyComPOA with 2 threads:
options.rxnNameList = EcCom.rxns(find(abs(result.flux) > 1e-2, 6));
options.savePOA = 'POA/EcComParallel'; % save with a new name
options.verbFlag = 3;
options.threads = 2;
options.Nstep = 5; % use a smaller number of steps for test
tic;
[POAtable1, fluxRange1] = SteadyComPOA(EcCom, options);
t3 = toc;
%%
% The parallelization code uses |spmd| and will redistribute jobs once any
% of the workers has finished to maximize the computational efficiency.
options.savePOA = 'POA/EcComSingeThread';
options.threads = 1;
tic;
[POAtable2, fluxRange2] = SteadyComPOA(EcCom, options);
t4 = toc;
dev = 0;
for i = 1:size(POAtable1, 1)
for j = i:size(POAtable1, 2)
dev = max(max(max(abs(POAtable1{i, j} - POAtable2{i, j}))));
dev = max(dev, max(max(abs(fluxRange1 - fluxRange2))));
end
end
fprintf('Maximum difference between the two solutions: %.4e\n', dev);
fprintf('\nSingle-thread computation: %.0f sec\nTwo-thread computation: %.0f sec\n', t4, t3);
%%
% The advantage will be more significant for more targets to analyzed and
% more threads used. Similar to |SteadyComFVA|, |SteadyComPOA| also supports continuation
% from previously interrupted computation by calling with the same |options.savePOA|.
%
%
%% REFERENCES
% _[1] _Chan SHJ, Simons MN, Maranas CD (2017) SteadyCom: Predicting microbial
% abundances while ensuring community stability. PLoS Comput Biol 13(5): e1005539.
% https://doi.org/10.1371/journal.pcbi.1005539
%
% _[2] _Khandelwal RA, Olivier BG, Röling WFM, Teusink B, Bruggeman FJ (2013)
% Community Flux Balance Analysis for Microbial Consortia at Balanced Growth.
% PLoS ONE 8(5): e64567. https://doi.org/10.1371/journal.pone.0064567
%
% __
4. Results display
6. Code Comments
%clean up environment variables
% clc
% clear
%initial environment(Execute after only one run, update link)
%initCobraToolbox;
%read model
% When testing the model, first build a community model for each model separately, modify the representation of biomass and save it again
modelFileName_1='AF.xls';
modelFileName_2='AP.xls';
modelFileName_3='AT.xls';
modelFileName_4='LB.xls';
models1= readCbModel(modelFileName_1);
models2= readCbModel(modelFileName_2);
models3= readCbModel(modelFileName_3);
models4= readCbModel(modelFileName_4);
models{1,1} = models1;
models{2,1} = models2;
models{3,1} = models3;
models{4,1} = models4;
%test
XP1=models1;
JY4=models4;
modelss{1,1} = XP1;
modelss{2,1} = JY4;
%The identification of biomass responses in the model is required to be bio and the reaction contains key metabolites. Each metabolite can be combined
%Building a multibacterial community model
[modelJoint] = createMultipleSpeciesModel(models);
%Obtain community metabolism information
nameTagsModels={'model1','model2','model3','model4'};%Species Response Name Prefix
[modelJoint.infoCom,modelJoint.indCom]=getMultiSpeciesModelId(modelJoint,nameTagsModels);
%Marker positions for biomass conversion Biomass production calculations
%spBm=[990,1834];
% spBm=zeros(size(nameTagsModels));
for i=1:size(nameTagsModels,2)
spBm(1,i)=find(strncmp(modelJoint.rxns, strcat(nameTagsModels(i),'_biomass'), 10));
spBm1(1,i)=strcat(nameTagsModels(i),'_biomass')
end
modelJoint.indCom.spBm=spBm;
modelJoint.infoCom.spBm=spBm1;
% modelJoint.infoCom=modelJoint.indCom;
%Analyzing Community Model Interactions
%Change all substance caps
modelJoint.ub(modelJoint.indCom.EXsp(:)~=0) = 1000;
%Display model upper and lower limits
printUptakeBoundCom(modelJoint, 1);
% Options
%Growth rate fitting, the growth rate of the first bacteria is 0.1 The second fungus is 0.2
% options.GRfx = [1,0.1;2,0.3];
% %determine whether to grow
% options.feasCrit=1;
% options.GRguess=[0.2,0.8];
options = struct();
options.GRguess = 0.1; % initial guess for max. growth rate initial growth rate
options.GRtol = 1e-6; % tolerance for final growth rate
options.algorithm = 1; % use the default algorithm (simple guessing for bounds, followed by matlab fzero)
options.algorithm = 2; % use the simple guessing algorithm
options.guessMethod=1;
% options.BMcon=[100;100;100;50];
% options.BMgdw=[1;2;1;2];
[sol, result, LP, LPminNorm, indLP] =SteadyCom(modelJoint,options);
[sol, result, LP, LPminNorm, indLP] =SteadyCom(modelJoint);
%%
options.optBMpercent = 100;
n = size(modelJoint.S, 2); % number of reactions in the model
% options.rxnNameList is the list of reactions subject to FVA. Can be reaction names or indices.
% Use n + j for the biomass variable of the j-th organism. Alternatively, use {'X_j'}
% for biomass variable of the j-th organism or {'X_Ec1'} for Ec1 (the abbreviation in EcCom.infoCom.spAbbr)
options.rxnNameList = spBm1;
options.optGRpercent = [1:10:99, 0:10:100]; % perform FVA at various percentages of the maximum growth rate, 89, 89.1, 89.2, ..., 100
[fvaComMin,fvaComMax] = SteadyComFVA(modelJoint, options);
[minFlux, maxFlux, minFD, maxFD, GRvector, result, LP] = SteadyComFVA(modelJoint);
%%
%Modify the objective function
modelJoint = changeObjective(modelJoint,{'model1_biomass','model2_biomass','model3_biomass'},3); %product of solution
%sol = optimizeCbModel(modelCom,'max','one');%To solve the maximum conversion rate of the target product, set the solving parameters,
[POAtable, fluxRange, Stat, GRvector] = SteadyComPOA(modelJoint);
% [sol, result, LP, LPminNorm, indLP] = SteadyComCplex(modelJoint);
%% Test a single model
modelCom = changeObjective(modelJoint,{'model1_biomass','model2_biomass','model3_biomass'},3); %The product of the solution, set the target product to be observed, M11_4 number the substrate
modelCom = changeObjective(modelJoint,{'model1_biomass'},1);
sol = optimizeCbModel(modelJoint,'max','one');%To solve the maximum conversion rate of the target product, set the solving parameters,
sol.obj=sol.f;