Skip to content

Commit e183067

Browse files
IVANDOMENZAINIVANDOMENZAIN
committed
feat: add setHamsMedium function
1 parent c57de43 commit e183067

File tree

1 file changed

+142
-0
lines changed

1 file changed

+142
-0
lines changed

‎ComplementaryScripts/Simulation/setHamsMedium.m‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
function exchModel = setHamsMedium(model,Csource,irrev,measuredMets,fluxes)
2+
% setHamsMedium
3+
%
4+
% Set a Hams culture medium for a humanGEM based model. This function works
5+
% for either standard or enzyme constrained-GEMs
6+
%
7+
% model An ihuman-based GEM
8+
% Csource (string) metName for the main carbon source
9+
% irrev (logical) TRUE if model comes in an irreversible format.
10+
% Default = false
11+
% measuredMets (cell) met IDs for measured compounds. Optional
12+
% fluxes (vector) Measured fluxes [mmol/gDw h], always as positive
13+
% quantities. Optional
14+
%
15+
% exchModel (struct) Model with Ham's media constraints
16+
%
17+
% Ivan Domenzain. Last edited: 2019-11-29
18+
19+
if nargin<5
20+
fluxes = [];
21+
if nargin<4
22+
measuredMets = [];
23+
if nargin<3
24+
irrev = false;
25+
end
26+
end
27+
end
28+
%Remove unconstrained field, if available
29+
if isfield(model,'unconstrained')
30+
model = rmfield(model,'unconstrained');
31+
end
32+
%Check if boundary metabolites are present, if so then remove them
33+
boundaryIndx = find(strcmpi(model.compNames,'Boundary'));
34+
boundary = find(model.metComps==boundaryIndx);
35+
if ~isempty(boundary)
36+
model = removeMets(model,boundary,false,false,false,true);
37+
end
38+
%Locate carbon source uptake reaction
39+
Cindxs = find(strcmpi(model.metNames,Csource));
40+
Cindxs = Cindxs(find(model.metComps(Cindxs)==1));
41+
Csource = model.mets(Cindxs);
42+
%Ham's media composition
43+
mediaComps ={ Csource{1} ...
44+
'm01365s' ... %arginine[s]
45+
'm02125s' ... %histidine[s]
46+
'm02471s' ... %methionine[s]
47+
'm02724s' ... phenylalanine[s]
48+
'm03089s' ... tryptophan[s]
49+
'm03101s' ... tyrosine[s]
50+
'm01307s' ... alanine[s]
51+
'm01986s' ... glycine[s]
52+
'm02896s' ... serine[s]
53+
'm02993s' ... threonine[s]
54+
'm01370s' ... aspartate[s]
55+
'm01974s' ... glutamate[s]
56+
'm01369s' ... asparagine[s]
57+
'm01975s' ... glutamine[s]
58+
'm02184s' ... isoleucine[s]
59+
'm02360s' ... leucine[s]
60+
'm02770s' ... proline[s]
61+
'm03135s' ... valine[s]
62+
'm01628s' ... cysteine[s]
63+
'm02982s' ... thiamin[s]
64+
'm02159s' ... hypoxanthine[s]
65+
'm01830s' ... folate[s]
66+
'm01401s' ... biotin[s]
67+
'm02680s' ... pantothenate[s]
68+
'm01513s' ... choline[s]
69+
'm02171s' ... inositol[s]
70+
'm02583s' ... nicotinamide[s]
71+
'm02817s' ... pyridoxine[s]
72+
'm02842s' ... riboflavin[s]
73+
'm02996s' ... %thymidine[s]
74+
'm02630s' ... %Oxygen
75+
'm02040s' ... %Water
76+
'm02519s' ... %sodium
77+
'm02200s' ... %K+
78+
'm02751s' ... %Pi
79+
'm02946s' ... %Sulfate
80+
'm01413s' ... % calcium
81+
'm01821s' ... %Fe2+
82+
'm01822s' ... %Fe3+
83+
'm02046s' ... %HCO3-
84+
'm01450s' ... %cholesterol
85+
'm01596s' ... %CO2
86+
'm02039s'}; %H+
87+
%Default flux bounds
88+
fluxBounds = ones(1,length(mediaComps))*1000;
89+
%Check if provided mets are part of media's formulation
90+
if ~isempty(measuredMets)
91+
[iA,iB] = ismember(measuredMets,mediaComps);
92+
else
93+
iA = 0;
94+
end
95+
96+
if ~irrev
97+
%Modify fluxBounds with the correspondent provided flux measurements
98+
if sum(iA)>0
99+
fluxBounds(iB(iA)) = fluxes(iA);
100+
end
101+
%Set the right flux direction
102+
fluxBounds = -1*fluxBounds;
103+
%Allow all exchanges (Makes sure that all secretions are open
104+
[exchModel,~] = setExchangeBounds(model);
105+
%Close uptakes for mets not present in the media
106+
[exchModel,~] = setExchangeBounds(exchModel,mediaComps,fluxBounds,1000,true);
107+
else
108+
if sum(iA)>0
109+
fluxBounds(iB(iA)) = fluxes(iA);
110+
end
111+
[exchRxns,exchIndxs] = getExchangeRxns(model);
112+
%Exclude protein pool exchange
113+
exchIndxs = exchIndxs(1:end-1);
114+
exchRxns = exchRxns(1:end-1);
115+
%Differentiate between uptakes and production reactions
116+
uptkIndxs = exchIndxs(find(contains(exchRxns,'_REV')));
117+
prodIndxs = exchIndxs(find(~contains(exchRxns,'_REV')));
118+
%Open all production reactions
119+
exchModel = setParam(model,'ub',prodIndxs,1000);
120+
%close all uptakes
121+
exchModel = setParam(exchModel,'ub',uptkIndxs,0);
122+
%Open uptake of media components one by one
123+
for i=1:length(mediaComps)
124+
%Get metabolite indx
125+
metIndx = find(strcmpi(model.mets,mediaComps{i}));
126+
%Get rxns for metabolite
127+
metRxns = find(model.S(metIndx,:));
128+
%Get the uptake reaction for the metabolite
129+
metExchRxn = intersect(metRxns,uptkIndxs);
130+
%Open it!
131+
exchModel.ub(metExchRxn) = fluxBounds(i);
132+
end
133+
end
134+
%Check if model is feasible
135+
sol = solveLP(exchModel);
136+
if ~isempty(sol.x)
137+
disp('Constrained model is feasible')
138+
else
139+
disp('Constrained model is unfeasible')
140+
exchModel = [];
141+
end
142+
end

0 commit comments

Comments
 (0)