/
single_level_ko.gms
executable file
·336 lines (247 loc) · 6.37 KB
/
single_level_ko.gms
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
*************** OptForce for FA synthesis in E. coli ****************
*** OptForce set for all chain lengths
****************************************************************************
$INLINECOM /* */
$set myroot1 input_files/iAF1260/iAF1260
$set myroot2 input_files/
$set pdt ac
$set value 154.285
$set max_val 25.249
$set biom 0.09305
$set condition aero
$set reg with
options
limrow = 10000
limcol = 10000
optCR = 1E-9
optCA = 0.0
iterlim = 100000
decimals = 8
reslim = 100000
work = 50000000
sysout = off
solprint = on;
*******************************************
Sets
*** Metabolites
i
$include "%myroot1%_cmp.txt"
*** Reactions
j
$include "%myroot1%_rxnnames.txt"
*** fluxes providing nutrients
exchange(j)
$include "%myroot1%_source_M9.txt"
*** Reactions shut off during anaerobic glucose condition
anaerobic(j)
$include "%myroot1%_offanaero.txt"
*** Reactions shut off during aerobic glucose condition
offaeroglu(j)
$include "%myroot1%_offaeroglu.txt"
*** Rxns encoded by invivo essential genes
essential(j)
$include "%myroot2%Find_essential_%condition%_%reg%.txt"
*** Product
biom(j) / 'Ec_biomass_iAF1260_WT_59p81M' /
*** rxns with no gene associated
nogene(j)
$include "%myroot1%_nogene_rxns.txt"
*** Exchange rxns
excg(j)
$include "%myroot1%_ex_pp_rxns.txt"
$ONEMPTY
*** Blocked rxns
blocked(j)
$include "%myroot2%Find_blocked_%condition%_%reg%.txt"
intervention(j)
index /1*1000/
;
*****************************************************************************************
Parameters
*** stoichiometry matrix
s(i,j)
$include "%myroot1%_sij.txt"
*** rxn type
rxntype(j)
$include "%myroot1%_rxntype.txt"
*** flux bounds
LB(j) , UB(j)
*** matrix to store Interventions
matrix(index,j)
;
matrix(index,j) = 0;
********************************
Variables
*** objectives
z
zprimal
zdual
*** flux values
v(j)
*** duals for Sij
lambda(i)
*** for linearizing y(j)*delta(j)
w1l(j), w1u(j)
;
positive variables
*** duals for bounds
delta1l(j), delta1u(j)
;
binary variables
y0(j)
;
**************** INITIALIZING PARAMETRIC VARIABLES AND SETS ****************************
*** Scalars
scalar counter /1/;
scalar M /2000/;
scalar c /0.0001/;
scalar k /0/;
scalar n /0/;
****************************************************************************************
Equations
primal1
primal7
primal8
primal
dual
dual1
dual2
outer
outer1
outer3
outer4
outer5
outer6
outer7
outer8
outer9
outer10
outer11
outer24
;
*** outer
outer.. z =e= 0.035*v('EX_%pdt%(e)') + 0.965*v('Ec_biomass_iAF1260_WT_59p81M');
outer1.. sum(j, y0(j)) =e= k;
outer3.. v('Ec_biomass_iAF1260_WT_59p81M') =e= sum(j, delta1u(j)*UB(j) - w1u(j)*UB(j)) - sum(j, delta1l(j)*LB(j) - w1l(j)*LB(j));
outer24(index).. sum(j, matrix(index,j) * y0(j)) =l= k-1;
outer4(j).. w1l(j) =l= M*y0(j);
outer5(j).. w1l(j) =g= -M*y0(j);
outer6(j).. w1l(j) =l= delta1l(j) + M*(1-y0(j) );
outer7(j).. w1l(j) =g= delta1l(j) - M*(1-y0(j) );
outer8(j).. w1u(j) =l= M*y0(j);
outer9(j).. w1u(j) =g= -M*y0(j);
outer10(j).. w1u(j) =l= delta1u(j) + M*(1-y0(j) );
outer11(j).. w1u(j) =g= delta1u(j) - M*(1-y0(j) );
*** primal
primal.. zprimal =e= v('Ec_biomass_iAF1260_WT_59p81M');
primal1(i).. sum(j, (S(i,j)*v(j))) =e= 0;
primal7(j).. v(j) =g= LB(j)*(1-y0(j) );
primal8(j).. v(j) =l= UB(j)*(1-y0(j) );
*** dual
dual.. zdual =e= sum(j, delta1u(j)*UB(j) - delta1u(j)*y0(j)*UB(j)) - sum(j, delta1l(j)*LB(j) - delta1l(j)*y0(j)*LB(j));
dual1(j)$biom(j).. sum(i,lambda(i)*S(i,j)) + delta1u(j) - delta1l(j) =e= 1;
dual2(j)$(not biom(j)).. sum(i,lambda(i)*S(i,j)) + delta1u(j) - delta1l(j) =e= 0;
***************************************************************
*** Get bounds for Fluxes
scalar max /1000/;
LB(j)$(rxntype(j) = 0) = 0;
UB(j)$(rxntype(j) = 0) = max;
LB(j)$(rxntype(j) = 2) = 0;
UB(j)$(rxntype(j) = 2) = 0;
LB(j)$(rxntype(j) = 1) = -max;
UB(j)$(rxntype(j) = 1) = max;
LB(j)$(rxntype(j) = 4) = 0;
UB(j)$(rxntype(j) = 4) = max;
LB(j)$(exchange(j)) = -max;
*** setting core biomass flux to 0
LB('Ec_biomass_iAF1260_core_59p81M') = 0;
UB('Ec_biomass_iAF1260_core_59p81M') = 0;
*** setting LB wild biomass flux to 10% of theoretical maximum
LB('Ec_biomass_iAF1260_WT_59p81M') = %biom%;
***shutting off reactions in aerobic glucose conditions
*LB(j)$anaerobic(j) = 0;
*UB(j)$anaerobic(j) = 0;
LB(j)$offaeroglu(j) = 0;
UB(j)$offaeroglu(j) = 0;
*** putting ATP Maintenance value
LB('ATPM') = 8.39;
UB('ATPM') = 8.39;
*** Glucose and oxygen uptake
LB('Ex_glc(e)') = -10;
LB('EX_o2(e)') = -20;
*** fix binaries
y0.fx(j)$essential(j) = 0;
y0.fx(j)$blocked(j) = 0;
y0.fx(j)$exchange(j) = 0;
y0.fx(j)$nogene(j) = 0;
y0.fx(j)$excg(j) = 0;
y0.fx('ATPM') = 0;
y0.fx('Ec_biomass_iAF1260_WT_59p81M') = 0;
*** Fixing certain interventions
y0.fx('PGCD') = 0;
y0.fx('PSERT') = 0;
y0.fx('ETOHtex_f') = 0;
********************************
model bilevel
/
primal1
primal7
primal8
*dual1
*dual2
outer
outer1
*outer3
*outer4
*outer5
*outer6
*outer7
*outer8
*outer9
*outer10
outer24
/
;
options iterlim = 1000000;
bilevel.optfile = 1;
scalar cutoff /0/;
scalar cutoff1 /0/;
scalar yield;
file forced /sing_level_ko.txt/;
put forced;
put '****************************************'//;
for (k = 2 to 2 by 1,
counter = 1;
n = 0;
matrix(index,j) = 0;
put 'Number of Interventions (k) =' k//;
put '****************************************'//;
while( counter = 1,
solve bilevel using mip maximizing z;
n = n + 1;
if(n eq 1, cutoff = z.l);
if(bilevel.modelstat eq 1,
if( z.l <= cutoff1 + .001,
counter = 0;
put "Objective value/better solution not found"/;
cutoff1 = cutoff;
);
put 'Iteration No. =' n//;
put 'v(%pdt%) = 'v.l('EX_%pdt%(e)'):0:8/;
put 'v(biomass) = 'v.l('Ec_biomass_iAF1260_WT_59p81M'):0:8/;
put /'knockouts :'//;
loop (j$(y0.l(j) = 1),
put "'"j.tl:0"' "v.l(j):0:8/;
);
put /'********************'//;
matrix(index, j)$(ord(index) = n and (y0.l(j) = 1)) = 1;
);
if(bilevel.modelstat ne 1 or n >= 6,
counter = 0;
put "No feasible solution achieved/max n reached, modelstat : "bilevel.modelstat//;
cutoff1 = cutoff;
);
);
put /'****************************************'/;
);
putclose forced;