ManualINITTests

These tests are using the Human-GEM model.

0001 %These tests are using the Human-GEM model.
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%
0004 % TC0901 - Test that all reactions in the minimal model to be
0005 %          used by the MILP can carry flux > 0.1
0006 %%%%%%%%%%%%%%%%%%%%%%%%%%
0007 
0008 %In the MILP, the flux of reactions with negative rxn scores cannot exceed 100
0009 %At the same time, the flux of positive reactions need to be >= 0.1 to be able
0010 %to turn them on. Likewise, we force the flux on essential reactions to be >= 0.1.
0011 %Therefore, it is important that all reactions can get a flux > 0.1, while all
0012 %reactions are limited to a flux <= 100. We can test this in a similar way to the
0013 %haveFlux function in RAVEN, so that code is copied here and modified.
0014 
0015 %This code is partly copied from the function haveFlux
0016 %The prepData takes a couple of hours to generate. If you have it somewhere,
0017 %better to load it
0018 %cd C:\Work\MatlabCode\components\human-GEM\Human-GEMftINIT\Human-GEM %replace this with the root of the Human-GEM repo.
0019 %ihuman = importYaml('model/Human-GEM.yml');
0020 %save('model/Human-GEM.mat', 'ihuman');
0021 %load('model/Human-GEM.mat')
0022 %prepData = prepHumanModelForftINIT(ihuman, false);
0023 prepData2 = load('model/PrepData2.mat').prepData;
0024 
0025 upperLimit = 100;
0026 lowerLimit = 0.1;
0027 
0028 model = prepData2.minModel;
0029 model = prepData.minModel;
0030 
0031 
0032 
0033 
0034 model.lb(model.lb==-inf)=-upperLimit;
0035 model.ub(model.ub==inf)=upperLimit;
0036 model.lb(model.lb==-1000)=-upperLimit;
0037 model.ub(model.ub==1000)=upperLimit;
0038 
0039 %check
0040 unique(model.lb) %ok, only 0 and -100
0041 unique(model.ub) %ok, only 100
0042 
0043 %First make a loop where we optimize for all
0044 onFwd = false(numel(model.rxns,1));
0045 onRev = false(numel(model.rxns,1));
0046 
0047 iter = 1;
0048 while true
0049     disp([num2str(iter) ': ' num2str(sum(onFwd))])
0050     model.c=ones(numel(model.c),1);
0051     model.c(onFwd) = 0; %don't include the ones already verified
0052     
0053     sol=solveLP(model);
0054     if isempty(sol.x)
0055         disp('Failed');
0056         break;
0057     end
0058     
0059     currRes = sol.x > lowerLimit;
0060     if sum(currRes & ~onFwd) == 0
0061         break;%we didn't turn any more rxns on
0062     end
0063     
0064     onFwd = onFwd | currRes;
0065     iter = iter + 1;
0066 end
0067 
0068 sum(~onFwd)%zero, perfect
0069 
0070 %Now the other direction (reversible):
0071 onRev = false(numel(model.rxns,1));
0072 
0073 iter = 1;
0074 while true
0075     disp([num2str(iter) ': ' num2str(sum(onRev))])
0076     model.c=-ones(numel(model.c),1);
0077     model.c(~model.rev) = 0;
0078     model.c(onRev) = 0; %don't include the ones already verified
0079     
0080     sol=solveLP(model);
0081     if isempty(sol.x)
0082         disp('Failed');
0083         break;
0084     end
0085     
0086     currRes = sol.x < -lowerLimit;
0087     if sum(currRes & ~onRev) == 0
0088         break;%we didn't turn any more rxns on
0089     end
0090     
0091     onRev = onRev | currRes;
0092     iter = iter + 1;
0093 end
0094 
0095 sum(~onRev(model.rev == 1))%
0096 
0097 onRev = false(numel(model.rxns,1));
0098 
0099 indLeft = find(model.rev == 1 & ~onRev);
0100 
0101 constructEquations(model, model.rxns(indLeft))%no scaling here...
0102 
0103 %go trough the last one by one
0104 for i = 1:numel(indLeft)
0105     model.c=zeros(numel(model.c),1);
0106     model.c(indLeft(i)) = -1;
0107     
0108     sol=solveLP(model);
0109     if isempty(sol.x)
0110         disp('Failed');
0111         break;
0112     end
0113     
0114     currRes = sol.x(indLeft(i)) < -lowerLimit;
0115     disp([num2str(indLeft(i)) ': ' num2str(currRes)])
0116     onRev(indLeft(i)) = currRes;
0117     onRev = onRev | currRes;
0118     iter = iter + 1;
0119 end
0120 
0121 sum(~onRev(model.rev == 1))%0, perfect. This means that we can use the limits 0.1 and 100 in the MILP!
0122 
0123