-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulationCrabtree.m
63 lines (49 loc) · 1.78 KB
/
simulationCrabtree.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
%% simulationCrabtree
% Timing: ~ 700 s
tic;
soplexpath = '/Users/cheyu/build/bin/soplex'; % change this to the soplex path on your PC
% load model
load('CofactorYeast.mat');
load('enzymedata.mat');
%% 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);
model = changeRxnBounds(model,'r_1631',0,'b');% acetaldehyde production
%% Set optimization
rxnID = 'r_1714'; %minimize glucose uptake rate
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
mu_list = [0.02:0.02:0.36 0.379];
fluxes = zeros(length(model.rxns),length(mu_list));
for i = 1:length(mu_list)
mu = mu_list(i);
model_tmp = changeRxnBounds(model,'r_2111',mu,'b');
disp(['mu = ' num2str(mu)]);
fileName = writeLP(model_tmp,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');
[sol_obj,sol_status,sol_full] = readSoplexResult('Simulation.lp.out',model_tmp);
disp(['solution status: ' sol_status]);
if strcmp(sol_status,'optimal')
fluxes(:,i) = sol_full;
end
end
cd Results/;
save('sC_fluxes.mat','fluxes');
cd ../;
clear;
toc;