Skip to content

Commit

Permalink
Fix: get flux target at different alphas (#342)
Browse files Browse the repository at this point in the history 
* fix: when to do protein usage minimization

* refactor: remove unnecesary inputs

* fix: alpha for WT flux

* chore: updateGECKOdoc

---------

Co-authored-by: Eduard Kerkhoven <[email protected]>
  • Loading branch information
@ae-tafur @edkerk
ae-tafur and edkerk authored Jul 17, 2023
1 parent 856872c commit f29852a
Show file tree
Hide file tree
Showing 4 changed files with 78 additions and 120 deletions.
6 changes: 3 additions & 3 deletions doc/src/geckomat/utilities/ecFSEOF.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../../matlabicon.gif)">
<li><a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter)">getFluxTarget</a> getFluxTarget</li></ul>
<li><a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)">getFluxTarget</a> getFluxTarget</li></ul>
This function is called by:
<ul style="list-style-image:url(../../../matlabicon.gif)">
</ul>
Expand Down Expand Up @@ -186,15 +186,15 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0098 tolerance = 1e-4;
0099
0100 <span class="comment">% Simulate WT (100% growth) and forced (X% growth and the rest towards product):</span>
0101 flux_WT = <a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter)">getFluxTarget</a>(ecModel,targetRxn,csRxn);
0101 flux_WT = <a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)">getFluxTarget</a>(ecModel,params.bioRxn,targetRxn,1);
0102
0103 <span class="comment">% Set values which are under solver tolerance</span>
0104 flux_WT(flux_WT &lt; 0 &amp; ~enzRxns) = 0;
0105 flux_WT(flux_WT &gt; 0 &amp; enzRxns) = 0;
0106
0107 progressbar(<span class="string">'Flux Scanning with Enforced Objective Function'</span>)
0108 <span class="keyword">for</span> i = 1:length(alpha)
0109 flux_MAX = <a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter)">getFluxTarget</a>(ecModel,targetRxn,csRxn,alpha(i));
0109 flux_MAX = <a href="getFluxTarget.html" class="code" title="function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)">getFluxTarget</a>(ecModel,params.bioRxn,targetRxn,alpha(i));
0110 <span class="comment">% Set values which are under solver tolerance</span>
0111 flux_MAX(flux_MAX &lt; 0 &amp; ~enzRxns) = 0;
0112 flux_MAX(flux_MAX &gt; 0 &amp; enzRxns) = 0;
Expand Down
129 changes: 53 additions & 76 deletions doc/src/geckomat/utilities/getFluxTarget.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ <h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../
<div class="box"><strong>getFluxTarget</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="box"><strong>function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter) </strong></div>
<div class="box"><strong>function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre class="comment"> getFluxTarget
Expand All @@ -37,15 +37,11 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="

Input:
model an ecModel in GECKO 3 format (with ecModel.ec structure).
bioRxn rxn ID for the biomass reaction.
targetRxn rxn ID for the production target reaction, a exchange
reaction is recommended.
csRxn rxn ID for the main carbon source uptake reaction.
alpha scalling factor for desired suboptimal growth.
(Optional, default 0.95)
tolerance numerical tolerance for fixing bounds
(Optional, default 1e-4)
modelAdapter a loaded model adapter. (Optional, will otherwise use
the default model adapter)

Output:
minFlux vector of minimum flux rates at minimum target production,
Expand All @@ -54,7 +50,7 @@ <h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="
corresponding to model.rxns

Usage:
[maxFlux, minFlux] = simulateGrowth(ecModel,targetRxn,csRxn,alpha,tol)</pre></div>
[maxFlux, minFlux] = getFluxTarget(ecModel,bioRxn,targetRxn,alpha)</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
Expand All @@ -69,7 +65,7 @@ <h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter)</a>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)</a>
0002 <span class="comment">% getFluxTarget</span>
0003 <span class="comment">%</span>
0004 <span class="comment">% Function that performs a series of LP optimizations on an ecModel,</span>
Expand All @@ -80,77 +76,58 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0009 <span class="comment">%</span>
0010 <span class="comment">% Input:</span>
0011 <span class="comment">% model an ecModel in GECKO 3 format (with ecModel.ec structure).</span>
0012 <span class="comment">% targetRxn rxn ID for the production target reaction, a exchange</span>
0013 <span class="comment">% reaction is recommended.</span>
0014 <span class="comment">% csRxn rxn ID for the main carbon source uptake reaction.</span>
0012 <span class="comment">% bioRxn rxn ID for the biomass reaction.</span>
0013 <span class="comment">% targetRxn rxn ID for the production target reaction, a exchange</span>
0014 <span class="comment">% reaction is recommended.</span>
0015 <span class="comment">% alpha scalling factor for desired suboptimal growth.</span>
0016 <span class="comment">% (Optional, default 0.95)</span>
0017 <span class="comment">% tolerance numerical tolerance for fixing bounds</span>
0018 <span class="comment">% (Optional, default 1e-4)</span>
0019 <span class="comment">% modelAdapter a loaded model adapter. (Optional, will otherwise use</span>
0020 <span class="comment">% the default model adapter)</span>
0021 <span class="comment">%</span>
0022 <span class="comment">% Output:</span>
0023 <span class="comment">% minFlux vector of minimum flux rates at minimum target production,</span>
0024 <span class="comment">% corresponding to model.rxns</span>
0025 <span class="comment">% maxFlux vector of maximum flux rates at maximum target production,</span>
0026 <span class="comment">% corresponding to model.rxns</span>
0027 <span class="comment">%</span>
0028 <span class="comment">% Usage:</span>
0029 <span class="comment">% [maxFlux, minFlux] = simulateGrowth(ecModel,targetRxn,csRxn,alpha,tol)</span>
0017 <span class="comment">%</span>
0018 <span class="comment">% Output:</span>
0019 <span class="comment">% minFlux vector of minimum flux rates at minimum target production,</span>
0020 <span class="comment">% corresponding to model.rxns</span>
0021 <span class="comment">% maxFlux vector of maximum flux rates at maximum target production,</span>
0022 <span class="comment">% corresponding to model.rxns</span>
0023 <span class="comment">%</span>
0024 <span class="comment">% Usage:</span>
0025 <span class="comment">% [maxFlux, minFlux] = getFluxTarget(ecModel,bioRxn,targetRxn,alpha)</span>
0026
0027 <span class="keyword">if</span> nargin &lt; 4 || isempty(alpha)
0028 alpha = 0.95;
0029 <span class="keyword">end</span>
0030
0031 <span class="keyword">if</span> nargin &lt; 6 || isempty(modelAdapter)
0032 modelAdapter = ModelAdapterManager.getDefault();
0033 <span class="keyword">if</span> isempty(modelAdapter)
0034 error(<span class="string">'Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.'</span>)
0035 <span class="keyword">end</span>
0036 <span class="keyword">end</span>
0037 params = modelAdapter.getParameters();
0038
0039 <span class="keyword">if</span> nargin &lt; 5 || isempty(tolerance)
0040 tolerance = 1e-4;
0041 <span class="keyword">end</span>
0042
0043 <span class="keyword">if</span> nargin &lt; 4 || isempty(alpha)
0044 alpha = 0.95;
0045 <span class="keyword">end</span>
0046
0047 <span class="comment">% Get relevant rxn indexes</span>
0048 targetRxnIdx = getIndexes(model, targetRxn,<span class="string">'rxns'</span>);
0049 csRxnIdx = getIndexes(model, csRxn,<span class="string">'rxns'</span>);
0050 bioRxnIdx = getIndexes(model, params.bioRxn,<span class="string">'rxns'</span>);
0051
0052 <span class="comment">% Maximize growth and fix carbon source and suboptimal growth</span>
0053 model = setParam(model, <span class="string">'obj'</span>, params.bioRxn, 1);
0054 sol = solveLP(model);
0055 model = setParam(model, <span class="string">'var'</span>, csRxn, sol.x(csRxnIdx), tolerance);
0056 model = setParam(model, <span class="string">'lb'</span>, params.bioRxn, sol.x(bioRxnIdx) * (1-tolerance) * alpha);
0057
0058 <span class="comment">% If minimum fluxes are required get them</span>
0059 minFlux = [];
0060 <span class="keyword">if</span> nargout == 2
0061 <span class="comment">% Minimize target</span>
0062 model = setParam(model, <span class="string">'obj'</span>, targetRxn, -1);
0063 sol = solveLP(model);
0064
0065 <span class="comment">% Now fix min value for target and minimize proteins usage</span>
0066 model.lb(targetRxnIdx) = sol.x(targetRxnIdx) * (1-tolerance);
0067 model = setParam(model, <span class="string">'obj'</span>, <span class="string">'prot_pool_exchange'</span>, 1);
0068 minSol = solveLP(model,1);
0069 minFlux = minSol.x;
0070 <span class="keyword">end</span>
0071
0072 <span class="comment">% Maximize target</span>
0073 model = setParam(model, <span class="string">'obj'</span>, targetRxn, 1);
0074 sol = solveLP(model);
0075
0076 <span class="comment">% Now fix max value for target and minimize proteins usage</span>
0077 model.lb(targetRxnIdx) = sol.x(targetRxnIdx) * (1-tolerance);
0078 model = setParam(model, <span class="string">'obj'</span>, <span class="string">'prot_pool_exchange'</span>, 1);
0079 maxSol = solveLP(model);
0080 maxFlux = maxSol.x;
0081
0082 <span class="keyword">end</span></pre></div>
0031 <span class="comment">% Fix a suboptimal alpha if equal to 1</span>
0032 <span class="keyword">if</span> alpha == 1
0033 alpha = 0.99;
0034 <span class="keyword">end</span>
0035
0036 <span class="comment">% Get relevant rxn indexes</span>
0037 bioRxnIdx = getIndexes(model, bioRxn,<span class="string">'rxns'</span>);
0038 poolIdx = strcmpi(model.rxns, <span class="string">'prot_pool_exchange'</span>);
0039
0040 <span class="comment">% Maximize growth and fix carbon source and suboptimal growth</span>
0041 model = setParam(model, <span class="string">'obj'</span>, bioRxn, 1);
0042 sol = solveLP(model);
0043 model = setParam(model, <span class="string">'lb'</span>, bioRxn, sol.x(bioRxnIdx) * alpha);
0044
0045 <span class="comment">% Minimize prot_pool_exchange and fix</span>
0046 model = setParam(model, <span class="string">'obj'</span>, <span class="string">'prot_pool_exchange'</span>, 1);
0047 sol = solveLP(model);
0048 model = setParam(model, <span class="string">'lb'</span>, <span class="string">'prot_pool_exchange'</span>, sol.x(poolIdx) * 1.01);
0049
0050 <span class="comment">% If minimum fluxes are required get them</span>
0051 minFlux = [];
0052 <span class="keyword">if</span> nargout == 2
0053 <span class="comment">% Minimize target</span>
0054 model = setParam(model, <span class="string">'obj'</span>, targetRxn, -1);
0055 minSol = solveLP(model);
0056 minFlux = minSol.x;
0057 <span class="keyword">end</span>
0058
0059 <span class="comment">% Maximize target</span>
0060 model = setParam(model, <span class="string">'obj'</span>, targetRxn, 1);
0061 maxSol = solveLP(model);
0062 maxFlux = maxSol.x;
0063 <span class="keyword">end</span></pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
4 changes: 2 additions & 2 deletions src/geckomat/utilities/ecFSEOF.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,15 +98,15 @@
tolerance = 1e-4;

% Simulate WT (100% growth) and forced (X% growth and the rest towards product):
flux_WT = getFluxTarget(ecModel,targetRxn,csRxn);
flux_WT = getFluxTarget(ecModel,params.bioRxn,targetRxn,1);

% Set values which are under solver tolerance
flux_WT(flux_WT < 0 & ~enzRxns) = 0;
flux_WT(flux_WT > 0 & enzRxns) = 0;

progressbar('Flux Scanning with Enforced Objective Function')
for i = 1:length(alpha)
flux_MAX = getFluxTarget(ecModel,targetRxn,csRxn,alpha(i));
flux_MAX = getFluxTarget(ecModel,params.bioRxn,targetRxn,alpha(i));
% Set values which are under solver tolerance
flux_MAX(flux_MAX < 0 & ~enzRxns) = 0;
flux_MAX(flux_MAX > 0 & enzRxns) = 0;
Expand Down
59 changes: 20 additions & 39 deletions src/geckomat/utilities/getFluxTarget.m
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [maxFlux, minFlux] = getFluxTarget(model,targetRxn,csRxn,alpha,tolerance,modelAdapter)
function [maxFlux, minFlux] = getFluxTarget(model,bioRxn,targetRxn,alpha)
% getFluxTarget
%
% Function that performs a series of LP optimizations on an ecModel,
Expand All @@ -9,15 +9,11 @@
%
% Input:
% model an ecModel in GECKO 3 format (with ecModel.ec structure).
% bioRxn rxn ID for the biomass reaction.
% targetRxn rxn ID for the production target reaction, a exchange
% reaction is recommended.
% csRxn rxn ID for the main carbon source uptake reaction.
% alpha scalling factor for desired suboptimal growth.
% (Optional, default 0.95)
% tolerance numerical tolerance for fixing bounds
% (Optional, default 1e-4)
% modelAdapter a loaded model adapter. (Optional, will otherwise use
% the default model adapter)
%
% Output:
% minFlux vector of minimum flux rates at minimum target production,
Expand All @@ -26,57 +22,42 @@
% corresponding to model.rxns
%
% Usage:
% [maxFlux, minFlux] = simulateGrowth(ecModel,targetRxn,csRxn,alpha,tol)

if nargin < 6 || isempty(modelAdapter)
modelAdapter = ModelAdapterManager.getDefault();
if isempty(modelAdapter)
error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
end
end
params = modelAdapter.getParameters();

if nargin < 5 || isempty(tolerance)
tolerance = 1e-4;
end
% [maxFlux, minFlux] = getFluxTarget(ecModel,bioRxn,targetRxn,alpha)

if nargin < 4 || isempty(alpha)
alpha = 0.95;
end

% Fix a suboptimal alpha if equal to 1
if alpha == 1
alpha = 0.99;
end

% Get relevant rxn indexes
targetRxnIdx = getIndexes(model, targetRxn,'rxns');
csRxnIdx = getIndexes(model, csRxn,'rxns');
bioRxnIdx = getIndexes(model, params.bioRxn,'rxns');
bioRxnIdx = getIndexes(model, bioRxn,'rxns');
poolIdx = strcmpi(model.rxns, 'prot_pool_exchange');

% Maximize growth and fix carbon source and suboptimal growth
model = setParam(model, 'obj', params.bioRxn, 1);
sol = solveLP(model);
model = setParam(model, 'var', csRxn, sol.x(csRxnIdx), tolerance);
model = setParam(model, 'lb', params.bioRxn, sol.x(bioRxnIdx) * (1-tolerance) * alpha);
% Maximize growth and fix carbon source and suboptimal growth
model = setParam(model, 'obj', bioRxn, 1);
sol = solveLP(model);
model = setParam(model, 'lb', bioRxn, sol.x(bioRxnIdx) * alpha);

% Minimize prot_pool_exchange and fix
model = setParam(model, 'obj', 'prot_pool_exchange', 1);
sol = solveLP(model);
model = setParam(model, 'lb', 'prot_pool_exchange', sol.x(poolIdx) * 1.01);

% If minimum fluxes are required get them
minFlux = [];
if nargout == 2
% Minimize target
model = setParam(model, 'obj', targetRxn, -1);
sol = solveLP(model);

% Now fix min value for target and minimize proteins usage
model.lb(targetRxnIdx) = sol.x(targetRxnIdx) * (1-tolerance);
model = setParam(model, 'obj', 'prot_pool_exchange', 1);
minSol = solveLP(model,1);
minSol = solveLP(model);
minFlux = minSol.x;
end

% Maximize target
model = setParam(model, 'obj', targetRxn, 1);
sol = solveLP(model);

% Now fix max value for target and minimize proteins usage
model.lb(targetRxnIdx) = sol.x(targetRxnIdx) * (1-tolerance);
model = setParam(model, 'obj', 'prot_pool_exchange', 1);
maxSol = solveLP(model);
maxFlux = maxSol.x;

end

0 comments on commit f29852a

Please sign in to comment.