Extract a thermodynamically consistent subnetwork from a given model

  1. Identify the largest subset of a model that admits a thermodynamically consistent flux
  2. Specify a random subset of active/inactive reactions and present/absent metabolites
  3. Compute the smallest thermodynamically consistent subnetwork containing a list of present metabolites and active reactions, and not containing a list of absent metabolites and inactive reactions
[solverOK,solverInstalled]=changeCobraSolver('ibm_cplex','all');
> changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path. > ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2021a on your operating system. > changeCobraSolver: Solver for LP problems has been set to ibm_cplex. > changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path. > ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2021a on your operating system. > changeCobraSolver: Solver for MILP problems has been set to ibm_cplex. > changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path. > ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2021a on your operating system. > changeCobraSolver: Solver for QP problems has been set to ibm_cplex. > changeCobraSolver: IBM ILOG CPLEX interface added to MATLAB path. > ibm_cplex (version 1210) is compatible and fully tested with MATLAB R2021a on your operating system. > changeCobraSolver: Solver for MIQP problems has been set to ibm_cplex. > changeCobraSolver: Solver ibm_cplex not supported for problems of type EP. No solver set for this problemtype > changeCobraSolver: Solver ibm_cplex not supported for problems of type NLP. No solver set for this problemtype
%[solverOK,solverInstalled]=changeCobraSolver('gurobi','all');
%[solverOK,solverInstalled]=changeCobraSolver('ibm_cplex','QP');

Load model

modelToLoad='circularToy';
modelToLoad='ecoli_core';
modelToLoad='modelRecon3MitoOpen';
modelToLoad='Recon3DModel';
%modelToLoad='iDopa';
Load a model
driver_thermoModelLoad
Model loaded: Recon3DModel lower bounds greater than zero Internal stochiometric nullspace computed in 4.5238 seconds.

Stoichiometric consistency

if ~isfield(model,'SConsistentRxnBool') || ~isfield(model,'SConsistentMetBool')
massBalanceCheck=0;
%massBalanceCheck=1;
printLevel=2;
[SConsistentMetBool, SConsistentRxnBool, SInConsistentMetBool, SInConsistentRxnBool, unknownSConsistencyMetBool, unknownSConsistencyRxnBool, model,stoichConsistModel]...
= findStoichConsistentSubset(model, massBalanceCheck, printLevel);
else
%Extract stoich consistent submodel
if any(~model.SConsistentMetBool)
rxnRemoveMethod='inclusive';%maintains stoichiometric consistency
[stoichConsistModel, rxnRemoveList] = removeMetabolites(model, model.mets(~model.SConsistentMetBool),rxnRemoveMethod);
SConsistentRxnBool2=~ismember(model.rxns,rxnRemoveList);
if ~all(model.SConsistentRxnBool==SConsistentRxnBool2)
error('inconsistent reaction removal')
end
try
stoichConsistModel = removeUnusedGenes(stoichConsistModel);
catch ME
disp(ME.message)
end
else
stoichConsistModel = model;
end
end
 
[nMet,nRxn]=size(stoichConsistModel.S)
nMet = 5835
nRxn = 10600

Flux consistency

fluxConsistentParam.method='fastcc';%can handle additional constraints
fluxConsistentParam.printLevel=1;
[~,~,~,~,stoichConsistModel]= findFluxConsistentSubset(stoichConsistModel,fluxConsistentParam);
Extract flux consistent submodel
if any(~stoichConsistModel.fluxConsistentRxnBool)
rxnRemoveList = stoichConsistModel.rxns(~stoichConsistModel.fluxConsistentRxnBool);
stoichFluxConsistModel = removeRxns(stoichConsistModel, rxnRemoveList,'metRemoveMethod','exclusive','ctrsRemoveMethod','inclusive');
try
stoichFluxConsistModel = removeUnusedGenes(stoichFluxConsistModel);
catch ME
disp(ME.message)
end
else
stoichFluxConsistModel = stoichConsistModel;
end
[nMet,nRxn]=size(stoichFluxConsistModel.S)
nMet = 5835
nRxn = 10600

Thermodynamic consistency

%save('debug_prior_to_findThermoConsistentFluxSubset.mat')
%return
param.printLevel = 1;
param.relaxBounds=0;
param.acceptRepairedFlux=1;
[thermoFluxConsistentMetBool,thermoFluxConsistentRxnBool,stoichFluxConsistModel,stoichFluxThermoConsistModel] = findThermoConsistentFluxSubset(stoichFluxConsistModel,param);
--- findThermoFluxConsistentSubset START ---- printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 epsilon: 1.0000e-06 formulation: 'pqzw' iterationMethod: 'random' nMax: 20 warmStartMethod: 'random' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 0 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 0.1 = beta, the global weight on one-norm of internal reaction rate. -5 = min(g0), the local weight on zero-norm of internal reaction rate. -0 = max(g0), the local weight on zero-norm of internal reaction rate. 0 = min(h0), the local weight on zero-norm of metabolite production rate. 0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 0 min cardinality variables: NaN mean(c(p)) NaN min(c(p)) NaN max(c(p)) 1 lambda0 NaN min(k) NaN max(k) 1 lambda1 NaN min(o(p)) NaN max(o(p)) 5303 max cardinality variables: -0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 5 min(d) 5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 22879 cardinality free variables: 0.077 mean(c(r)) -0 min(c(r)) 0.1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 7.6657e+05 -8.8e+06 -2.2e+04 1.6e+03 0 0 0 -23650 -2.3e+04 0 NaN 0 1.236540e+00 2 0.75 184.22 -8.5e+02 -2.3e+04 1.1e+03 0 0 0 -23840 -2.4e+04 0 NaN 0 1.592884e+00 3 1.12 98.541 -4.8e+02 -2.3e+04 7.3e+02 0 0 0 -23900 -2.4e+04 0 NaN 0 1.325164e+00 4 1.69 61.895 -2.6e+02 -2.3e+04 4.9e+02 0 0 0 -23910 -2.4e+04 0 NaN 0 1.459285e+00 5 2.53 40.001 -1.7e+02 -2.4e+04 3.3e+02 0 0 0 -23915 -2.4e+04 0 NaN 0 1.535346e+00 6 3.80 26.092 -1.3e+02 -2.4e+04 2.2e+02 0 0 0 -23935 -2.4e+04 0 NaN 0 1.639884e+00 7 5.70 18.183 -73 -2.4e+04 1.5e+02 0 0 0 -23940 -2.4e+04 0 NaN 0 1.398942e+00 8 8.54 11.708 -68 -2.4e+04 97 0 0 0 -23955 -2.4e+04 0 NaN 0 1.285520e+00 9 12.81 7.942 -34 -2.4e+04 65 0 0 0 -23965 -2.4e+04 0 NaN 0 1.246812e+00 10 19.22 5.2426 -36 -2.4e+04 44 0 0 0 -23990 -2.4e+04 0 NaN 0 1.474974e+00 11 28.83 3.5411 -48 -2.4e+04 29 0 0 0 -24010 -2.4e+04 0 NaN 0 1.207763e+00 12 43.25 2.3919 -17 -2.4e+04 20 0 0 0 -24020 -2.4e+04 0 NaN 0 1.119197e+00 13 64.87 1.6485 -13 -2.4e+04 14 0 0 0 -24020 -2.4e+04 0 NaN 0 1.071447e+00 14 97.31 1.0952 -4.2 -2.4e+04 9.3 0 0 0 -24020 -2.4e+04 0 NaN 0 1.335563e+00 15 145.96 0.76091 -2.9 -2.4e+04 6.5 0 0 0 -24040 -2.4e+04 0 NaN 0 1.300141e+00 16 218.95 0.51853 -42 -2.4e+04 4.7 0 0 0 -24060 -2.4e+04 0 NaN 0 1.308127e+00 17 328.42 0.53726 -14 -2.4e+04 3.5 0 0 0 -24090 -2.4e+04 0 NaN 0 1.286361e+00 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. cycleFreeFlux: No solution found. Debugging relaxation etc... full: [] obj: [] rcost: [] dual: [] slack: [] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 0 origStat: 3 origStatText: 'Model has been proved infeasible' time: 0.0300 basis: [] itn obj obj_old err(obj) err(x) card(v) card(r) card(p) card(q) 0 14.072 66705 66691 165.16 3913 0 543 541 1 1391 14.072 1376.9 48.563 2734 0 692 717 2 944.42 1391 446.54 62.844 3799 0 652 666 3 1526.3 944.42 581.91 64.636 3026 0 712 702 4 828.59 1526.3 697.73 63.09 3157 0 574 605 5 992.49 828.59 163.9 55.391 3334 0 716 682 6 1074.4 992.49 81.866 59.058 3674 0 687 657 7 1300.5 1074.4 226.18 63.35 3238 0 700 708 8 1040.9 1300.5 259.63 62.911 3773 0 685 639 9 1327.8 1040.9 286.85 63.278 3013 0 703 667 10 814.59 1327.8 513.18 60.13 3064 0 614 603 11 1027.5 814.59 212.94 55.715 2717 0 628 614 12 677.72 1027.5 349.8 53.573 3435 0 602 618 13 1061.2 677.72 383.44 54.092 3505 0 734 693 14 1272.7 1061.2 211.54 62.794 3509 0 693 670 15 1197.8 1272.7 74.901 64.628 3152 0 669 689 16 971.72 1197.8 226.09 60.517 2983 0 585 622 17 903.71 971.72 68.006 56.208 2923 0 612 658 18 834.71 903.71 69 54.078 3806 0 671 642 19 1348.3 834.71 513.59 60.702 3075 0 720 659 itn obj obj_old err(obj) err(x) card(v) card(r) card(p) card(q) Relaxed model is feasible. Statistics: 599 lower bound relaxation(s) 576 upper bound relaxation(s) 0 steady state relaxation(s) ... done.
ans = 9.1346e-07
ans = 0
ans = 0
solution = struct with fields:
stat: 1 v: [19391×1 double] r: [23417×1 double] p: [19391×1 double] q: [19391×1 double]
solutionRelaxed1 = struct with fields:
full: [19391×1 double] obj: 26.3611 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 1.2415 basis: []
solutionRelaxed2 = struct with fields:
full: [19391×1 double] obj: 0 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 0.6092 basis: []
solutionRelaxed3 = struct with fields:
full: [19391×1 double] obj: 0 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 0.3797 basis: []
cycleFreeFlux: No solution found, try using a different solver. cycleFreeFlux: No solution found. Debugging relaxation etc... full: [] obj: [] rcost: [] dual: [] slack: [] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 0 origStat: 3 origStatText: 'Model has been proved infeasible' time: 0.0305 basis: [] itn obj obj_old err(obj) err(x) card(v) card(r) card(p) card(q) 0 14.072 66355 66341 164.43 3913 0 543 541 1 1391 14.072 1376.9 48.563 2734 0 692 717 2 944.42 1391 446.54 62.844 3799 0 652 666 3 1526.3 944.42 581.91 64.636 3026 0 712 702 4 828.59 1526.3 697.73 63.09 3157 0 574 605 5 992.49 828.59 163.9 55.391 3334 0 716 682 6 1074.4 992.49 81.866 59.058 3674 0 687 657 7 1300.5 1074.4 226.18 63.35 3238 0 700 708 8 1040.9 1300.5 259.63 62.911 3773 0 685 639 9 1327.8 1040.9 286.85 63.278 3013 0 703 667 10 814.59 1327.8 513.18 60.13 3064 0 614 603 11 1027.5 814.59 212.94 55.715 2717 0 628 614 12 677.72 1027.5 349.8 53.573 3435 0 602 618 13 1061.2 677.72 383.44 54.092 3505 0 734 693 14 1272.7 1061.2 211.54 62.794 3509 0 693 670 15 1197.8 1272.7 74.901 64.628 3152 0 669 689 16 971.72 1197.8 226.09 60.517 2983 0 585 622 17 903.71 971.72 68.006 56.208 2923 0 612 658 18 834.71 903.71 69 54.078 3806 0 671 642 19 1348.3 834.71 513.59 60.702 3075 0 720 659 itn obj obj_old err(obj) err(x) card(v) card(r) card(p) card(q) Relaxed model is feasible. Statistics: 599 lower bound relaxation(s) 576 upper bound relaxation(s) 0 steady state relaxation(s) ... done.
ans = 9.1346e-07
ans = 0
ans = 0
solution = struct with fields:
stat: 1 v: [19391×1 double] r: [23417×1 double] p: [19391×1 double] q: [19391×1 double]
solutionRelaxed1 = struct with fields:
full: [19391×1 double] obj: 26.3611 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 1.3163 basis: []
solutionRelaxed2 = struct with fields:
full: [19391×1 double] obj: 0 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 0.5761 basis: []
solutionRelaxed3 = struct with fields:
full: [19391×1 double] obj: 0 rcost: [19391×1 double] dual: [23417×1 double] slack: [23417×1 double] solver: 'ibm_cplex' algorithm: 'Automatic' stat: 1 origStat: 1 origStatText: 'Optimal solution found' time: 0.3650 basis: []
Warning: cycleFreeFlux did not solve, trying v2QNTy
76.81% thermodynamically feasible internal fluxes (checked by v2QNty method). iter card(y) nz %feas int.nz. tot %feas int.nz. tot %feas ext.nz. iteration formulation 1 5303 7882 0.68 0.50 0.80 random pqzw 2 2376 3075 1.00 0.69 0.93 random pqzw 3 1427 1708 1.00 0.77 0.98 random pqzw 4 1022 1296 1.00 0.83 0.99 random pqzw 5 769 794 1.00 0.86 1.00 random pqzw 6 606 587 1.00 0.88 1.00 random pqzw 7 506 432 1.00 0.90 1.00 random pqzw 8 449 377 1.00 0.91 1.00 random pqzw 9 389 324 1.00 0.92 1.00 random pqzw 10 348 301 1.00 0.93 1.00 random pqzw 11 298 332 1.00 0.94 1.00 random pqzw 12 254 384 1.00 0.95 1.00 random pqzw 13 216 209 1.00 0.96 1.00 random pqzw 14 191 200 1.00 0.97 1.00 random pqzw 15 153 201 1.00 0.97 1.00 random pqzw 16 151 162 1.00 0.97 1.00 random pqzw 17 106 139 1.00 0.98 1.00 random pqzw 18 102 161 1.00 0.98 1.00 random pqzw 19 78 127 1.00 0.98 1.00 random pqzw 20 91 128 1.00 0.98 1.00 random pqzw iter card(y) nz %feas int.nz. tot %feas int.nz. tot %feas ext.nz. iteration formulation findThermoConsistentFluxSubset terminating early: n = nMax = 20 --- findThermoFluxConsistentSubset END ----
Size of the largest flux, stoich and thermo consistent submodel
[nMet,nRxn]=size(stoichFluxThermoConsistModel.S)
save(['~/work/sbgCloud/programModelling/projects/thermoModel/results/thermoKernel/' modelToLoad '_stoichFluxThermoConsistModel.mat'],'stoichFluxThermoConsistModel')
%modelToLoad='Recon3DModel';
load(['~/work/sbgCloud/programModelling/projects/thermoModel/results/thermoKernel/' modelToLoad '_stoichFluxThermoConsistModel.mat'],'stoichFluxThermoConsistModel')

Nullspace

Nullspace is necessary for backup check of thermodynamic consistency using thermoFlux2QNty
[stoichFluxThermoConsistModel,rankK,nnzK,timeTaken] = internalNullspace(stoichFluxThermoConsistModel);
rankK
rankK = 5485

Data to define a thermodynamically consistent subnetwork

Setup random data to select a random subset
param.n=200;
[rankMetConnectivity,rankMetInd,rankConnectivity] = rankMetabolicConnectivity(stoichFluxThermoConsistModel,param);
[nMet,nRxn]=size(stoichFluxThermoConsistModel.S);
rxnWeights=rand(nRxn,1)-0.5;
rxnWeights(~stoichFluxThermoConsistModel.SConsistentRxnBool)=0;
 
coreRxnBool=rxnWeights<-0.45;
removeRxnBool=rxnWeights>0.45;
if 0
rxnWeights(rxnWeights>0.4)=1;
rxnWeights(rxnWeights<-0.4)=-1;
end
rxnWeights(rxnWeights>=-0.4 & rxnWeights<=0.4)=0;
hist(rxnWeights,40)
metWeights=rand(nMet,1)-0.5;
metWeights(rankMetInd(1:200))=0;
coreMetBool=metWeights<-0.45;
removeMetBool=metWeights>0.45;
if 0
metWeights(metWeights>0.4)=1;
metWeights(metWeights<-0.4)=-1;
end
metWeights(metWeights>=-0.4 & metWeights<=0.4)=0;
hist(metWeights,40)

nlt=length(coreRxnBool);
activeInactiveRxn=zeros(nlt,1);
activeInactiveRxn(coreRxnBool)=1;
activeInactiveRxn(removeRxnBool)=-1;
hist(activeInactiveRxn)
mlt=length(coreMetBool);
presentAbsentMet=zeros(mlt,1);
presentAbsentMet(coreMetBool)=1;
presentAbsentMet(removeMetBool)=-1;
if 0
activeInactiveRxn(:)=0;
presentAbsentMet(:)=0;
end
param.normalizeZeroNormWeights=0;
 
hist(presentAbsentMet)

Compute the smallest thermodynamically consistent subnetwork given a list of present/absent metabolites and active/inactive reactions

[tissueModel, thermoModelMetBool, thermoModelRxnBool] = thermoKernel(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights,param);
--- thermoKernel START ---- thermoKernel parameters: printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 n: 200 normalizeZeroNormWeights: 0 formulation: 'pqzwrs' epsilon: 1.0000e-06 removeOrphanGenes: 1 nbMaxIteration: 30 nMax: 20 iterationMethod: 'greedyRandom' warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0.5 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0.5 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 1393 min cardinality variables: 0 mean(c(p)) -0 min(c(p)) -0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 1912 max cardinality variables: 0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 35242 cardinality free variables: 0.48 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6992e+05 -8.5e+07 7.9 11 1.64389 0.11 0 -14.9047 -3 0 NaN 0 6.493390e-01 2 0.75 4.5121 -1.7 6.2 12 2.0597 0.16 0 -15.8473 -5.7 0 NaN 0 5.738400e-01 3 1.12 2.8094 -1.8 4.4 12 1.65769 0.25 0 -17.5105 -7.8 0 NaN 0 7.604760e-01 4 1.69 3.0116 -2 2.4 12 1.23553 0.43 0 -18.8108 -10 0 NaN 0 4.662190e-01 5 2.53 3.0252 -2.3 0.041 13 1.23553 0.45 0 -19.6814 -14 0 NaN 0 7.602410e-01 6 3.80 2.7007 -3.3 -3.3 13 1.23553 0.47 0 -24.7884 -17 0 NaN 0 8.377580e-01 7 5.70 2.0612 -4.9 -8.2 14 1.23553 0.87 0 -25.6958 -23 0 NaN 0 9.495780e-01 8 8.54 1.9653 -2.8 -11 13 1.23553 0.9 0 -27.0323 -25 0 NaN 0 6.459370e-01 9 12.81 1.3137 -3.8 -15 12 1.23553 0.95 0 -30.4979 -28 0 NaN 0 6.463110e-01 10 19.22 1.2354 -4 -19 11 1.23553 1 0 -32.3232 -31 0 NaN 0 7.743710e-01 11 28.83 0.809 -2.2 -21 10 1.23553 1.1 0 -32.3232 -32 0 NaN 0 8.423850e-01 12 43.25 0.71837 -0.57 -21 9.6 1.23553 1.2 0 -32.3232 -32 0 NaN 0 6.718740e-01 13 64.87 0.46434 -0.55 -22 9.2 1.23553 1.2 0 -33.2682 -32 0 NaN 0 6.453390e-01 14 97.31 0.33824 -1.1 -23 8.9 1.23553 1.2 0 -33.2682 -33 0 NaN 0 5.601560e-01 15 145.96 0.21621 -0.14 -23 8.8 1.23553 1.2 0 -33.2682 -33 0 NaN 0 7.039990e-01 16 218.95 0.2235 -0.088 -23 8.7 1.23553 1.2 0 -33.2682 -33 0 NaN 0 9.558880e-01 17 328.42 0.10113 -0.058 -23 8.6 1.23553 1.2 0 -33.2682 -33 0 NaN 0 8.566870e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0.5 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0.5 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 1393 min cardinality variables: 0 mean(c(p)) -0 min(c(p)) -0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 1912 max cardinality variables: 0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 35242 cardinality free variables: 0.48 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6862e+05 -8.4e+07 8.1 11 3.74292 0.11 0 -13.2993 -2.7 0 NaN 0 5.372190e-01 2 0.75 2.3033 -1 7.1 10 2.06605 0.16 0 -12.3698 -3.1 0 NaN 0 4.644000e-01 3 1.12 2.5538 -1.1 6 10 2.06605 0.26 0 -13.2447 -4.4 0 NaN 0 5.029130e-01 4 1.69 2.8214 -1.1 4.8 11 1.23553 0.43 0 -14.5616 -6.1 0 NaN 0 4.557920e-01 5 2.53 3.0167 -2.9 2 12 1.23553 0.45 0 -17.1114 -10 0 NaN 0 8.096440e-01 6 3.80 2.0716 -2.1 -0.084 12 1.23553 0.47 0 -18.8392 -12 0 NaN 0 6.200020e-01 7 5.70 2.5196 -3.8 -3.9 13 0.832228 0.47 0 -20.6793 -17 0 NaN 0 8.419070e-01 8 8.54 1.9952 -3.6 -7.5 13 0.832228 0.5 0 -22.5903 -21 0 NaN 0 6.341520e-01 9 12.81 1.589 -3 -10 12 0.832228 0.55 0 -23.0752 -23 0 NaN 0 5.179680e-01 10 19.22 1.0613 -1.5 -12 10 0.832228 0.62 0 -23.0752 -23 0 NaN 0 5.395840e-01 11 28.83 0.69458 -0.65 -13 9.7 0.832228 0.73 0 -23.0752 -23 0 NaN 0 6.755270e-01 12 43.25 0.49103 -0.32 -13 9.3 0.832228 0.83 0 -23.0752 -23 0 NaN 0 7.312430e-01 13 64.87 0.24192 -0.24 -13 9 0.832228 0.83 0 -23.0752 -23 0 NaN 0 6.665880e-01 14 97.31 0.21816 -0.55 -14 8.9 0.832228 0.83 0 -23.4777 -23 0 NaN 0 6.377100e-01 15 145.96 0.17636 -0.075 -14 8.8 0.832228 0.83 0 -23.4777 -23 0 NaN 0 9.225070e-01 16 218.95 0.097141 -0.048 -14 8.7 0.832228 0.83 0 -23.4777 -23 0 NaN 0 6.015810e-01 17 328.42 0.16623 -0.031 -14 8.7 0.832228 0.83 0 -23.4777 -23 0 NaN 0 6.092270e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation 1 178 1.00 0.02 121 1.29 0.03 greedyRandom 2 208 1.00 0.04 155 1.19 0.10 greedyRandom 3 335 1.00 0.07 249 1.20 0.22 greedyRandom 4 424 1.00 0.11 301 1.28 0.35 greedyRandom 5 393 1.00 0.20 272 1.30 0.42 greedyRandom 6 505 1.00 0.29 354 1.29 0.51 greedyRandom 7 280 1.00 0.34 208 1.28 0.54 greedyRandom 8 378 1.00 0.41 281 1.30 0.58 greedyRandom 9 395 1.00 0.45 291 1.31 0.63 greedyRandom 10 272 1.00 0.50 211 1.28 0.65 greedyRandom 11 281 1.00 0.54 203 1.36 0.67 greedyRandom 12 315 1.00 0.59 228 1.36 0.68 greedyRandom 13 189 1.00 0.61 134 1.29 0.69 greedyRandom 14 358 1.00 0.65 258 1.40 0.70 greedyRandom 15 299 1.00 0.68 218 1.28 0.72 greedyRandom 16 166 1.00 0.69 116 1.30 0.72 greedyRandom 17 275 1.00 0.70 191 1.45 0.73 greedyRandom 18 249 1.00 0.71 175 1.45 0.73 greedyRandom 19 150 1.00 0.72 107 1.38 0.74 greedyRandom 20 235 1.00 0.73 176 1.36 0.74 greedyRandom iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation thermoKernel terminating early: n = nMax = 20
--- thermoKernel END ----
[nMet,nRxn]=size(tissueModel.S)
nMet = 2203
nRxn = 2986
Compare the target versus predicted model
%plotThermoCoreStats(activeInactiveRxn, presentAbsentMet, thermoModelMetBool, thermoModelRxnBool);
plotThermoKernelExtractStats(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights, thermoModelMetBool, thermoModelRxnBool)

Save weights

rxnWeightsTmp=rxnWeights;
metWeightsTmp=metWeights;

Submodel with just metabolites specified

metWeights=metWeightsTmp;
rxnWeights(:)=0;
[tissueModel, thermoModelMetBool, thermoModelRxnBool] = thermoKernel(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights,param);
--- thermoKernel START ---- thermoKernel parameters: printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 n: 200 normalizeZeroNormWeights: 0 formulation: 'pqzwrs' epsilon: 1.0000e-06 removeOrphanGenes: 1 nbMaxIteration: 30 nMax: 20 iterationMethod: 'greedyRandom' warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. 0 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0.5 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 539 min cardinality variables: 0 mean(c(p)) 0 min(c(p)) 0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 1088 max cardinality variables: 0 mean(c(q)) 0 min(c(q)) 0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 36920 cardinality free variables: 0.46 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6648e+05 -8.4e+07 7.9 10 2.08529 0.087 0 -7.85827 -2.2 0 NaN 0 4.709090e-01 2 0.75 2.4394 -0.69 7.2 9.4 1.67053 0.13 0 -8.28783 -2.4 0 NaN 0 4.895420e-01 3 1.12 1.334 -0.47 6.7 9.3 1.23553 0.24 0 -8.28783 -2.8 0 NaN 0 4.830180e-01 4 1.69 3.7656 -1 5.7 12 1.23553 0.36 0 -8.28783 -6.2 0 NaN 0 5.411140e-01 5 2.53 3.2162 -1.7 4 11 1.23553 0.45 0 -8.73752 -7.9 0 NaN 0 4.433150e-01 6 3.80 2.1001 -1.8 2.1 11 1.23553 0.47 0 -10.0028 -9.1 0 NaN 0 4.702590e-01 7 5.70 1.3982 -1.2 0.95 10 1.23553 0.51 0 -10.0028 -9.8 0 NaN 0 5.574460e-01 8 8.54 0.78275 -1.5 -0.51 9.8 0.832228 0.5 0 -12.29 -11 0 NaN 0 5.730400e-01 9 12.81 0.52353 -1.8 -2.3 9.5 0.832228 0.55 0 -12.29 -12 0 NaN 0 5.270800e-01 10 19.22 0.37265 -0.29 -2.6 9.1 0.832228 0.62 0 -12.29 -12 0 NaN 0 4.930840e-01 11 28.83 0.33297 -0.11 -2.7 8.9 0.832228 0.73 0 -12.29 -12 0 NaN 0 5.038720e-01 12 43.25 0.19468 -0.5 -3.2 8.7 0.832228 0.83 0 -13.2195 -13 0 NaN 0 4.956000e-01 13 64.87 0.099656 -0.91 -4.1 8.7 0.832228 0.83 0 -15.0362 -14 0 NaN 0 7.131520e-01 14 97.31 0.11652 -1.5 -5.6 8.6 0.832228 0.83 0 -15.0362 -15 0 NaN 0 5.487360e-01 15 145.96 0.12539 -0.03 -5.6 8.6 0.832228 0.83 0 -15.0362 -15 0 NaN 0 6.832550e-01 16 218.95 0.067216 -0.018 -5.6 8.6 0.832228 0.83 0 -15.0362 -15 0 NaN 0 5.578550e-01 17 328.42 0.13964 -0.012 -5.6 8.6 0.832228 0.83 0 -15.0362 -15 0 NaN 0 5.004120e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. 0 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0.5 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 539 min cardinality variables: 0 mean(c(p)) 0 min(c(p)) 0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 1088 max cardinality variables: 0 mean(c(q)) 0 min(c(q)) 0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 36920 cardinality free variables: 0.46 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.7116e+05 -8.5e+07 7.9 10 2.08529 0.079 0 -7.12189 -2.2 0 NaN 0 4.543470e-01 2 0.75 2.4188 -0.68 7.2 9.5 1.23553 0.11 0 -7.55145 -2.4 0 NaN 0 5.571920e-01 3 1.12 2.6941 -0.72 6.5 10 1.23553 0.24 0 -8.35329 -3.7 0 NaN 0 4.783050e-01 4 1.69 4.9777 -0.99 5.5 11 1.23553 0.36 0 -8.35329 -6 0 NaN 0 5.466920e-01 5 2.53 4.4276 -1.7 3.8 11 1.23553 0.45 0 -8.80298 -7.9 0 NaN 0 4.477940e-01 6 3.80 2.3513 -2.2 1.6 11 1.23553 0.47 0 -10.5203 -9.6 0 NaN 0 4.773480e-01 7 5.70 1.6188 -1.2 0.4 10 1.23553 0.51 0 -10.5203 -10 0 NaN 0 5.048410e-01 8 8.54 0.63996 -0.59 -0.19 9.6 1.23553 0.56 0 -10.5203 -10 0 NaN 0 5.939500e-01 9 12.81 0.55001 -0.9 -1.1 9.4 0.832228 0.55 0 -12.3555 -11 0 NaN 0 6.417610e-01 10 19.22 0.4508 -1.5 -2.6 9.1 0.832228 0.62 0 -12.3555 -12 0 NaN 0 5.281920e-01 11 28.83 0.42169 -0.11 -2.7 8.9 0.832228 0.73 0 -12.3555 -12 0 NaN 0 5.437710e-01 12 43.25 0.24485 -0.032 -2.8 8.7 0.832228 0.83 0 -12.3555 -12 0 NaN 0 5.730080e-01 13 64.87 0.19244 -0.075 -2.9 8.7 0.832228 0.83 0 -12.3555 -12 0 NaN 0 5.207950e-01 14 97.31 0.10717 -0.042 -2.9 8.6 0.832228 0.83 0 -12.3555 -12 0 NaN 0 6.306880e-01 15 145.96 0.083129 -0.023 -2.9 8.6 0.832228 0.83 0 -12.3555 -12 0 NaN 0 5.274500e-01 16 218.95 0.070749 -0.015 -2.9 8.6 0.832228 0.83 0 -12.3555 -12 0 NaN 0 5.263160e-01 17 328.42 0.14835 -0.0098 -2.9 8.6 0.832228 0.83 0 -12.3555 -12 0 NaN 0 4.673910e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation 1 137 1.00 NaN 95 1.21 0.03 greedyRandom 2 186 1.00 NaN 139 1.17 0.08 greedyRandom 3 297 1.00 NaN 226 1.19 0.20 greedyRandom 4 467 1.00 NaN 349 1.30 0.37 greedyRandom 5 203 1.00 NaN 149 1.23 0.43 greedyRandom 6 282 1.00 NaN 213 1.24 0.51 greedyRandom 7 285 1.00 NaN 201 1.32 0.56 greedyRandom 8 213 1.00 NaN 156 1.30 0.60 greedyRandom 9 168 1.00 NaN 124 1.27 0.62 greedyRandom 10 132 1.00 NaN 97 1.20 0.63 greedyRandom 11 248 1.00 NaN 184 1.33 0.66 greedyRandom 12 111 1.00 NaN 78 1.22 0.66 greedyRandom 13 200 1.00 NaN 148 1.36 0.68 greedyRandom 14 112 1.00 NaN 77 1.25 0.68 greedyRandom 15 984 1.00 NaN 709 1.26 0.70 random 16 1076 1.00 NaN 789 1.24 0.74 random 17 1024 1.00 NaN 744 1.27 0.74 random 18 1000 1.00 NaN 726 1.26 0.75 random 19 974 1.00 NaN 732 1.22 0.75 random 20 965 1.00 NaN 718 1.24 0.76 random iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation thermoKernel terminating early: n = nMax = 20
--- thermoKernel END ----
[nMet,nRxn]=size(tissueModel.S)
nMet = 1870
nRxn = 2635
Compare the target versus predicted model
plotThermoKernelExtractStats(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights, thermoModelMetBool, thermoModelRxnBool)

Submodel with just reactions specified

rxnWeights=rxnWeightsTmp;
metWeights(:)=0;
[tissueModel, thermoModelMetBool, thermoModelRxnBool] = thermoKernel(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights,param);
--- thermoKernel START ---- thermoKernel parameters: printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 n: 200 normalizeZeroNormWeights: 0 formulation: 'pqzwrs' epsilon: 1.0000e-06 removeOrphanGenes: 1 nbMaxIteration: 30 nMax: 20 iterationMethod: 'greedyRandom' warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0.5 = max(g0), the local weight on zero-norm of internal reaction rate. 0 = min(h0), the local weight on zero-norm of metabolite production rate. 0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 854 min cardinality variables: -0 mean(c(p)) -0 min(c(p)) -0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 824 max cardinality variables: -0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 36869 cardinality free variables: 0.46 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6945e+05 -8.5e+07 8.4 8.4 0.408362 0.02 0 -3.08478 -0.07 0 NaN 0 6.376060e-01 2 0.75 0.082294 -0.032 8.4 8.4 0.408362 0.028 0 -3.50665 -0.11 0 NaN 0 5.303260e-01 3 1.12 0.25722 -0.046 8.3 8.5 0.408362 0.03 0 -3.50665 -0.18 0 NaN 0 4.786180e-01 4 1.69 0.41187 -0.093 8.2 8.5 0.422161 0.0076 0 -3.50665 -0.3 0 NaN 0 4.327600e-01 5 2.53 5.4552e-14 -0.15 8.1 8.5 0.422161 0.011 0 -3.50665 -0.45 0 NaN 0 4.586200e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0.5 = max(g0), the local weight on zero-norm of internal reaction rate. 0 = min(h0), the local weight on zero-norm of metabolite production rate. 0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 854 min cardinality variables: -0 mean(c(p)) -0 min(c(p)) -0 max(c(p)) 1 lambda0 0.4 min(k) 0.5 max(k) 0 lambda1 0 min(o(p)) 0 max(o(p)) 824 max cardinality variables: -0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 36869 cardinality free variables: 0.46 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6798e+05 -8.5e+07 8.4 8.4 0.830523 0.02 0 -3.50665 -0.07 0 NaN 0 4.300630e-01 2 0.75 0.11192 -0.032 8.4 8.4 0.408362 0.028 0 -3.50665 -0.11 0 NaN 0 4.447940e-01 3 1.12 0.25722 -0.046 8.3 8.5 0.408362 0.03 0 -3.50665 -0.18 0 NaN 0 5.177530e-01 4 1.69 0.41187 -0.093 8.2 8.5 0.422161 0.0076 0 -3.50665 -0.3 0 NaN 0 4.119660e-01 5 2.53 5.4552e-14 -0.15 8.1 8.5 0.422161 0.011 0 -3.50665 -0.45 0 NaN 0 6.149510e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation 1 100 1.00 0.01 72 1.15 NaN greedyRandom 2 99 1.00 0.01 70 1.21 NaN greedyRandom 3 113 1.00 0.02 79 1.15 NaN greedyRandom 4 113 1.00 0.03 83 1.12 NaN greedyRandom 5 181 1.00 0.07 134 1.14 NaN greedyRandom 6 240 1.00 0.14 183 1.14 NaN greedyRandom 7 192 1.00 0.20 128 1.38 NaN greedyRandom 8 184 1.00 0.24 130 1.31 NaN greedyRandom 9 276 1.00 0.30 192 1.38 NaN greedyRandom 10 267 1.00 0.35 193 1.28 NaN greedyRandom 11 309 1.00 0.42 227 1.29 NaN greedyRandom 12 293 1.00 0.49 209 1.44 NaN greedyRandom 13 295 1.00 0.54 205 1.39 NaN greedyRandom 14 214 1.00 0.57 173 1.25 NaN greedyRandom 15 340 1.00 0.62 249 1.37 NaN greedyRandom 16 284 1.00 0.66 203 1.37 NaN greedyRandom 17 177 1.00 0.67 122 1.42 NaN greedyRandom 18 116 1.00 0.67 87 1.20 NaN greedyRandom 19 147 1.00 0.68 101 1.40 NaN greedyRandom 20 199 1.00 0.70 139 1.42 NaN greedyRandom iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation thermoKernel terminating early: n = nMax = 20
--- thermoKernel END ----
 
[nMet,nRxn]=size(tissueModel.S)
nMet = 1441
nRxn = 1924
Compare the target versus predicted model
plotThermoKernelExtractStats(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeights, presentAbsentMet, metWeights, thermoModelMetBool, thermoModelRxnBool)
 

Submodel with just active metabolites specified

metWeightsRed=metWeightsTmp;
rxnWeightsRed=rxnWeightsTmp*0;
metWeightsRed(metWeightsRed>=0)=0;
[tissueModel, thermoModelMetBool, thermoModelRxnBool] = thermoKernel(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeightsRed, presentAbsentMet, metWeightsRed,param);
--- thermoKernel START ---- thermoKernel parameters: printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 n: 200 normalizeZeroNormWeights: 0 formulation: 'pqzwrs' epsilon: 1.0000e-06 removeOrphanGenes: 1 nbMaxIteration: 30 nMax: 20 iterationMethod: 'greedyRandom' warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. 0 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 0 min cardinality variables: NaN mean(c(p)) NaN min(c(p)) NaN max(c(p)) 1 lambda0 NaN min(k) NaN max(k) 0 lambda1 NaN min(o(p)) NaN max(o(p)) 1088 max cardinality variables: 0 mean(c(q)) 0 min(c(q)) 0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 37459 cardinality free variables: 0.45 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6587e+05 -8.5e+07 7.4 13 0 0 0 -14.7771 -5.9 0 NaN 0 5.243140e-01 2 0.75 5.8717 -2.6 4.8 14 0 0 0 -14.3228 -9 0 NaN 0 4.814740e-01 3 1.12 3.0027 -2.1 2.7 13 0 0 0 -15.2166 -10 0 NaN 0 4.797330e-01 4 1.69 5.6907 -2.2 0.55 13 0 0 0 -17.0067 -13 0 NaN 0 5.129960e-01 5 2.53 5.2054 -4.5 -3.9 14 0 0 0 -20.1726 -18 0 NaN 0 5.961270e-01 6 3.80 2.462 -3.4 -7.3 13 0 0 0 -21.5946 -20 0 NaN 0 6.517570e-01 7 5.70 1.606 -2.7 -10 12 0 0 0 -22.502 -22 0 NaN 0 6.102790e-01 8 8.54 1.0881 -1.7 -12 11 0 0 0 -22.502 -22 0 NaN 0 5.072370e-01 9 12.81 0.80985 -0.89 -13 10 0 0 0 -22.502 -23 0 NaN 0 4.409440e-01 10 19.22 0.49523 -0.58 -13 9.4 0 0 0 -22.502 -23 0 NaN 0 5.111380e-01 11 28.83 0.38668 -0.36 -13 9 0 0 0 -22.502 -23 0 NaN 0 5.276610e-01 12 43.25 0.28276 -0.38 -14 8.8 0 0 0 -23.4718 -23 0 NaN 0 5.855910e-01 13 64.87 0.21412 -0.94 -15 8.7 0 0 0 -23.4718 -23 0 NaN 0 5.123230e-01 14 97.31 0.16629 -0.083 -15 8.6 0 0 0 -23.4718 -23 0 NaN 0 5.380850e-01 15 145.96 0.084169 -0.26 -15 8.5 0 0 0 -25.2885 -24 0 NaN 0 6.800160e-01 16 218.95 0.12117 -1.6 -17 8.5 0 0 0 -25.2885 -25 0 NaN 0 5.618230e-01 17 328.42 0.16767 -0.14 -17 8.5 0 0 0 -25.7429 -25 0 NaN 0 6.236500e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. 0 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0.5 = min(h0), the local weight on zero-norm of metabolite production rate. 0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 0 min cardinality variables: NaN mean(c(p)) NaN min(c(p)) NaN max(c(p)) 1 lambda0 NaN min(k) NaN max(k) 0 lambda1 NaN min(o(p)) NaN max(o(p)) 1088 max cardinality variables: 0 mean(c(q)) 0 min(c(q)) 0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 37459 cardinality free variables: 0.45 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6895e+05 -8.5e+07 8.1 9 0 0 0 -6.13179 -0.87 0 NaN 0 4.684520e-01 2 0.75 1.691 -0.42 7.7 8.8 0 0 0 -6.56136 -1.1 0 NaN 0 3.991660e-01 3 1.12 2.8196 -0.42 7.3 10 0 0 0 -6.56136 -2.7 0 NaN 0 4.194870e-01 4 1.69 3.3125 -0.91 6.4 11 0 0 0 -6.56136 -4.4 0 NaN 0 4.419690e-01 5 2.53 3.3113 -2.3 4.1 11 0 0 0 -7.84509 -6.9 0 NaN 0 4.484160e-01 6 3.80 1.408 -1 3 11 0 0 0 -8.27635 -7.5 0 NaN 0 5.097870e-01 7 5.70 1.4232 -1.5 1.5 10 0 0 0 -8.72837 -8.5 0 NaN 0 4.473800e-01 8 8.54 0.75486 -0.82 0.71 9.6 0 0 0 -9.63582 -8.9 0 NaN 0 4.262520e-01 9 12.81 0.5787 -1.1 -0.43 9.2 0 0 0 -9.63582 -9.6 0 NaN 0 4.615360e-01 10 19.22 0.45877 -0.29 -0.72 8.9 0 0 0 -9.63582 -9.6 0 NaN 0 5.839770e-01 11 28.83 0.19759 -0.19 -0.91 8.7 0 0 0 -9.63582 -9.6 0 NaN 0 4.733700e-01 12 43.25 0.17445 -0.12 -1 8.6 0 0 0 -9.63582 -9.6 0 NaN 0 4.586240e-01 13 64.87 0.1615 -0.44 -1.5 8.5 0 0 0 -11.4525 -10 0 NaN 0 6.736550e-01 14 97.31 0.14097 -1.5 -2.9 8.5 0 0 0 -11.4525 -11 0 NaN 0 4.498600e-01 15 145.96 0.12588 -0.026 -3 8.5 0 0 0 -11.4525 -11 0 NaN 0 4.918090e-01 16 218.95 0.086138 -0.015 -3 8.5 0 0 0 -11.4525 -11 0 NaN 0 4.800920e-01 17 328.42 0.074785 -0.01 -3 8.5 0 0 0 -11.4525 -11 0 NaN 0 4.616160e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation 1 122 1.00 NaN 86 1.16 0.02 greedyRandom 2 230 1.00 NaN 172 1.16 0.10 greedyRandom 3 231 1.00 NaN 183 1.18 0.18 greedyRandom 4 352 1.00 NaN 263 1.22 0.31 greedyRandom 5 283 1.00 NaN 206 1.32 0.40 greedyRandom 6 274 1.00 NaN 199 1.37 0.47 greedyRandom 7 275 1.00 NaN 201 1.34 0.54 greedyRandom 8 175 1.00 NaN 139 1.19 0.56 greedyRandom 9 130 1.00 NaN 97 1.20 0.58 greedyRandom 10 307 1.00 NaN 235 1.34 0.62 greedyRandom 11 243 1.00 NaN 183 1.31 0.65 greedyRandom 12 120 1.00 NaN 84 1.29 0.65 greedyRandom 13 218 1.00 NaN 158 1.32 0.67 greedyRandom 14 113 1.00 NaN 80 1.24 0.67 greedyRandom 15 131 1.00 NaN 99 1.30 0.68 greedyRandom 16 129 1.00 NaN 92 1.30 0.69 greedyRandom 17 164 1.00 NaN 116 1.30 0.69 greedyRandom 18 166 1.00 NaN 123 1.26 0.70 greedyRandom 19 165 1.00 NaN 117 1.35 0.71 greedyRandom 20 141 1.00 NaN 95 1.36 0.72 greedyRandom iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation thermoKernel terminating early: n = nMax = 20
--- thermoKernel END ----
[nMet,nRxn]=size(tissueModel.S)
nMet = 1448
nRxn = 1753
Compare the target versus predicted model
plotThermoKernelExtractStats(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeightsRed, presentAbsentMet, metWeightsRed, thermoModelMetBool, thermoModelRxnBool)

Submodel with just active reactions specified

rxnWeightsRed=rxnWeightsTmp;
metWeightsRed=metWeightsTmp*0;
rxnWeightsRed(rxnWeightsRed>=0)=0;
[tissueModel, thermoModelMetBool, thermoModelRxnBool] = thermoKernel(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeightsRed, presentAbsentMet, metWeightsRed,param);
--- thermoKernel START ---- thermoKernel parameters: printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 n: 200 normalizeZeroNormWeights: 0 formulation: 'pqzwrs' epsilon: 1.0000e-06 removeOrphanGenes: 1 nbMaxIteration: 30 nMax: 20 iterationMethod: 'greedyRandom' warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0 = min(h0), the local weight on zero-norm of metabolite production rate. -0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 0 min cardinality variables: NaN mean(c(p)) NaN min(c(p)) NaN max(c(p)) 1 lambda0 NaN min(k) NaN max(k) 0 lambda1 NaN min(o(p)) NaN max(o(p)) 824 max cardinality variables: -0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 37723 cardinality free variables: 0.45 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6384e+05 -8.4e+07 8.4 8.5 0 0 0 -3.11866 -0.037 0 NaN 0 4.160000e-01 2 0.75 1.1692 -0.091 8.3 8.4 0 0 0 -3.54899 -0.1 0 NaN 0 5.045740e-01 3 1.12 0.15183 -0.062 8.3 8.4 0 0 0 -3.50665 -0.17 0 NaN 0 4.256620e-01 4 1.69 0.092443 -0.084 8.2 8.4 0 0 0 -3.50665 -0.26 0 NaN 0 4.319150e-01 5 2.53 0.44167 -0.14 8.1 8.5 0 0 0 -3.50665 -0.44 0 NaN 0 4.320100e-01 6 3.80 0.31473 -0.25 7.8 8.5 0 0 0 -3.95605 -0.74 0 NaN 0 4.248850e-01 7 5.70 1.8987 -0.76 7 9.6 0 0 0 -4.90384 -2.5 0 NaN 0 5.142230e-01 8 8.54 1.0204 -1.2 5.8 9.4 0 0 0 -5.39475 -3.5 0 NaN 0 4.927650e-01 9 12.81 1.4536 -1.1 4.8 10 0 0 0 -5.83297 -5.2 0 NaN 0 4.585050e-01 10 19.22 0.67982 -0.85 3.9 9.8 0 0 0 -5.83297 -5.8 0 NaN 0 4.072770e-01 11 28.83 0.60377 -0.54 3.4 9.2 0 0 0 -5.83297 -5.8 0 NaN 0 3.689780e-01 12 43.25 0.45629 -0.35 3 8.9 0 0 0 -5.83297 -5.8 0 NaN 0 4.147150e-01 13 64.87 0.57354 -0.22 2.8 8.7 0 0 0 -5.83297 -5.8 0 NaN 0 4.833740e-01 14 97.31 0.4286 -0.13 2.7 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.474080e-01 15 145.96 0.145 -0.036 2.7 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.640280e-01 16 218.95 0.11648 -0.019 2.6 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 3.695540e-01 17 328.42 0.083385 -0.012 2.6 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.384230e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). warmStartMethod: 'random' formulation: 'pqzwrs' thetaMultiplier: 1.5000 theta: 0.5000 regularizeOuter: 1 epsilon: 1.0000e-06 printLevel: 1 relaxBounds: 0 acceptRepairedFlux: 1 thermoConsistencyMethod: 'cycleFreeFlux' bigNum: 10000 debug: 0 optCardThermo objective data: 1 = beta, the global weight on one-norm of internal reaction rate. -0.5 = min(g0), the local weight on zero-norm of internal reaction rate. 0 = max(g0), the local weight on zero-norm of internal reaction rate. -0 = min(h0), the local weight on zero-norm of metabolite production rate. -0 = max(h0), the local weight on zero-norm of metabolite production rate. optimizeCardinality objective data: 0 min cardinality variables: NaN mean(c(p)) NaN min(c(p)) NaN max(c(p)) 1 lambda0 NaN min(k) NaN max(k) 0 lambda1 NaN min(o(p)) NaN max(o(p)) 824 max cardinality variables: -0 mean(c(q)) -0 min(c(q)) -0 max(c(q)) 1 delta0 0.4 min(d) 0.5 max(d) 0 delta1 0 min(o(q)) 0 max(o(q)) 37723 cardinality free variables: 0.45 mean(c(r)) -0 min(c(r)) 1 max(c(r)) 0 alpha1 0 min(o(r)) 0 max(o(r)) itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec 1 0.50 9.6834e+05 -8.4e+07 8.4 8.5 0 0 0 -3.54053 -0.037 0 NaN 0 4.197890e-01 2 0.75 1.251 -0.098 8.3 8.4 0 0 0 -3.97087 -0.11 0 NaN 0 4.360080e-01 3 1.12 0.091425 -0.056 8.3 8.4 0 0 0 -3.50665 -0.17 0 NaN 0 5.071690e-01 4 1.69 0.092443 -0.084 8.2 8.4 0 0 0 -3.50665 -0.26 0 NaN 0 4.041080e-01 5 2.53 0.44167 -0.14 8.1 8.5 0 0 0 -3.50665 -0.44 0 NaN 0 4.477540e-01 6 3.80 0.31473 -0.25 7.8 8.5 0 0 0 -3.95605 -0.74 0 NaN 0 4.119260e-01 7 5.70 1.8987 -0.76 7 9.6 0 0 0 -4.90384 -2.5 0 NaN 0 4.538870e-01 8 8.54 1.0204 -1.2 5.8 9.4 0 0 0 -5.39475 -3.5 0 NaN 0 5.025660e-01 9 12.81 1.4536 -1.1 4.8 10 0 0 0 -5.83297 -5.2 0 NaN 0 4.340600e-01 10 19.22 0.67982 -0.85 3.9 9.8 0 0 0 -5.83297 -5.8 0 NaN 0 4.349220e-01 11 28.83 0.60377 -0.54 3.4 9.2 0 0 0 -5.83297 -5.8 0 NaN 0 3.765820e-01 12 43.25 0.45629 -0.35 3 8.9 0 0 0 -5.83297 -5.8 0 NaN 0 4.053390e-01 13 64.87 0.57354 -0.22 2.8 8.7 0 0 0 -5.83297 -5.8 0 NaN 0 4.645910e-01 14 97.31 0.4286 -0.13 2.7 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.914500e-01 15 145.96 0.145 -0.036 2.7 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.420920e-01 16 218.95 0.11648 -0.019 2.6 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 3.801850e-01 17 328.42 0.083385 -0.012 2.6 8.5 0 0 0 -5.83297 -5.8 0 NaN 0 4.669970e-01 itn theta ||dx|| del_obj obj linear ||x||0 a(x) ||x||1 ||y||0 a(y) ||y||1 c(x,y) ||z||1 sec Optimise cardinality reached the stopping criterion. Finished. 100.00% thermodynamically feasible internal fluxes (checked by cycleFreeFlux method). iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation 1 111 1.00 0.02 74 1.24 NaN greedyRandom 2 110 1.00 0.02 76 1.24 NaN greedyRandom 3 132 1.00 0.03 93 1.22 NaN greedyRandom 4 127 1.00 0.04 94 1.18 NaN greedyRandom 5 215 1.00 0.10 152 1.24 NaN greedyRandom 6 190 1.00 0.15 138 1.23 NaN greedyRandom 7 247 1.00 0.21 175 1.29 NaN greedyRandom 8 268 1.00 0.28 193 1.34 NaN greedyRandom 9 300 1.00 0.36 217 1.27 NaN greedyRandom 10 284 1.00 0.43 201 1.42 NaN greedyRandom 11 234 1.00 0.47 162 1.37 NaN greedyRandom 12 325 1.00 0.53 235 1.39 NaN greedyRandom 13 217 1.00 0.56 155 1.39 NaN greedyRandom 14 223 1.00 0.59 158 1.29 NaN greedyRandom 15 197 1.00 0.60 137 1.37 NaN greedyRandom 16 136 1.00 0.61 94 1.29 NaN greedyRandom 17 189 1.00 0.63 135 1.39 NaN greedyRandom 18 222 1.00 0.65 148 1.47 NaN greedyRandom 19 316 1.00 0.68 219 1.42 NaN greedyRandom 20 155 1.00 0.69 114 1.34 NaN greedyRandom iter. nz.flux.%it.feas.int.flux. %feas.inc.flux. nz.prod. %it.feas.nz.prod. %feas.inc.prod. formulation thermoKernel terminating early: n = nMax = 20
--- thermoKernel END ----
 
[nMet,nRxn]=size(tissueModel.S)
nMet = 1449
nRxn = 1915
Compare the target versus predicted model
plotThermoKernelExtractStats(stoichFluxThermoConsistModel, activeInactiveRxn, rxnWeightsRed, presentAbsentMet, metWeightsRed, thermoModelMetBool, thermoModelRxnBool)