-
Notifications
You must be signed in to change notification settings - Fork 2
/
fluxVariability_new.m
802 lines (733 loc) · 35.5 KB
/
fluxVariability_new.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
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
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
function [minFlux, maxFlux, Vmin, Vmax] = fluxVariability_new(model, varargin)
% Performs flux variablity analysis
%
% USAGE:
%
% [minFlux, maxFlux] = fluxVariability(model, optPercentage, osenseStr, rxnNameList, printLevel, allowLoops)
% [minFlux, maxFlux, Vmin, Vmax] = fluxVariability(model, optPercentage, osenseStr, rxnNameList, printLevel, allowLoops, method, solverParams, advind, threads, heuristics, useMtFVA)
% [...] = fluxVariability(model, ..., 'name', value, ..., solverParams)
% [...] = fluxVariability(model, ..., paramStruct)
%
% INPUT:
% model: COBRA model structure
%
% OPTIONAL INPUTS:
% (support name-value argument inputs or a single [function + solver] parameter structure input)
% optPercentage: Only consider solutions that give you at least a certain
% percentage of the optimal solution (Default = 100
% or optimal solutions only)
% osenseStr: Objective sense, 'min' or 'max' (default)
% rxnNameList: List of reactions for which FVA is performed
% (Default = all reactions in the model)
% printLevel: Verbose level (default: 0). 1 to show a progress bar. 2 to print results for each reaction
% allowLoops: Whether loops are allowed in solution or which method to block loops.
%
% * 1 (or true) : loops allowed (default)
% * 0 (or false): loops not allowed. Use to find loopless solutions
% * 'original' : original loopless FVA
% * 'fastSNP' : loopless FVA with with Fast-SNP preprocessing of nullspace
% * 'LLC-NS' : localized loopless FVA using information from nullsapce
% * 'LLC-EFM' : localized loopless FVA using information from EFMs.
% Require CalculateFluxModes.m from EFMtool to calculate EFMs.
% method: when Vmin and Vmax are in the output, the flux vector can be (Default = 2-norm):
%
% * 'FBA' : standards FBA solution
% * '0-norm' : minimzes the vector 0-norm
% * '1-norm' : minimizes the vector 1-norm
% * '2-norm' : minimizes the vector 2-norm
% * 'minOrigSol' : minimizes the euclidean distance of each vector to the original solution vector
%
% solverParams: solver-specific parameter structure. Can also be inputted as the first or last arguement
% if using name-value argument inputs (with or without the keyword 'solverParams').
% Can also be inputted as part of a parameter structure together with other function parameters
%
% advind: switch to use the solution basis
%
% - 0 : default
% - 1 : uses the original problem solution basis as advanced basis
%
% threads: number of threads used for the analysis
% * 1, 2, 3, ...: number of threads
% * 0: defaulted number of threads for the parallel computing toolbox
% (default to be 1 if no parpool is activited, otherwise use the existing parpool)
%
% heuristics: level of heuristics to accelerate FVA.
% 0: no heuristics (default if rxnNameList has < 5 reactions)
% 1: solve max-sum-flux and min-sum-flux LPs to get reactions which already hit the bounds
% 2: solve additionally a single LP to find all blocked irreversible reactions (default if rxnNameList has >= 5 reactions)
%
% useMtFVA: run FVA multi-threaded via an external JVM with CPLEX as solver
% does not return Vmin, Vmax; requires allowLoops = true and method = 'FBA'
%
% - 0 : default, do not use mtFVA
% - 1 : use mtFVA
%
% paramStruct: one single parameter structure including any of the inputs above and the solver-specific parameter
%
% OUTPUTS:
% minFlux: Minimum flux for each reaction
% maxFlux: Maximum flux for each reaction
%
% OPTIONAL OUTPUT:
% Vmin: Matrix of column flux vectors, where each column is a separate minimization.
% Vmax: Matrix of column flux vectors, where each column is a separate maximization.
%
% EXAMPLES:
% FVA for all rxns at 100% max value of the objective function:
% [minFlux, maxFlux] = fluxVariability(model);
% Loopless FVA for rxns in `rxnNames` at <= 90% min value of the objective function, print results for each reaction:
% [minFlux, maxFlux] = fluxVariability(model, 90, 'min', rxnNames, 2, 0);
% Same as the 1st example, but also return the corresponding flux distributions with 2-norm minimized:
% [minFlux, maxFlux, Vmin, Vmax] = fluxVariability(model, [], [], [], 0, 1, '2-norm');
% Name-value inputs, with Cobra LP parameter `feasTol` and solver-specific (gurobi) parameter `Presolve`:
% [minFlux, maxFlux] = fluxVariability(model, 'optPercentage', 99, 'allowLoops', 0, 'threads', 0, 'feasTol', 1e-8, struct('Presolve', 0));
% Single parameter structure input including function, Cobra LP and solver parameters:
% [minFlux, maxFlux] = fluxVariability(model, struct('optPercentage', 99, 'allowLoops', 'original', 'threads', 0, 'feasTol', 1e-8, 'Presolve', 0));
%
% .. Authors:
% - Markus Herrgard 8/21/06 Original code.
% - Ronan Fleming 01/20/10 Take the extremal flux from the flux vector,
% not from the objective since this is invariant
% to the value and sign of the coefficient
% - Ronan Fleming 27/09/10 Vmin, Vmax
% - Marouen Ben Guebila 22/02/2017 Vmin,Vmax method
global CBT_LP_PARAMS
optArgin = { 'optPercentage', 'osenseStr', 'rxnNameList', 'printLevel', 'allowLoops', 'method', 'solverParams', 'advind', 'threads', 'heuristics', 'useMtFVA'};
defaultValues = {100, getObjectiveSense(model), model.rxns, 0, true, '2-norm', struct(), 0, [], [], 0};
validator = {@(x) isscalar(x) & isnumeric(x) & x >= 0 & x <= 100, ... % optPercentage
@(x) strcmp(x, 'max') | strcmp(x, 'min'), ... % osenseStr
@(x) ischar(x) | iscellstr(x), ... % rxnNameList
@(x) isscalar(x) & (islogical(x) | isnumeric(x)) & x >= 0, ... % printLevel
@(x) isscalar(x) | ischar(x), ... % allowLoops
@(x) ischar(x), ... % method
@isstruct, ... % solverParams
@(x) true, ... % advind
@(x) isscalar(x) & (islogical(x) | isnumeric(x)), ... % threads
@(x) isscalar(x) & (islogical(x) | isnumeric(x)) ... % heuristics
@(x) isscalar(x) & (islogical(x) | isnumeric(x)) ... % useMtFVA
};
% get all potentially supplied COBRA parameter names
problemTypes = {'LP', 'MILP', 'QP', 'MIQP'};
[funParams, cobraParams, solverVarargin] = parseCobraVarargin(varargin, optArgin, defaultValues, validator, problemTypes, 'solverParams', true);
if length(funParams)==10
% solverParams not output as a function parameter since it is individually handled and embedded in solverVarargin
[optPercentage, osenseStr, rxnNameList, printLevel, allowLoops, method, ...
advind, threads, heuristics, useMtFVA] = deal(funParams{:});
elseif length(funParams)==11
% default empty solverParams output as a function parameter
[optPercentage, osenseStr, rxnNameList, printLevel, allowLoops, method, ...
solverParams, advind, threads, heuristics, useMtFVA] = deal(funParams{:});
fields = fieldnames(solverParams);
if ~isempty(fields)
disp(fields)
warning('expecting default empty solverParams, but it is not empty')
end
else
error('Failure to parse inputs')
end
allowLoopsError = false;
loopMethod = '';
if isscalar(allowLoops) && ~ischar(allowLoops)
if allowLoops
loopMethod = 'none';
else
% default using LLCs when the method is not specified
loopMethod = 'LLC-NS';
end
elseif ischar(allowLoops)
for str = {'none', 'original', 'fastSNP', 'LLC-NS', 'LLC-EFM'}
if strncmpi(allowLoops, str{:}, length(allowLoops))
loopMethod = str{:};
break
end
end
allowLoops = strcmp(loopMethod, 'none');
if isempty(loopMethod)
allowLoopsError = true;
end
else
allowLoopsError = true;
end
if allowLoopsError
error('"allowLoops" must be one of the following: 1 (usual FVA), 0 (default using ''LLC-NS''), ''original'', ''fastSNP'', ''LLC-NS'' or ''LLC-EFM''')
end
if ischar(rxnNameList)
rxnNameList = {rxnNameList};
end
%Stop if there are reactions, which are not part of the model
if any(~ismember(rxnNameList,model.rxns))
presence = ismember(rxnNameList,model.rxns);
error('There were reactions in the rxnList which are not part of the model:\n%s\n',strjoin(rxnNameList(~presence),'\n'));
end
if useMtFVA && (nargout > 2 || ~allowLoops || ~strcmp(method,'FBA'))
error('mtFVA only supports the FBA method and neither supports loopless contraints nor Vmin/Vmax');
end
% Set up the problem size
[~, nRxns] = size(model.S);
% LP solution tolerance
[minNorm, tol] = deal(0, 1e-6);
if exist('CBT_LP_PARAMS', 'var') && isfield(CBT_LP_PARAMS, 'objTol')
tol = CBT_LP_PARAMS.objTol;
end
if nargout >= 3 && ~useMtFVA
minNorm = 1;
end
% Return if minNorm is minOrigSol but allowloops is set to false
if ~allowLoops && minNorm && strcmp(method,'minOrigSol')
error(['minOrigSol is meant for finding a minimally adjusted solution from an FBA solution. ', ...
'Cannot return solutions if allowLoops is set to false. ', ...
'If you want solutions without loops please set method to ''FBA'', ''2-norm'', ''1-norm'' or ''0-norm''.']);
end
% Determine constraints for the correct space (0-100% of the full space)
if sum(model.f ~= 0) > 0
hasObjective = true;
else
hasObjective = false;
end
if ~isfield(model,'b')
model.b = zeros(size(model.S,1),1);
end
% Set up the general problem
LPproblem = buildLPproblemFromModel(model);
loopInfo = struct('method', loopMethod, 'printLevel', printLevel);
% Solve initial (normal) LP
if allowLoops
tempSolution = solveCobraLP(LPproblem, solverVarargin.LP{:});
else
% Both Fast-SNP and solving an MILP return a minimal feasible nullspace
if printLevel
switch loopMethod
case 'original'
% find the usual internal nullspace (Schellenberger et al., 2009)
fprintf('Use loop law MILP implementation by Schellenberger et al., 2009\n');
case 'fastSNP'
% find a minimal feasible nullspace Fast-SNP (Saa and Nielson, 2016)
fprintf('Reduce complexity by nullspace preprocessing (Fast-SNP)\n')
otherwise
% find a minimal feasible nullspace by one single MILP (Chan et al., 2017)
% then implement localized loopless constraints
fprintf('Reduce complexity by nullspace preprocessing and implementing localized loopless constraints (LLCs)\n')
end
end
[MILPproblem, loopInfo] = addLoopLawConstraints(LPproblem, model, 1:nRxns, [], [], loopInfo);
if ~strncmp(loopMethod, 'LLC', 3)
tempSolution = solveCobraMILP(MILPproblem, solverVarargin.MILP{:});
else
% preprocessing for LLCs
[solveLP, MILPproblem, loopInfo] = processingLLCs('preprocess', loopInfo, LPproblem, model, nRxns, osenseStr, MILPproblem);
if solveLP
tempSolution = solveCobraLP(LPproblem, solverVarargin.LP{:});
else
tempSolution = solveCobraMILP(MILPproblem, solverVarargin.MILP{:});
end
end
end
if tempSolution.stat == 1
if strcmp(osenseStr,'max')
objValue = floor(tempSolution.obj / tol) * tol * optPercentage / 100;
else
objValue = ceil(tempSolution.obj / tol) * tol * optPercentage / 100;
end
else
error('The FVA could not be run because the model is infeasible or unbounded')
end
%set the objective
if hasObjective
LPproblem.A = [LPproblem.A; columnVector(LPproblem.f)'];
LPproblem.b = [LPproblem.b; objValue];
if strcmp(osenseStr, 'max')
LPproblem.csense(end+1) = 'G';
else
LPproblem.csense(end+1) = 'L';
end
model = addCOBRAConstraints(model, model.rxns(model.f ~= 0), objValue, 'dsense', LPproblem.csense(end));
end
% get the initial basis
if advind == 1
LPproblem.basis = tempSolution.basis;
end
LPproblem.S = LPproblem.A; % needed for sparse optimisation
%Thats not true. The Euclidean norm does not get rid of loops if the
%objective reaction is part of the loop.
% if length(minNorm)> 1 || minNorm > 0
% %minimizing the Euclidean norm gets rid of the loops, so there
% %is no need for a second slower MILP approach
% allowLoops=1;
% end
v = ver;
PCT = 'Parallel Computing Toolbox';
parpoolOn = false;
if any(strcmp(PCT, {v.Name})) && license('test', 'Distrib_Computing_Toolbox')
try
parpoolOn = ~isempty(gcp('nocreate'));
PCT_status = 1;
catch
PCT_status = 0;
end
else
PCT_status = 0; % Parallel Computing Toolbox not found.
end
parallelJob = false;
if isempty(threads) && PCT_status % no input for threads and parallel toolbox exists
% do parallel job if a parpool already exists, otherwise not
parallelJob = parpoolOn;
elseif ~isempty(threads) && threads ~= 1 && PCT_status % explicit input for parallel job
parallelJob = true;
if ~parpoolOn
% create a parpool if no existing one
if threads <= 0 % <= 0 for the default number of workers
parpool;
elseif threads > 1 % otherwise specified number of workers
parpool(threads);
end
end
end
minFlux = model.lb(findRxnIDs(model, rxnNameList));
maxFlux = model.ub(findRxnIDs(model, rxnNameList));
% preCompMaxSols(i) = k => the k-th heuristic solution solves max v_i
% k = 1: no heuristic solution
% k = 2: reactions that hit bounds during max-sum-flux heuristic
% k = 3: reactions that hit bounds during min-sum-flux heuristic
% k = 4: reactions that hit bounds during blocked-irreversible-reaction heuristic
% k = 5: blocked reactions found during blocked-irreversible-reaction heuristic
[preCompMaxSols, preCompMinSols] = deal(ones(numel(rxnNameList), 1));
heuristicSolutions = cell(5, 1);
[maxSolved, minSolved] = deal(false(numel(rxnNameList), 1));
[Vmax, Vmin] = deal([]);
if minNorm
[Vmax, Vmin] = deal(zeros(nRxns, numel(rxnNameList)));
end
% set the level of heuristics
if isempty(heuristics)
if numel(rxnNameList) >= 5
heuristics = 2;
else
heuristics = 0;
end
if printLevel > 0
fprintf(['> The level of heuristics has been set to ' num2str(heuristics) '.\n'])
end
end
% turn heuristics off for Mosek (don't solve the heuristic problem)
if strcmp(cobraParams.LP.solver, 'mosek')
heuristics = 0;
warning([' > The level of heuristics has been set to 0 for the Mosek solver']);
end
% each cell in rxnNameList must be a reaction at this point, otherwise there would be error earlier
Order = findRxnIDs(model, rxnNameList);
if heuristics > 0
%We will calculate a min and max sum flux solution.
%This solution will (hopefully) provide multiple solutions for individual
%reactions.
QuickProblem = LPproblem;
QuickProblem.f(:) = 0;
QuickProblem.f(Order) = 1;
if any(strcmp(loopMethod, {'original', 'fastSNP'}))
% Skip this when using localized loopless constraints (LLCs) for two reasons:
% 1. With loopless constraints, one does not expect to have many reactions hitting the bounds
% 2. LLCs invoke a subset of all loopless constraints and the associated binary
% variables. Maximize everything at once will require invoking almost all binary variables
QuickProblem = addLoopLawConstraints(QuickProblem, model, 1:nRxns, [], [], loopInfo);
end
%Maximise all reactions
QuickProblem.osense = -1;
quickSolultionFound = false;
switch loopMethod
case {'original', 'fastSNP'}
% Set a short time limit or do not this at all for MILP with loop law
% because if the model and the number of reactions in the objective is large,
% it is non-trivial for solvers to find the optimal solution under the loop law.
idTimeLimit = strcmp(solverVarargin.MILP, 'timeLimit');
if any(idTimeLimit)
idTimeLimit(find(idTimeLimit) + 1) = true;
end
sol = solveCobraMILP(QuickProblem, solverVarargin.MILP{~idTimeLimit}, 'timeLimit', 10);
if sol.stat == 1 || (sol.stat == 3 && ~isempty(sol.full))
% accept if there is a feasible solution for the MILP
quickSolultionFound = true;
end
case 'none'
sol = solveCobraLP(QuickProblem, solverVarargin.LP{:});
if sol.stat == 1 || checkSolFeas(QuickProblem, sol) <= cobraParams.LP.feasTol
quickSolultionFound = true;
end
end
% If we reach this point, we can be certain, that there is a solution, i.e.
% if the stat is not 1, we have to check all reactions.
if quickSolultionFound % accept if there is a feasible solution for the MILP
% Obtain fluxes at their boundaries
maxSolved = model.ub(Order) == sol.full(Order);
minSolved = model.lb(Order) == sol.full(Order);
% If preCompMaxSols/preCompMinSols is non-empty, no need to solve LP again
sol.heuristics = 'maxSumFlux';
heuristicSolutions{2} = sol;
[preCompMaxSols(maxSolved), preCompMinSols(minSolved)] = deal(2);
end
%Minimise reactions
QuickProblem.osense = 1;
quickSolultionFound = false;
switch loopMethod
case {'original', 'fastSNP'}
% Set a short time limit or do not this at all for MILP with loop law
% because if the model and the number of reactions in the objective is large,
% it is non-trivial for solvers to find the optimal solution under the loop law.
sol = solveCobraMILP(QuickProblem, solverVarargin.MILP{~idTimeLimit}, 'timeLimit', 10);
if sol.stat == 1 || (sol.stat == 3 && ~isempty(sol.full))
% accept if there is a feasible solution for the MILP
quickSolultionFound = true;
end
case 'none'
sol = solveCobraLP(QuickProblem, solverVarargin.LP{:});
if sol.stat == 1 || checkSolFeas(QuickProblem, sol) <= cobraParams.LP.feasTol
quickSolultionFound = true;
end
end
if quickSolultionFound
%Again obtain fluxes at their boundaries
maxSolved = maxSolved | (model.ub(Order) == sol.full(Order));
minSolved = minSolved | (model.lb(Order) == sol.full(Order));
% If preCompMaxSols/preCompMinSols is non-empty, no need to solve LP again
sol.heuristics = 'minSumFlux';
heuristicSolutions{3} = sol;
[preCompMaxSols(model.ub(Order) == sol.full(Order)), preCompMinSols(model.lb(Order) == sol.full(Order))] = deal(3);
end
end
% generate the loopless problem beforehand instead of during each loop
[MILPproblem, LPproblemLLC] = deal([]);
switch loopMethod
case {'original', 'fastSNP'}
% always run MILP. Can directly replace the LP
LPproblem = addLoopLawConstraints(LPproblem, model, 1:nRxns, [], [], loopInfo);
case {'LLC-NS', 'LLC-EFM'}
% keep both LP and MILP. Also update the constraint indices in loopInfo
[MILPproblem, loopInfo] = addLoopLawConstraints(LPproblem, model, 1:nRxns, [], [], loopInfo);
% store the default RHS without problem-specific LLCs
loopInfo.rhs0 = MILPproblem.b;
LPproblemLLC = LPproblem;
end
if heuristics > 1
% solve one LP to find all blocked irreversible reactions
if any((model.lb(Order) >= 0 | model.ub(Order) <= 0) & ~(minSolved & maxSolved))
[~, sol] = findBlockedIrrRxns(model, [], solverVarargin.LP{:});
% for irreversible reactions with fluxes < feasTol, we can safely say that they are blocked
sol.full(abs(sol.full) < cobraParams.LP.feasTol) = 0;
rxnBlocked = (sol.full(Order) == 0) & (model.lb(Order) >= 0 | model.ub(Order) <= 0);
if allowLoops
rxnHitUB = (sol.full(Order) == model.ub(Order)) & ~rxnBlocked & ~maxSolved;
rxnHitLB = (sol.full(Order) == model.lb(Order)) & ~rxnBlocked & ~minSolved;
else
% the reactions that hit bounds in the LP may not hit bounds under the loopless condition
[rxnHitUB, rxnHitLB] = deal(false(numel(rxnNameList), 1));
end
sol.heuristics = 'hitBounds';
if any(rxnHitUB) || any(rxnHitLB)
% store the heuristic solutions
heuristicSolutions{4} = sol;
[preCompMaxSols(rxnHitUB), preCompMinSols(rxnHitLB)] = deal(4);
end
if any(rxnBlocked)
if minNorm
% if flux distributions are to be returned, get one for one of the blocked reactions. It works for all blocked reactions
allowLoopsI = allowLoops;
% For LLCs, solve LP if the problem constraints do not necessitate the loop law and the target reaction has its forward diretion in loops
if strncmpi(loopMethod, 'LLC', 3)
[allowLoopsI, MILPproblem] = processingLLCs('update', loopInfo, 'max', MILPproblem, zeros(nRxns, 1));
if allowLoopsI
% solving LP is sufficient
LPproblem = LPproblemLLC;
else
% need to solve MILP
LPproblem = MILPproblem;
end
end
[~, V] = calcSolForEntry(model, Order(find(rxnBlocked, 1)), LPproblem, method, allowLoopsI, minNorm, solverVarargin, sol, 1);
sol.fluxMinNorm = V;
end
% store the heuristic solutions
sol.heuristics = 'blockedIrr';
heuristicSolutions{5} = sol;
[preCompMaxSols(rxnBlocked), preCompMinSols(rxnBlocked)] = deal(5);
end
[minSolved, maxSolved] = deal(minSolved | rxnBlocked | rxnHitLB, maxSolved | rxnBlocked | rxnHitUB);
end
end
if useMtFVA
for i = 1:length(rxnNameList)
% retrieve max/min values from heuristic solutions
if minSolved(i)
minFlux(i) = calcSolForEntry([], Order(i) ,LPproblem, method, allowLoops, minNorm, [], heuristicSolutions{preCompMinSols(i)}, 1);
end
if maxSolved(i)
maxFlux(i) = calcSolForEntry([], Order(i) ,LPproblem, method, allowLoops, minNorm, [], heuristicSolutions{preCompMaxSols(i)}, -1);
end
end
if any(~minSolved | ~maxSolved)
[fvalb, fvaub]= mtFVA(LPproblem, [columnVector(Order(~maxSolved)); columnVector(-Order(~minSolved))], solverVarargin.LP{1});
minFlux(~minSolved) = fvalb(Order(~minSolved));
maxFlux(~maxSolved) = fvaub(Order(~maxSolved));
end
return
end
if ~parallelJob % single-thread FVA
if printLevel == 1
showprogress(0,'Single-thread flux variability analysis in progress ...');
elseif printLevel > 1
fprintf('%4s\t%4s\t%10s\t%9s\t%9s\n','No','Perc','Name','Min','Max');
end
% Do this to keep the progress printing in place
for i = 1:length(rxnNameList)
rxnID = findRxnIDs(model, rxnNameList(i));
objVector = sparse(rxnID, 1, 1, nRxns, 1);
%Calc minimums
allowLoopsI = allowLoops;
% For LLCs, solve LP if the problem constraints do not necessitate the
% loop law and the target reaction has its reverse diretion in loops
if (~minSolved(i) || minNorm) && strncmpi(loopMethod, 'LLC', 3)
[allowLoopsI, MILPproblem] = processingLLCs('update', loopInfo, 'min', MILPproblem, objVector);
if allowLoopsI
% solving LP is sufficient
LPproblem = LPproblemLLC;
else
% need to solve MILP
LPproblem = MILPproblem;
end
end
[minFlux(i), V] = calcSolForEntry(model, rxnID ,LPproblem, method, allowLoopsI, minNorm, solverVarargin, heuristicSolutions{preCompMinSols(i)}, 1);
% store the flux distribution
if minNorm
Vmin(:, i) = V;
end
%Calc maximums
allowLoopsI = allowLoops;
% For LLCs, solve LP if the problem constraints do not necessitate the
% loop law and the target reaction has its forward diretion in loops
if (~maxSolved(i) || minNorm) && strncmpi(loopMethod, 'LLC', 3)
[allowLoopsI, MILPproblem] = processingLLCs('update', loopInfo, 'max', MILPproblem, objVector);
if allowLoopsI
% solving LP is sufficient
LPproblem = LPproblemLLC;
else
% need to solve MILP
LPproblem = MILPproblem;
end
end
[maxFlux(i), V] = calcSolForEntry(model, rxnID ,LPproblem, method, allowLoopsI, minNorm, solverVarargin, heuristicSolutions{preCompMaxSols(i)}, -1);
% store the flux distribution
if minNorm
Vmax(:, i) = V;
end
if printLevel == 1
showprogress(i/length(rxnNameList));
end
if printLevel > 1
fprintf('%4d\t%4.0f\t%10s\t%9.3f\t%9.3f\n',i,100*i/length(rxnNameList),rxnNameList{i},minFlux(i),maxFlux(i));
end
end
else % parallel job. pretty much does the same thing.
environment = getEnvironment();
if printLevel == 1
fprintf('Parallel flux variability analysis in progress ...\n');
elseif printLevel > 1
fprintf('%4s\t%10s\t%9s\t%9s\n','No','Name','Min','Max');
end
parfor i = 1:length(rxnNameList)
restoreEnvironment(environment,0);
parLPproblem = LPproblem;
rxnID = findRxnIDs(model, rxnNameList(i));
objVector = sparse(rxnID, 1, 1, nRxns, 1);
%Calc minimums
allowLoopsI = allowLoops;
% For LLCs, solve LP if the problem constraints do not necessitate the
% loop law and the target reaction has its reverse diretion in loops
if (~minSolved(i) || minNorm) && strncmpi(loopMethod, 'LLC', 3)
[allowLoopsI, parMILPproblem] = processingLLCs('update', loopInfo, 'min', MILPproblem, objVector);
if allowLoopsI
% solving LP is sufficient
parLPproblem = LPproblemLLC;
else
% need to solve MILP
parLPproblem = parMILPproblem;
end
end
[minFlux(i), V] = calcSolForEntry(model, rxnID ,parLPproblem, method, allowLoopsI, minNorm, solverVarargin, heuristicSolutions{preCompMinSols(i)}, 1);
% store the flux distribution
if minNorm
Vmin(:, i) = V;
end
%Calc maximums
allowLoopsI = allowLoops;
% For LLCs, solve LP if the problem constraints do not necessitate the
% loop law and the target reaction has its forward diretion in loops
if (~maxSolved(i) || minNorm) && strncmpi(loopMethod, 'LLC', 3)
[allowLoopsI, parMILPproblem] = processingLLCs('update', loopInfo, 'max', MILPproblem, objVector);
if allowLoopsI
% solving LP is sufficient
parLPproblem = LPproblemLLC;
else
% need to solve MILP
parLPproblem = parMILPproblem;
end
end
[maxFlux(i), V] = calcSolForEntry(model, rxnID ,parLPproblem, method, allowLoopsI, minNorm, solverVarargin, heuristicSolutions{preCompMaxSols(i)}, -1);
% store the flux distribution
if minNorm
Vmax(:, i) = V;
end
if printLevel > 1
fprintf('%4d\t%10s\t%9.3f\t%9.3f\n', i, rxnNameList{i}, minFlux(i), maxFlux(i));
end
end
end
maxFlux = columnVector(maxFlux);
minFlux = columnVector(minFlux);
end
function [Flux, V] = calcSolForEntry(model, rxnID, LPproblem, method, allowLoops, minNorm, solverVarargin, sol, osense)
%Set the correct objective
LPproblem.osense = sign(osense);
LPproblem.f(:) = 0;
LPproblem.f(rxnID) = 1;
if isempty(sol)
if allowLoops
% solve LP
LPsolution = solveCobraLP(LPproblem, solverVarargin.LP{:});
else
% solve MILP
LPsolution = solveCobraMILP(LPproblem, solverVarargin.MILP{:});
end
% take the maximum flux from the flux vector, not from the obj -Ronan
% A solution is possible, so the only problem should be if its
% unbounded and if it is unbounded, the max flux is infinity.
if LPsolution.stat == 2
Flux = -LPproblem.osense * inf;
elseif LPsolution.stat == 1
Flux = getObjectiveFlux(LPsolution, LPproblem);
else
error(sprintf(['A Solution could not be found!\nThis should not be possible but can happen',...
'if the used solver cannot properly handle unboundedness, or if there are numerical issues.\n',...
'Please try to use a different solver.\n']))
end
else
% use pre-computed solutions from heuristics
Flux = sol.full(rxnID);
if minNorm
LPsolution = sol;
end
end
V = [];
if minNorm
if ~isempty(sol) && strcmp(sol.heuristics, 'blockedIrr')
% use the solution calculated during heuristics
V = sol.fluxMinNorm;
elseif allowLoops
V = getMinNorm(LPproblem, LPsolution, numel(model.rxns), Flux, model, method, solverVarargin);
else
V = getMinNormWoLoops(LPproblem, LPsolution, numel(model.rxns), Flux, method, solverVarargin);
end
end
end
function V = getMinNorm(LPproblem, LPsolution, nRxns, cFlux, model, method, solverVarargin)
% get the Flux distribution for the specified min norm.
% update LPproblem to fix objective function value for 1-norm and
% 0-norm to work
LPproblem.lb(LPproblem.f ~= 0) = cFlux - 1e-12;
LPproblem.ub(LPproblem.f ~= 0) = cFlux + 1e-12;
switch method
case '2-norm'
LPproblem.f(:)=0;
%Minimise Euclidean norm using quadratic programming
LPproblem.F = [speye(nRxns, nRxns), sparse(nRxns, size(LPproblem.A, 2) - nRxns);...
sparse(size(LPproblem.A, 2) - nRxns, size(LPproblem.A, 2))];
LPproblem.osense = 1;
%quadratic optimization
solution = solveCobraQP(LPproblem, solverVarargin.QP{:});
V = solution.full(1:nRxns, 1);
case '1-norm'
V = sparseFBA(LPproblem, 'min', 0, 0, 'l1');
case '0-norm'
V = sparseFBA(LPproblem, 'min', 0, 0);
case 'FBA'
V= LPsolution.full(1:nRxns);
case 'minOrigSol'
% we take the original model, and constrain the objective reaction
% accordingly.
LPproblemMOMA = model;
LPproblemMOMA.lb(LPproblem.f(1:nRxns)~=0) = cFlux - 1e-11;
LPproblemMOMA.ub(LPproblem.f(1:nRxns)~=0) = cFlux + 1e-11;
momaSolution = linearMOMA(model,LPproblemMOMA);
V = momaSolution.x;
end
end
function V = getMinNormWoLoops(MILPproblem, MILPsolution, nRxns, cFlux, method, solverVarargin)
% It will be great if sparseFBA can somehow support MILP problems
[m, n] = size(MILPproblem.A);
% too small gap between lb and ub may cause numerical difficulty for solvers
MILPproblem.lb(MILPproblem.f ~= 0) = floor(cFlux * 1e9) / 1e9;
MILPproblem.ub(MILPproblem.f ~= 0) = ceil(cFlux * 1e9) / 1e9;
switch method
case '2-norm'
MILPproblem.f(:)=0;
%Minimise Euclidean norm using quadratic programming
MILPproblem.F = [speye(nRxns,nRxns), sparse(nRxns, n - nRxns); ...
sparse(n - nRxns, n)];
MILPproblem.osense = 1;
% supplying a known initial solution has a much higher chance for the solver to return a solution
MILPproblem.x0 = MILPsolution.full;
%quadratic optimization
solution = solveCobraMIQP(MILPproblem, solverVarargin.MIQP{:});
V = solution.full(1:nRxns);
case '1-norm'
MILPproblem.A = [MILPproblem.A, sparse(m, nRxns); ... original problem
sparse(1:nRxns, 1:nRxns, 1, nRxns, n), -speye(nRxns); ... v - |v| <= 0
sparse(1:nRxns, 1:nRxns, -1, nRxns, n), -speye(nRxns)]; % -v - |v| <= 0
MILPproblem.b = [MILPproblem.b; zeros(nRxns * 2, 1)];
MILPproblem.csense = [MILPproblem.csense(:); repmat('L', nRxns * 2, 1)];
MILPproblem.f = [zeros(n, 1); ones(nRxns, 1)];
MILPproblem.osense = 1;
MILPproblem.vartype = [MILPproblem.vartype(:); repmat('C', nRxns, 1)];
MILPproblem.lb = [MILPproblem.lb; zeros(nRxns, 1)];
MILPproblem.ub = [MILPproblem.ub; max(abs([MILPproblem.lb(1:nRxns), MILPproblem.ub(1:nRxns)]), [], 2)];
% supplying a known initial solution has a much higher chance for the solver to return a solution
MILPproblem.x0 = [MILPsolution.full; abs(MILPsolution.full(1:nRxns))];
solution = solveCobraMILP(MILPproblem, solverVarargin.MILP{:});
V = solution.full(1:nRxns);
case '0-norm'
% use binary switch
MILPproblem.A = [MILPproblem.A, sparse(m, nRxns); ... original problem
sparse(1:nRxns, 1:nRxns, 1, nRxns, n), sparse(1:nRxns, 1:nRxns, -MILPproblem.ub(1:nRxns), nRxns, nRxns); ... v - ub*a <= 0
sparse(1:nRxns, 1:nRxns, -1, nRxns, n), sparse(1:nRxns, 1:nRxns, MILPproblem.lb(1:nRxns), nRxns, nRxns)]; % -v + lb*a <= 0
MILPproblem.b = [MILPproblem.b; zeros(nRxns * 2, 1)];
MILPproblem.f = [zeros(n, 1); ones(nRxns, 1)];
MILPproblem.csense = [MILPproblem.csense(:); repmat('L', nRxns * 2, 1)];
MILPproblem.osense = 1;
MILPproblem.vartype = [MILPproblem.vartype(:); repmat('B', nRxns, 1)];
MILPproblem.lb = [MILPproblem.lb; zeros(nRxns, 1)];
MILPproblem.ub = [MILPproblem.ub; ones(nRxns, 1)];
intTol = getCobraSolverParams('MILP', 'intTol');
% supplying a known initial solution has a much higher chance for the solver to return a solution
V = MILPsolution.full(1:nRxns);
V(V > 0 & MILPproblem.ub(1:nRxns) > 0 & V ./ MILPproblem.ub(1:nRxns) <= intTol) = 0;
V(V < 0 & MILPproblem.lb(1:nRxns) < 0 & V ./ MILPproblem.lb(1:nRxns) <= intTol) = 0;
MILPproblem.x0 = [MILPsolution.full; V ~= 0];
solution = solveCobraMILP(MILPproblem, solverVarargin.MILP{:});
V = solution.full(1:nRxns);
case 'FBA'
V= MILPsolution.full(1:nRxns);
case 'minOrigSol'
warning('method ''minOrigSol'' not supported with ''allowLoops'' turned on. Return the ''FBA'' solution');
V= MILPsolution.full(1:nRxns);
end
end
function flux = getObjectiveFlux(LPsolution,LPproblem)
% Determine the current flux based on an LPsolution, the original LPproblem
% The LPproblem is used to retrieve the current objective position.
% min indicates, whether the minimum or maximum is requested, the
% upper/lower bounds are used, if the value is exceeding them
Index = LPproblem.f~=0;
if LPsolution.full(Index)<LPproblem.lb(Index) %takes out tolerance issues
flux = LPproblem.lb(Index);
elseif LPsolution.full(Index)>LPproblem.ub(Index)
flux = LPproblem.ub(Index);
else
flux = LPsolution.full(Index);
end
end