-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulationIronFVA.m
70 lines (54 loc) · 2.12 KB
/
simulationIronFVA.m
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
%% simulationIronFVA
%
% Timing: ~ 1800 s
tic;
load('CofactorYeast.mat');
load('enzymedata.mat');
soplexpath = '/Users/cheyu/build/bin/soplex'; % change this to the soplex path on your PC
%% Set model
% set medium
model = setMedia(model,1);% minimal media (Delft media)
% set carbon source
model = changeRxnBounds(model,'r_1714',-1000,'l');% glucose
% set oxygen
model = changeRxnBounds(model,'r_1992',-1000,'l');
% block reactions
model = blockRxns(model);
%% Set optimization
rxnID = 'dilute_dummy';
osenseStr = 'Maximize';
tot_protein = 0.46; %g/gCDW, estimated from the original GEM.
f_modeled_protein = extractModeledprotein(model,'r_4041','s_3717[c]'); %g/gProtein
% r_4041 is pseudo_biomass_rxn_id in the GEM
% s_3717[c] is protein id
f = tot_protein * f_modeled_protein;
f_mito = 0.1;
clear tot_protein f_modeled_protein;
factor_k_withoutcofator = 0;
%% Solve LPs
% reference
[mu_ref,flux_ref] = searchMaxgrowth(model,f,f_mito,osenseStr,rxnID,enzymedata,factor_k_withoutcofator,1e-6,soplexpath);
q_fe_ref = flux_ref(strcmp(model.rxns,'r_1861'),1);
%% FVA
rxnID = 'r_1861';
mu = mu_ref;
model = changeRxnBounds(model,'r_2111',mu,'b');
osenseStr = 'Maximize';
disp(osenseStr);
fileName = writeLP(model,mu,f,f_mito,osenseStr,rxnID,enzymedata,factor_k_withoutcofator);
command = sprintf([soplexpath,' -s0 -g5 -t300 -f1e-20 -o1e-20 -x -q -c --int:readmode=1 --int:solvemode=2 --int:checkmode=2 --real:fpfeastol=1e-9 --real:fpopttol=1e-9 %s > %s.out %s'],fileName,fileName);
system(command,'-echo');
[q_fe_max,~,flux_max] = readSoplexResult('Simulation.lp.out',model);
osenseStr = 'Minimize';
disp(osenseStr);
fileName = writeLP(model,mu,f,f_mito,osenseStr,rxnID,enzymedata,factor_k_withoutcofator);
command = sprintf([soplexpath,' -s0 -g5 -t300 -f1e-20 -o1e-20 -x -q -c --int:readmode=1 --int:solvemode=2 --int:checkmode=2 --real:fpfeastol=1e-9 --real:fpopttol=1e-9 %s > %s.out %s'],fileName,fileName);
system(command,'-echo');
[q_fe_min,~,flux_min] = readSoplexResult('Simulation.lp.out',model);
sIFVA_res.flux_ref = flux_ref;
sIFVA_res.flux_max = flux_max;
sIFVA_res.flux_min = flux_min;
cd Results/;
save('sIFVA_res.mat','sIFVA_res');
cd ../;
toc;