Entropic flux balance analysis

Author: Ronan Fleming, NUI Galway, Leiden University

INTRODUCTION

PROCEDURE

Set the directory for results
basePath='~/work/sbgCloud';
resultsDirectory=[basePath '/programModelling/projects/thermoModel/results/entropicFBA'];
Set the model to load
if ~exist('model','var')
modelToUse='iDN';
%modelToUse='ecoli_core';
%modelToUse = 'fork';
modelToUse = 'single';
%modelToUse='recon3';
switch modelToUse
case 'single'
driver_singleModel
case 'fork'
driver_forkModel
 
case 'ecoli_core'
dataDirectory=[basePath '/data'];
load([dataDirectory '/models/published/e_coli_core.mat']);
if ~isfield(model,'SConsistentMetBool') || ~isfield(model,'SConsistentRxnBool')
massBalanceCheck=0;
[~, ~, ~, ~, ~, ~, model, ~] = findStoichConsistentSubset(model, massBalanceCheck, 0);
end
model.lb(ismember(model.rxns,'BIOMASS_Ecoli_core_w_GAM'))=0.1;
case 'recon3'
% genericModelName = 'Recon3DModel_301_thermoFeas_BBF.mat';
% inputFolder = ['~' filesep 'work' filesep 'sbgCloud' filesep 'programReconstruction' filesep 'projects' filesep 'exoMetDN' filesep 'data' filesep 'xomics'];
% load([inputFolder filesep genericModelName])
genericModelName = 'Recon3DModel_301_xomics_input.mat';
inputFolder = ['~' filesep 'work' filesep 'sbgCloud' filesep 'programExperimental' filesep 'projects' filesep 'xomics' filesep 'data' filesep 'Recon3D_301'];
load([inputFolder filesep genericModelName])
%load([dataDirectory '/models/unpublished/2016_07_13_Recon3_SConsistencyCheck.mat']);
%model=modelRecon3model;
 
%change onto uMol/gDW/hr
model.lb = model.lb;
model.ub = model.ub;
case 'iDN'
%load('~/work/sbgCloud/programReconstruction/projects/exoMetDN/results/codeResults/iDN1/iDopaNeuro1/iDopaNeuro1.mat')
load('~/work/sbgCloud/programExperimental/projects/tracerBased/data/models/iDopaNeuro1.mat')
model = iDopaNeuro1;
clear Model
clear modelGenerationReport
%problematic reactions
% {'34DHPHEt'} {'3, 4-Dihydroxy-L-Phenylalanine Transport' } -0.53635 0 {'34dhphe[e] <=> 34dhphe[c] '}
% {'URAt' } {'Uracil Transport via Faciliated Diffusion'} -0.0013626 0 {'ura[e] <=> ura[c] ' }
model.lb(ismember(model.rxns,'34DHPHEt'))=-1e4;
model.ub(ismember(model.rxns,'34DHPHEt'))= 1e4;
model.lb(ismember(model.rxns,'URAt'))=-1e4;
model.ub(ismember(model.rxns,'URAt'))= 1e4;
end
end
Examine the constraints on the model
printConstraints(model,-1e4,1e4,model.SIntRxnBool)
Reversible_Reaction Name lb ub equation ___________________ ______ ___ __ ___________________ {'AB'} {'AB'} -20 20 {'A[c] <=> B[c] '}
printConstraints(model,-1e4,1e4,~model.SIntRxnBool)
printCouplingConstraints(model,1);
Set the data parameters for the entropic flux balance analysis problem
Linear optimisation direction
model.osenseStr = 'min';
Reaction parameters
model.cf=0.1;
model.cr=0.1;
model.g=2;
if 0
quad = zeros(size(model.S,2),1);
quad(~model.SIntRxnBool)=1;
model.Q=diag(quad);
end

Reaction parameters

Setting the bounds on reaction direction is important, because certain bounds may be incompatible with the existence of a thermodynamically feasible steady state
param.internalNetFluxBounds:
'original' maintains direction and magnitude of net flux from model.lb & model.ub
'directional' maintains direction of net flux from model.lb & model.ub but not magnitude
'random' random net flux direction, replacing constraints from model.lb & model.ub
%param.internalNetFluxBounds='directional';
param.internalNetFluxBounds='original';
%param.internalNetFluxBounds='none';
Maximum permitted internal flux
param.maxUnidirectionalFlux=100;
%model.vfl = [0;0;0]

Method parameters

param.method='fluxes';
param.method='fluxConc';
%param.method='fluxConcNorm';

Metabolite parameters (only relevant if a concentration method chosen)

Standard chemical potential
model.u0=0.2;
Strictly positive weight on concentration entropy maximisation
model.f=1;
Scalar maximum permitted metabolite concentration
param.maxConc=1;

Solver parameters

Set the solver
%param.solver='pdco';
param.solver='mosek';
if strcmp(param.solver,'mosek')
%set default mosek parameters for this type of problem
param=mosekParamSetEFBA(param);
end
Minimum permitted internal flux required for exponential cone reformulation
if strcmp(param.solver,'mosek')
param.minUnidirectionalFlux=0;
end
Tolerances
param.feasTol=1e-7;
Printing level
param.printLevel=2;
param.debug = 1;
Run entropic flux balance analysis
[solution,modelOut] = entropicFluxBalanceAnalysis(model,param);
Input model is feasible according to optimizeCbModel. Using existing internal net flux bounds without modification. MOSEK Version 9.2.37 (Build date: 2021-2-8 10:03:08) Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com Platform: Linux/64-X86 Problem Name : Objective sense : min Type : LO (linear optimization problem) Constraints : 5 Cones : 0 Scalar variables : 6 Matrix variables : 0 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.01 Problem Name : Objective sense : min Type : LO (linear optimization problem) Constraints : 5 Cones : 0 Scalar variables : 6 Matrix variables : 0 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 2 Optimizer - Cones : 0 Optimizer - Scalar variables : 6 conic : 0 Optimizer - Semi-definite variables: 0 scalarized : 0 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 3 after factor : 3 Factor - dense dim. : 0 flops : 2.10e+01 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 1.7e+03 5.0e+01 1.0e+02 0.00e+00 -7.000000000e+00 0.000000000e+00 2.9e+02 0.01 1 7.0e+01 2.1e+00 4.3e+00 -9.59e-01 -6.931219434e+01 -4.974403441e+01 1.2e+01 0.02 2 1.8e+01 5.4e-01 1.1e+00 1.52e-01 -3.096425733e+02 -3.059711217e+02 3.2e+00 0.02 3 1.3e+00 3.8e-02 7.8e-02 1.15e+00 -3.772348253e+02 -3.770527078e+02 2.2e-01 0.02 4 2.3e-02 6.7e-04 1.4e-03 9.77e-01 -3.831095240e+02 -3.831046437e+02 3.9e-03 0.02 5 2.4e-06 7.0e-08 1.5e-07 1.00e+00 -3.831999904e+02 -3.831999899e+02 4.1e-07 0.02 6 2.3e-10 7.2e-12 6.0e-11 1.00e+00 -3.832000000e+02 -3.832000000e+02 4.1e-11 0.02 Basis identification started. Primal basis identification phase started. Primal basis identification phase terminated. Time: 0.00 Dual basis identification phase started. Dual basis identification phase terminated. Time: 0.00 Basis identification terminated. Time: 0.00 Optimizer terminated. Time: 0.03 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -3.8320000000e+02 nrm: 1e+02 Viol. con: 0e+00 var: 0e+00 Dual. obj: -3.8320000000e+02 nrm: 2e+00 Viol. con: 0e+00 var: 0e+00 Basic solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -3.8320000000e+02 nrm: 1e+02 Viol. con: 0e+00 var: 0e+00 Dual. obj: -3.8320000000e+02 nrm: 2e+00 Viol. con: 0e+00 var: 0e+00 Optimizer summary Optimizer - time: 0.03 Interior-point - iterations : 6 time: 0.02 Basis identification - time: 0.00 Primal - iterations : 0 time: 0.00 Dual - iterations : 2 time: 0.00 Clean primal - iterations : 0 time: 0.00 Clean dual - iterations : 0 time: 0.00 Simplex - time: 0.00 Primal simplex - iterations : 0 time: 0.00 Dual simplex - iterations : 0 time: 0.00 Mixed integer - relaxations: 0 time: 0.00 solveCobraEP: LP part of EPproblem is feasible according to solveCobraLP with mosek. MOSEK Version 9.2.37 (Build date: 2021-2-8 10:03:08) Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com Platform: Linux/64-X86 MOSEK Version 9.2.37 (Build date: 2021-2-8 10:03:08) Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com Platform: Linux/64-X86 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 23 Cones : 6 Scalar variables : 31 Matrix variables : 0 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 1 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 23 Cones : 6 Scalar variables : 31 Matrix variables : 0 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 3 Optimizer - Cones : 6 Optimizer - Scalar variables : 19 conic : 18 Optimizer - Semi-definite variables: 0 scalarized : 0 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 6 after factor : 6 Factor - dense dim. : 0 flops : 8.00e+01 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 9.8e+01 2.7e+00 9.7e+00 0.00e+00 9.869931814e-01 -7.745566259e+00 1.0e+00 0.01 1 1.9e+01 5.1e-01 4.0e+00 -9.64e-01 3.568212170e+01 3.144143894e+01 1.9e-01 0.02 2 8.2e+00 2.3e-01 2.0e+00 -6.80e-01 4.332290988e+01 4.273971862e+01 8.3e-02 0.02 3 3.2e+00 8.9e-02 4.6e-01 1.73e-01 1.822389114e+00 9.745048948e-01 3.3e-02 0.02 4 2.4e-01 6.7e-03 7.4e-03 7.74e-01 -6.336444307e+00 -6.473740492e+00 2.5e-03 0.02 5 3.5e-02 9.6e-04 4.2e-04 9.95e-01 -6.866398497e+00 -6.885043239e+00 3.6e-04 0.02 6 5.5e-03 1.5e-04 2.8e-05 9.99e-01 -7.050420803e+00 -7.053109702e+00 5.6e-05 0.02 7 1.5e-03 4.2e-05 4.3e-06 1.00e+00 -7.075586361e+00 -7.076259038e+00 1.6e-05 0.02 8 1.4e-05 3.9e-07 5.1e-09 1.00e+00 -7.079709735e+00 -7.079709825e+00 1.4e-07 0.02 9 1.6e-07 4.3e-09 5.9e-12 1.00e+00 -7.079838540e+00 -7.079838541e+00 1.6e-09 0.02 10 3.0e-09 8.2e-11 1.2e-14 1.01e+00 -7.079840671e+00 -7.079840672e+00 2.8e-11 0.02 11 7.5e-10 2.1e-11 1.6e-15 1.00e+00 -7.079840723e+00 -7.079840723e+00 1.6e-11 0.02 12 6.5e-10 1.8e-11 1.3e-15 1.00e+00 -7.079840703e+00 -7.079840703e+00 6.5e-12 0.02 13 6.4e-10 1.8e-11 1.3e-15 1.00e+00 -7.079840702e+00 -7.079840702e+00 6.0e-12 0.02 14 3.8e-10 5.7e-14 2.4e-19 1.00e+00 -7.079840712e+00 -7.079840712e+00 5.3e-13 0.02 15 3.8e-10 5.7e-14 2.4e-19 1.00e+00 -7.079840712e+00 -7.079840712e+00 5.3e-13 0.02 16 3.8e-10 5.7e-14 2.4e-19 1.00e+00 -7.079840712e+00 -7.079840712e+00 5.3e-13 0.03 Optimizer terminated. Time: 0.04 Interior-point solution summary Problem status : PRIMAL_AND_DUAL_FEASIBLE Solution status : OPTIMAL Primal. obj: -7.0798407116e+00 nrm: 1e+00 Viol. con: 0e+00 var: 0e+00 cones: 2e-09 Dual. obj: -7.0798407116e+00 nrm: 7e+00 Viol. con: 1e-13 var: 2e-13 cones: 0e+00 Optimizer summary Optimizer - time: 0.04 Interior-point - iterations : 17 time: 0.03 Basis identification - time: 0.00 Primal - iterations : 0 time: 0.00 Dual - iterations : 0 time: 0.00 Clean primal - iterations : 0 time: 0.00 Clean dual - iterations : 0 time: 0.00 Simplex - time: 0.00 Primal simplex - iterations : 0 time: 0.00 Dual simplex - iterations : 0 time: 0.00 Mixed integer - relaxations: 0 time: 0.00 > [mosek] Primal optimality condition in solveCobraEP satisfied. > [mosek] Dual optimality condition in solveCobraEP satisfied. Optimality conditions (biochemistry) 0 || N*(vf - vr) + B*ve - x + x0 - b ||_inf 1.7e-13 || k_vf - g + cf + ci + N'*y_N + y_v - z_vf + y_vt ||_inf 1.7e-13 || k_vr - g + cr - ci - N'*y_N - y_v - z_vr + y_vt ||_inf 3.2e-13 || k_x + u0 - y_N + z_dx + z_x + y_xt ||_inf 3.2e-13 || k_x0 + u0 + y_N - z_dx + z_x0 + y_x0t ||_inf 0 || k_1 + z_t ||_inf 2.4e-13 || k_e_vf - g - z_t_f ||_inf 2.4e-13 || k_e_vr - g - z_t_r ||_inf 3.5e-14 || k_tx - f - z_t_x ||_inf 3.5e-14 || k_tx0 - f - z_t_x0 ||_inf 2.7e-09 || t_f + vf*log(vf/(1'*(vf + vr))) ||_inf 2.7e-09 || t_r + vr*log(vr/(1'*(vf + vr))) ||_inf 3e-09 || t_x + x*log(x/(1'*x)) ||_inf 3e-09 || t_x0 + x0*log(x0/(1'*x0)) ||_inf Derived optimality conditions (fluxes) 1.6e-07 || k_vf - g.*log(vf) - g ||_inf 1.6e-07 || k_vr - g.*log(vr) - g ||_inf 1.6e-07 || g.*log(vf) + cf + ci + N'*y_N - z_vf ||_inf 1.6e-07 || g.*log(vr) + cr - ci - N'*y_N + - z_vr ||_inf Effects of internal bounds on net fluxes 0 || y_v ||_inf Effects of internal bounds on forward fluxes 1.9e-13 || z_vf ||_inf Effects of internal bounds on reverse fluxes 1.9e-13 || z_vr ||_inf Derived optimality conditions (concentrations) 1.5e-06 || sx - f.*log( x/ (1'*x)) - f ||_inf 1.5e-06 || sx0 - f.*log(x0/(1'*x0)) - f ||_inf 1.5e-06 || f.*log(x/(1'*x)) + f + u0 - y_N + z_dx + z_x + y_vt ||_inf 1.5e-06 || f.*log(x0/(1'*x0)) + f + u0 + y_N - z_dx + z_x0 + y_xt ||_inf 1.5e-06 || f.*log(x/(1'*x)) + u0 - y_N + z_dx + z_x ||_inf 1.5e-06 || f.*log(x0/(1'*x0)) + u0 + y_N - z_dx + z_x0 ||_inf Thermo conditions (fluxes) 0 || g.*log(vr/vf) - 2*N'*y_N ||_inf 0 || g.*log(vr/vf) + cr - cf - 2*ci - 2*N'*y_N - 2*y_v ||_inf 0 || g.*log(vr/vf) + cr - cf - 2*ci - 2*N'*y_N - 2*y_v - z_vr + z_vf ||_inf Thermo conditions (concentrations) 0 || f.*(log(x/x0)) - 2*y_N + 2*z_dx + z_x - z_x0 ||_inf Effects of internal bounds on change in concentrations 0 || z_dx ||_inf Effects of internal bounds on concentrations 2.3e-12 || z_x ||_inf Effects of internal bounds on initial concentrations 2.3e-12 || z_x0 ||_inf Inf min(slack) Inf max(slack)
model=modelOut;
Extract the fields of the solution structure
if solution.stat ~=1
disp('No solution')
return
end
[v,vf,vr,vt,y_N,z_v,z_dx,z_vf,z_vr,stat]=...
deal(solution.v,solution.vf,solution.vr,solution.vt,solution.y_N,solution.z_v,solution.z_dx,solution.z_vf,solution.z_vr,solution.stat);
cf = model.cf;
cr = model.cr;
 
g=deal(model.g);
if isfield(solution,'x')
[x, x0, z_x, z_x0, z_dx] = deal(solution.x, solution.x0, solution.z_x, solution.z_x0, solution.z_dx);
[u0,f]=deal(model.u0,model.f);
dx = x - x0;
end
Extract the internal (N) and exchange (B) stoichiometric matrices
N = model.S(:,model.SConsistentRxnBool);
B = model.S(:,~model.SConsistentRxnBool);
Extract the matrix of coupling constraints, if present
if isfield(model,'C')
C = model.C(:,model.SConsistentRxnBool);
y_C = solution.y_C;
end
Extract the objective on the internal reaction rates
e = ones(length(vf),1);
ci = solution.osense*model.c(model.SConsistentRxnBool);
Compute the nullspace of the internal stoichiometric matrix
if 0
%get nullspace of N
[Z,rankS]=getNullSpace(N,0);
else
Z=[];
end
if param.printLevel>2 || 1
v_detailed=1;
else
v_detailed=0;
end
Biochemical optimality conditions: steady state
if isfield(solution,'x')
fprintf('%8.2g %s\n',norm(N*(vf - vr) - x + x0 - model.b,inf), '|| N*(vf - vr) - x + x0 - b ||_inf');
else
fprintf('%8.2g %s\n',norm(model.S*v - model.b,inf), '|| S*v - b ||_inf');
end
0 || N*(vf - vr) - x + x0 - b ||_inf
Optionally, display detailed properties of flux related solution vector
if v_detailed
fprintf('%8.2g %s\n',sum(vf)+sum(vr),'sum(vf)+sum(vr) >= 0');
fprintf('%8.2g %s\n',norm(g.*reallog(vf),inf),'|| g*log(vf) ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr),inf),'|| g*log(vr) ||_inf');
fprintf('%8.2g %s\n',norm(cr,inf),'|| cr ||_inf');
fprintf('%8.2g %s\n',norm(cf,inf),'|| cf ||_inf');
fprintf('%8.2g %s\n',norm(ci,inf),'|| ci ||_inf');
fprintf('%8.2g %s\n',norm(N'*y_N,inf),'|| N''*y_N ||_inf');
fprintf('%8.2g %s\n',norm(z_v,inf),'|| z_v ||_inf');
fprintf('%8.2g %s\n',norm(z_vf,inf),'|| z_vf ||_inf');
fprintf('%8.2g %s\n',norm(z_vr,inf),'|| z_vr ||_inf');
if isfield(model,'C')
fprintf('%8.2g %s\n',norm(C'*y_C,inf),'|| C''*y_C ||_inf');
end
end
1.9 sum(vf)+sum(vr) >= 0
0.1 || g*log(vf) ||_inf
0.1 || g*log(vr) ||_inf
0.1 || cr ||_inf
0.1 || cf ||_inf
0 || ci ||_inf
0 || N'*y_N ||_inf
0 || z_v ||_inf
1.9e-13 || z_vf ||_inf
1.9e-13 || z_vr ||_inf
if isfield(model,'C')
fprintf('%8.2g %s\n',norm(g.*reallog(vf) + ci + cf + N'*y_N + C'*y_C + z_v + z_vf,inf), '|| g*log(vf) + ci + cf + N''*y_N + C''*y_C + z_v + z_vf ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr) - ci + cr - N'*y_N - C'*y_C - z_v + z_vr,inf),'|| g*log(vr) - ci + cr - N''*y_N - C''*y_C - z_v + z_vr ||_inf');
fprintf('%8.2g %s\n',norm( C'*y_C + z_v + z_vf,inf), '|| + C''*y_C + z_v + z_vf ||_inf');
fprintf('%8.2g %s\n',norm( - C'*y_C - z_v + z_vr,inf),'|| - C''*y_C - z_v + z_vr ||_inf');
else
fprintf('%8.2g %s\n',norm(g.*reallog(vf) + ci + cf + N'*y_N + z_v - z_vf,inf), '|| g*log(vf) + ci + cf + N''*y_N + z_v - z_vf ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr) - ci + cr - N'*y_N - z_v - z_vr,inf),'|| g*log(vr) - ci + cr - N''*y_N - z_v - z_vr ||_inf');
end
1.6e-07 || g*log(vf) + ci + cf + N'*y_N + z_v - z_vf ||_inf
1.6e-07 || g*log(vr) - ci + cr - N'*y_N - z_v - z_vr ||_inf
Derived biochemical optimality conditions (fluxes)
if isfield(model,'C')
res = g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*y_N - 2*C'*y_C - 2*z_v - z_vr + z_vf;
fprintf('%8.2g %s\n',norm(res,inf),'|| g*log(vr/vf) + cr - cf - 2*ci - 2*N''*y_N - 2*C''*y_C - 2*z_v - z_vr + z_vf ||_inf');
else
res = g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*y_N - 2*z_v - z_vr + z_vf;
fprintf('%8.2g %s\n',norm(res,inf),'|| g*log(vr/vf) + cr - cf - 2*ci - 2*N''*y_N - 2*z_v - z_vr + z_vf ||_inf');
end
0 || g*log(vr/vf) + cr - cf - 2*ci - 2*N'*y_N - 2*z_v - z_vr + z_vf ||_inf
Optionally, display detailed properties of concentration related solution vector
if isfield(solution,'x')
fprintf('\n%s\n','Optimality conditions (concentrations)')
if v_detailed
fprintf('%8.2g %s\n',norm(reallog(x),inf), '|| log(x)||_inf');
fprintf('%8.2g %s\n',norm(u0,inf), '|| u0 ||_inf');
fprintf('%8.2g %s\n',norm(y_N,inf), '|| y_N ||_inf');
fprintf('%8.2g %s\n',norm(z_dx,inf), '|| z_dx ||_inf');
fprintf('%8.2g %s\n',norm(z_x,inf), '|| z_x ||_inf');
fprintf('%8.2g %s\n',norm(z_x0,inf), '|| z_x0 ||_inf');
end
fprintf('%8.2g %s\n',norm(f.*reallog(x) + u0 - y_N + z_dx + z_x,inf), '|| f.*log(x) + u0 - y_N + z_dx + z_x ||_inf');
fprintf('%8.2g %s\n',norm(f.*reallog(x0) + u0 + y_N - z_dx + z_x0,inf),'|| f.*log(x0) + u0 + y_N - z_dx + z_x0 ||_inf');
 
fprintf('\n%s\n','Derived biochemical optimality conditions (concentrations)')
fprintf('%8.2g %s\n',norm(f.*reallog(x./x0) - 2*y_N + 2*z_dx + z_x - z_x0,inf),'|| f.*log(x/x0) - 2*y_N + 2*z_dx + z_x - z_x0 ||_inf');
fprintf('%8.2g %s\n',norm(f.*reallog(x.*x0) + 2*u0 + z_x + z_x0,inf),'|| f.*log(x.*x0) + 2*u0 + z_x + z_x0 ||_inf');
 
 
fprintf('\n%s\n','Derived biochemical optimality conditions (fluxes and concentrations)')
if isfield(model,'C')
fprintf('%8.2g %s\n',norm(g.*reallog(vf) + cf + ci + N'*(u0 + f.*reallog(x) + z_dx + z_x) + C'*y_C + z_v + z_vf,inf),'|| g*log(vf) + cf + ci + N''*(u0 + f*log(x) + z_dx + z_x) + C''*y_C + z_v + z_vf ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr) + cr - ci - N'*(u0 + f.*reallog(x) + z_dx + z_x) - C'*y_C - z_v + z_vr,inf), '|| g*log(vr) + cr - ci - N''*(u0 + f*log(x) + z_dx + z_x) - C''*y_C - z_v + z_vr ||_inf');
else
fprintf('%8.2g %s\n',norm(g.*reallog(vf) + cf + ci + N'*(u0 + f.*reallog(x) + z_dx + z_x) + z_v + z_vf,inf),'|| g*log(vf) + cf + ci + N''*(u0 + f*log(x) + z_dx + z_x) + z_v + z_vf ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr) + cr - ci - N'*(u0 + f.*reallog(x) + z_dx + z_x) - z_v + z_vr,inf),'|| g*log(vr) + cr - ci - N''*(u0 + f*log(x) + z_dx + z_x) - z_v + z_vr ||_inf');
end
 
fprintf('\n%s\n','Derived biochemical optimality conditions (fluxes and concentrations)')
if isfield(model,'C')
fprintf('%8.2g %s\n',norm(g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*(u0 + f.*reallog(x) + z_dx + z_x) - 2*C'*y_C - 2*z_v - z_vf + z_vr,inf), '|| g.*reallog(vr./vf) + cr - cf - 2*(ci + N''*(u0 + f*log(x) + z_dx + z_x) + C''*y_C + z_v) - z_vf + z_vr ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr./vf) + cr - cf - 2*ci + 2*N'*(u0 + f.*reallog(x0) - z_dx + z_x0) - 2*C'*y_C - 2*z_v - z_vf + z_vr,inf),'|| g.*reallog(vr./vf) + cr - cf - 2*(ci - N''*(u0 + f*log(x0) - z_dx + z_x0) + C''*y_C + z_v) - z_vf + z_vr ||_inf');
else
fprintf('%8.2g %s\n',norm(g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*(u0 + f.*reallog(x) + z_dx + z_x) - 2*z_v - z_vf + z_vr,inf), '|| g.*reallog(vr./vf) + cr - cf - 2*(ci + N''*(u0 + f*log(x) + z_dx + z_x) + z_v) - z_vf + z_vr ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr./vf) + cr - cf - 2*ci + 2*N'*(u0 + f.*reallog(x0) - z_dx + z_x0) - 2*z_v - z_vf + z_vr,inf),'|| g.*reallog(vr./vf) + cr - cf - 2*(ci - N''*(u0 + f*log(x0) - z_dx + z_x0) + z_v) - z_vf + z_vr ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr.*vf) + cr + cf - N'*(u0 + f.*reallog(x) + z_dx + z_x) - 2*z_v - z_vf + z_vr,inf), '|| g.*reallog(vr./vf) + cr - cf - 2*(ci + N''*(u0 + f*log(x) + z_dx + z_x) + z_v) - z_vf + z_vr ||_inf');
fprintf('%8.2g %s\n',norm(g.*reallog(vr./vf) + cr - cf - 2*ci + 2*N'*(u0 + f.*reallog(x0) - z_dx + z_x0) - 2*z_v - z_vf + z_vr,inf),'|| g.*reallog(vr./vf) + cr - cf - 2*(ci - N''*(u0 + f*log(x0) - z_dx + z_x0) + z_v) - z_vf + z_vr ||_inf');
end
end
Optimality conditions (concentrations)
0.2 || log(x)||_inf
0.2 || u0 ||_inf
0 || y_N ||_inf
0 || z_dx ||_inf
2.3e-12 || z_x ||_inf
2.3e-12 || z_x0 ||_inf
1.5e-06 || f.*log(x) + u0 - y_N + z_dx + z_x ||_inf
1.5e-06 || f.*log(x0) + u0 + y_N - z_dx + z_x0 ||_inf
Derived biochemical optimality conditions (concentrations)
0 || f.*log(x/x0) - 2*y_N + 2*z_dx + z_x - z_x0 ||_inf
3e-06 || f.*log(x.*x0) + 2*u0 + z_x + z_x0 ||_inf
fprintf('\n')
if isfield(solution,'v')
vfpad = NaN*ones(size(model.S,2),1);
vfpad(model.SConsistentRxnBool)=vf;
vrpad = NaN*ones(size(model.S,2),1);
vrpad(model.SConsistentRxnBool)=vr;
zvfpad = NaN*ones(size(model.S,2),1);
zvfpad(model.SConsistentRxnBool)=z_vr;
zvrpad = NaN*ones(size(model.S,2),1);
zvrpad(model.SConsistentRxnBool)=z_vr;
zvpad = NaN*ones(size(model.S,2),1);
zvpad(model.SConsistentRxnBool)=z_v;
T1 = table(model.lb, v,model.ub,zvpad,vfpad,zvfpad,vrpad,zvrpad,'VariableNames',{'vl','v','vu','z_v','vf','z_vf','vr','z_vr'});
if length(v)<=10
disp(T1)
end
end
vl v vu z_v vf z_vf vr z_vr ___ _ __ ___ _______ __________ _______ __________ -20 0 20 0 0.95123 -1.873e-13 0.95123 -1.873e-13
 
if isfield(solution,'x')
T2 = table(model.xl,solution.x,model.xu,z_x,model.x0l,x0,model.x0u,z_x0,model.dxl,dx,model.dxu,z_dx,...
'VariableNames',{'xl' 'x' 'xu' 'z_x' 'x0l' 'x0' 'x0u' 'z_x0' 'dxl' 'dx' 'dxu' 'z_dx'});
if length(solution.x)<=10
disp(T2)
end
end
xl x xu z_x x0l x0 x0u z_x0 dxl dx dxu z_dx __ _______ __ __________ ___ _______ ___ __________ ____ __ ___ ____ 0 0.81873 1 2.3487e-12 0 0.81873 1 2.3487e-12 -Inf 0 Inf 0 0 0.81873 1 2.3487e-12 0 0.81873 1 2.3487e-12 -Inf 0 Inf 0
%useful for debugging when tolerances not satifactory
l = model.lb(model.SConsistentRxnBool);
u = model.ub(model.SConsistentRxnBool);
vi = vf - vr;
logvrvf = reallog(vr./vf);
 
figure
histogram(log10(abs(z_v)))
title('distribution of abs(z_v)')
 
zBool = abs(z_v)>0.1;
if any(l(zBool)~=0)
figure
histogram(l(zBool))
title('distribution of internal lower bound when z_v substantial')
end
if any(u(zBool)~=inf)
figure
histogram(u(zBool))
title('distribution of internal upper bound when z_v substantial')
end
 
if any(zBool)
figure
histogram(res(zBool))
title('distribution of thermodynamic residual when z_v is substantial')
end
if any(zBool)
figure
plot(vi(zBool),logvrvf(zBool),'.')
xlabel('Small net flux')
ylabel('Change in chemical potential')
end
 
if any(~zBool)
%thermo residual is large even when z_v is small
figure
histogram(res(~zBool))
title('distribution of thermodynamic residual when z_v is small')
end
if any(~zBool)
figure
plot(vi(~zBool),logvrvf(~zBool),'.')
xlabel('Large net flux')
ylabel('Change in chemical potential')
end
rBool = abs(res)>0.001;
fwdBool = vi>0;
if any(rBool)
ab = z_vf - z_vr;
figure
histogram(ab(rBool))
title('distribution of z_vf - z_vr when residual is large')
end
if any(rBool)
%deviation from thermodynamic constraint occuring when net flux is
%substantial
figure
res1 = g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*y_N;
histogram(res1(rBool))
title('distribution of log(vr/vf) - 2*N''*y_N when residual is large')
end
if any(rBool)
figure
plot(vi(rBool),res(rBool),'.')
xlabel('Net flux')
ylabel('Large thermo residual')
end
if any(rBool)
figure
plot(vi(rBool & fwdBool),res(rBool & fwdBool),'.')
xlabel('Net flux')
ylabel('Large thermo residual')
end
if any(rBool)
figure
if isfield(model,'C')
res5 = g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*y_N - 2*C'*y_C - 2*z_v;
histogram(res5(rBool))
title('distribution of g.*log(vr/vf) + cr - cf - 2*ci - 2*N''*y_N - 2*C''*y_C - 2*z_v when residual is large')
else
res5 = g.*reallog(vr./vf) + cr - cf - 2*ci - 2*N'*y_N - 2*z_v;
histogram(res5(rBool))
title('distribution of g.*log(vr/vf) + cr - cf - 2*ci - 2*N''*y_N - 2*z_v when residual is large')
end
end
if any(~rBool)
figure
plot(res(~rBool),vi(~rBool),'.')
xlabel('Net flux')
ylabel('Small thermo residual')
end
if any(rBool)
figure
plot(vi(~rBool & fwdBool),res(~rBool & fwdBool),'.')
xlabel('Net flux')
ylabel('Small thermo residual')
end
if isequal(modelToUse,'single')
fprintf('%8.2g %s\n',norm(vf - exp(-(cf + cr - u0(1) - u0(2))/2)*x(1)), '|| vf - exp(-(cr + cf - u0 - u0)/2)*x1 ||')
end
0.046 || vf - exp(-(cr + cf - u0 - u0)/2)*x1 ||