Skip to content

Commit fa2adac

Browse files
rmtflemingrmtfleming
authored
Merge pull request opencobra#2330 from Kretaks/optEnvelope
optEnvelope
2 parents c8acdba + 6cdc246 commit fa2adac

File tree

8 files changed

+1042
-0
lines changed

8 files changed

+1042
-0
lines changed

‎src/design/optEnvelope/addEnv.m‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
function line = addEnv(origModel, biomass, desiredProduct, varargin)
2+
% addEnv adds envelope to figure
3+
% Algorithm is able to knock out genes as well as reactions to produce
4+
% production envelope
5+
%
6+
% USAGE:
7+
% line = addEnv(origModel, biomass, desiredProduct, 'KnockOuts', knockouts, 'colour', colour, 'prodMol', prodMol, 'subUptake', subUptake, 'molarSum', molarSum)
8+
%
9+
% INPUTS:
10+
% origModel COBRA model structure [struct]
11+
% biomass Reaction name of biomass [char]
12+
% desiredProduct Reaction name of desired product [char]
13+
%
14+
% OPTIONAL INPUTS:
15+
% KnockOuts (opt) List of knockouts for production envelope [cell array] (default: {})
16+
% colour (opt) Short name for colour of line to plot [anything that can be colour in matlab] (default: 'r' (red))
17+
% prodMol (opt) Molar mass of target product for yield plot [double]
18+
% subUptake (opt) Uptake of substrate for yield plot [double]
19+
% molarSum (opt) Molar mass of substrate for yield plot [double]
20+
%
21+
% OUTPUTS
22+
% line Line data for plot function
23+
%
24+
% NOTES
25+
% Sometimes last point of envelope drops to zero (might be rounding error)
26+
% but this function connects last points of lines so the graph creates
27+
% continuous line.
28+
% This algorithm only adds graph. It does not change labels.
29+
%
30+
% EXAMPLE:
31+
% line = addEnv(model, 'BIOMASS_Ecoli', 'EX_ac_e', {'GHMT2r','GND','SUCCt2r','SUCD4','AKGt2r','GLUt2r'}, 'm')
32+
%
33+
% AUTHORS:
34+
% Created by Kristaps Berzins 31/10/2022
35+
% Modified by Kristaps Berzins 30/09/2024
36+
37+
parser = inputParser();
38+
parser.addRequired('model', @(x) isstruct(x) && isfield(x, 'S') && isfield(origModel, 'rxns')...
39+
&& isfield(origModel, 'mets') && isfield(origModel, 'lb') && isfield(origModel, 'ub') && isfield(origModel, 'b')...
40+
&& isfield(origModel, 'c'))
41+
parser.addRequired('biomass', @(x) any(validatestring(x, origModel.rxns)))
42+
parser.addRequired('desiredProduct', @(x) any(validatestring(x, origModel.rxns)))
43+
parser.addOptional('KnockOuts', {}, @(x) iscell(x) && ismatrix(x))
44+
parser.addOptional('colour', 'r', @(x) any(validatecolor(x)))
45+
parser.addOptional('prodMol', [], @(x) isnumeric(x))
46+
parser.addOptional('subUptake', 10, @(x) isnumeric(x))
47+
parser.addOptional('molarSum', 180, @(x) isnumeric(x))
48+
49+
parser.parse(origModel, biomass, desiredProduct, varargin{:});
50+
origModel = parser.Results.model;
51+
biomass = parser.Results.biomass;
52+
desiredProduct = parser.Results.desiredProduct;
53+
KnockOuts = parser.Results.KnockOuts;
54+
colour = parser.Results.colour;
55+
prodMol = parser.Results.prodMol;
56+
subUptake = parser.Results.subUptake;
57+
molarSum = parser.Results.molarSum;
58+
59+
if isempty(prodMol)
60+
prodMolIs = false;
61+
else
62+
prodMolIs = true;
63+
end
64+
65+
model = origModel;
66+
67+
if any(ismember(model.rxns, KnockOuts))
68+
rxns = ismember(model.rxns, KnockOuts);
69+
model.ub(rxns) = 0;
70+
model.lb(rxns) = 0;
71+
elseif any(ismember(model.genes, KnockOuts))
72+
model = buildRxnGeneMat(model);
73+
[model, ~, ~] = deleteModelGenes(model, KnockOuts);
74+
%elseif %Enzymes
75+
end
76+
77+
solMin = optimizeCbModel(model, 'min');
78+
solMax = optimizeCbModel(model, 'max');
79+
controlFlux1 = linspace(solMin.f, solMax.f, 100)';
80+
if nnz(controlFlux1) == 0
81+
return;
82+
end
83+
model = changeObjective(model, desiredProduct);
84+
85+
for i = 1:numel(controlFlux1)
86+
model = changeRxnBounds(model, biomass, controlFlux1(i), 'b');
87+
s = optimizeCbModel(model, 'min'); Min1(i, 1) = s.f;
88+
if s.stat == 0
89+
model = changeRxnBounds(model, biomass, controlFlux1(i) - 0.0001 * controlFlux1(i), 'b');
90+
s = optimizeCbModel(model, 'min'); Min1(i, 1) = s.f;
91+
s = optimizeCbModel(model, 'max'); Max1(i, 1) = s.f;
92+
end
93+
s = optimizeCbModel(model, 'max'); Max1(i, 1) = s.f;
94+
if s.stat == 0
95+
model = changeRxnBounds(model, biomass, controlFlux1(i) - 0.0001 * controlFlux1(i), 'b');
96+
s= optimizeCbModel(model,'min');Min1(i,1)=s.f;
97+
s= optimizeCbModel(model,'max');Max1(i,1)=s.f;
98+
end
99+
end
100+
101+
if prodMolIs
102+
controlFlux1 = controlFlux1 / subUptake * 1000 / molarSum;
103+
Max1 = Max1 / molarSum * prodMol / subUptake;
104+
Min1 = Min1 / molarSum * prodMol / subUptake;
105+
end
106+
107+
hold on
108+
line = plot(controlFlux1, Max1, 'color', colour, 'LineWidth', 2);
109+
plot(controlFlux1, Min1, 'color', colour, 'LineWidth', 2)
110+
hold off
111+

‎src/design/optEnvelope/milpOEReinserts.m‎

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
function [knockouts] = milpOEReinserts(model, data, K, minP, numKO, toDel, timeLimit, printLevel)
2+
% This function is creating and calculating MILP to get certain amount
3+
% of knockouts to achieve best production envelope
4+
%
5+
% USAGE:
6+
% [knockouts] = milpOEReinserts(model, data, K, minP, numKO, toDel, timeLimit, printLevel)
7+
%
8+
% INPUTS:
9+
% model COBRA model structure in irreversible form [struct]
10+
% data Struct with information about:
11+
% * mainActive List of active reactions for main envelope [cell array]
12+
% K List of reactions that cannot be selected for knockout (reaction IDs) [double array]
13+
% minP Struct with information about biomass and desired product.
14+
% * bioID ID of biomass [double]
15+
% * proID ID of desired product[double]
16+
% numKO Number of knockouts to achieve [double]
17+
%
18+
% OPTIONAL INPUTS:
19+
% toDel Numeric variable that shows what to delete:
20+
% 0: reactions
21+
% 1: genes
22+
% 2: enzymes
23+
% timeLimit Time limit for gurobi optimization (in seconds) [double]
24+
% printLevel Print level for gurobi optimization [double]
25+
%
26+
% OUTPUTS:
27+
% knockouts List of reactions that when removed gives optimal envelope
28+
%
29+
% AUTHORS:
30+
% created by Kristaps Berzins 31/10/2022
31+
%
32+
% NOTES:
33+
% This function is not designed for stand-alone use. Should be used by
34+
% using optEnvelope.m with set numKO parameter
35+
36+
if nargin < 6
37+
toDel = 0;
38+
end
39+
if nargin < 7
40+
timeLimit = inf;
41+
end
42+
if nargin < 8
43+
printLevel = 0;
44+
end
45+
%%
46+
switch toDel
47+
case 0
48+
rxns = findRxnIDs(model, data.mainActive);
49+
rxns = ismember(model.rxns,[model.rxns(K);model.rxns(rxns)]);
50+
reactions = model.rxns(~rxns);
51+
52+
% MILP
53+
54+
[r,c] = size(model.S);
55+
56+
model.C_chemical=zeros(c,1);
57+
model.C_chemical(minP.proID)=1;
58+
59+
yInd=[];
60+
notYInd=[];
61+
reactionIDs = findRxnIDs(model, reactions);
62+
63+
for i=1:c
64+
if find(reactionIDs == i)
65+
yInd(end+1)=i;
66+
continue;
67+
else
68+
notYInd(end+1)=i;
69+
end
70+
end
71+
72+
notYInd=sort(notYInd)';
73+
yInd=sort(yInd)';
74+
75+
ySize=size(yInd,1);
76+
I = eye(c);
77+
A=[model.S;-model.S;I(notYInd,:);-I(notYInd,:); I(yInd,:); -I(yInd,:)];
78+
79+
[aSizeRow, vSize]=size(A);
80+
reinserts=zeros(c,1);
81+
reinserts(yInd)=1;
82+
83+
Ay1=diag(reinserts); %ub
84+
Ay1=Ay1*diag(model.ub);
85+
86+
Ay2=diag(reinserts); %lb
87+
Ay2=Ay2*diag(model.lb);
88+
89+
z1=find(Ay1);
90+
z2=find(Ay2);
91+
zSize=size([z1;z2],1);
92+
93+
Ay= [zeros(2*r+2*(vSize-ySize),ySize);
94+
-Ay1(yInd,yInd);
95+
Ay2(yInd,yInd);
96+
];
97+
98+
B = [zeros(2*r,1);
99+
model.ub(notYInd);
100+
-model.lb(notYInd);
101+
zeros(2*(ySize),1);
102+
];
103+
104+
C=model.c;
105+
106+
[A_w,Ay_w ,B_w,C_w, ~, ~, wSize] = seperateTransposeJoinOE(-A,-Ay,-B,-C,ySize,1,c,1000,zSize); %max(model.ub) 1000
107+
108+
awSizeRow=size(A_w,1);
109+
110+
Cjoined=[model.C_chemical; zeros(wSize+zSize,1); zeros(ySize,1)];
111+
112+
A2=-[
113+
-C', -C_w';
114+
A, sparse(aSizeRow, wSize+zSize);
115+
sparse(awSizeRow, vSize), A_w];
116+
117+
Ay2=-[
118+
zeros(1, ySize);
119+
Ay;
120+
Ay_w];
121+
122+
C2=Cjoined(1:vSize+wSize+zSize,:);
123+
124+
B2=-[
125+
0;
126+
B;
127+
B_w;
128+
];
129+
130+
z3=find(Ay2);
131+
zSizeOptKnock2=size(z3,1);
132+
133+
[A2_w,Ay2_w,B2_w,C2_w,lb2_w,ub2_w,uSize]=seperateTransposeJoinOE(A2,Ay2,B2,C2,ySize,1,vSize+wSize+zSize,1000,zSizeOptKnock2); %max(model.ub) 1000
134+
135+
[A2_wRow, A2_wCol]=size(A2_w);
136+
[ARow, ACol]=size(A);
137+
138+
A3=[
139+
A2_w, sparse(A2_wRow,ACol), Ay2_w; %dual constraints
140+
zeros(1,uSize+zSizeOptKnock2+vSize), -ones(1,ySize); %y sum constraints
141+
sparse(ARow,uSize+zSizeOptKnock2), A, Ay; %feasibility conatraint
142+
];
143+
144+
B3=[
145+
B2_w;
146+
numKO-ySize;
147+
B;
148+
];
149+
150+
C3=[C2_w;
151+
zeros(ACol,1);
152+
zeros(ySize,1);
153+
];
154+
155+
tmpL=lb2_w(1:uSize+zSizeOptKnock2);
156+
tmpH=ub2_w(1:uSize+zSizeOptKnock2);
157+
ysUpperBound=ones(ySize,1);
158+
lb3=[tmpL; model.lb; zeros(ySize,1)];
159+
ub3=[tmpH; model.ub; ysUpperBound];
160+
intVars=A2_wCol+ACol+1:A2_wCol+ACol+ySize;
161+
162+
MILP.c = C3;
163+
MILP.osense = -1; % max
164+
MILP.A = sparse(A3);
165+
MILP.b = B3;
166+
MILP.lb = lb3; MILP.ub = ub3;
167+
MILP.x0 = [];
168+
MILP.vartype = char(ones(1,length(C3)).*double('C'));
169+
MILP.vartype(intVars) = 'B';
170+
MILP.csense = char(ones(1,length(B3)).*double('L'));
171+
172+
milpSol = solveCobraMILP(MILP, 'timeLimit', timeLimit, 'printLevel', printLevel);
173+
174+
%% add only found reactions
175+
176+
rxns = reactions(milpSol.int == 1);
177+
knockouts = reactions(milpSol.int ~= 1);
178+
numel(rxns)
179+
numel(knockouts)
180+
%debug
181+
[irrev,~,~,~] = convertToIrreversible(model);
182+
addEnv(irrev, irrev.rxns(minP.bioID), irrev.rxns(minP.proID), knockouts, 'r');
183+
print('debug');
184+
case 1
185+
%Genes
186+
case 2
187+
%Enzymes
188+
end

0 commit comments

Comments
 (0)