From 22ee6516c02f177bf2f6170dfb76bfa1dfdb2d9a Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Tue, 4 Apr 2023 14:44:18 +0200
Subject: [PATCH 01/13] reverse: simulateGrowth for ecFSEOF

---
 .../ecYeastGEM/code => src/geckomat/utilities}/simulateGrowth.m   | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename {userData/ecYeastGEM/code => src/geckomat/utilities}/simulateGrowth.m (100%)

diff --git a/userData/ecYeastGEM/code/simulateGrowth.m b/src/geckomat/utilities/simulateGrowth.m
similarity index 100%
rename from userData/ecYeastGEM/code/simulateGrowth.m
rename to src/geckomat/utilities/simulateGrowth.m

From 85cc8affabcfdd31406949dbe56afae4488ce7db Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Tue, 4 Apr 2023 14:46:53 +0200
Subject: [PATCH 02/13] fix: consistency in terminology

---
 src/geckomat/utilities/{simulateGrowth.m => getFluxTarget.m} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename src/geckomat/utilities/{simulateGrowth.m => getFluxTarget.m} (100%)

diff --git a/src/geckomat/utilities/simulateGrowth.m b/src/geckomat/utilities/getFluxTarget.m
similarity index 100%
rename from src/geckomat/utilities/simulateGrowth.m
rename to src/geckomat/utilities/getFluxTarget.m

From 4e9fae17e856ff95facac0ffd9761c1ba9354569 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Tue, 4 Apr 2023 16:30:41 +0200
Subject: [PATCH 03/13] refactor: fluxes for a target

---
 src/geckomat/utilities/getFluxTarget.m | 135 ++++++++++++-------------
 1 file changed, 63 insertions(+), 72 deletions(-)

diff --git a/src/geckomat/utilities/getFluxTarget.m b/src/geckomat/utilities/getFluxTarget.m
index c6d0e41fb..0ab939975 100644
--- a/src/geckomat/utilities/getFluxTarget.m
+++ b/src/geckomat/utilities/getFluxTarget.m
@@ -1,83 +1,74 @@
-function flux = simulateGrowth(model,target,C_source,objCoeff,alpha,tol)
-%simulateGrowth
+function [minFlux, maxFlux] = getFluxTarget(ecModel,target,cSource,alpha,tolerance,modelAdapter)
+% getFluxTarget
 %
 %   Function that performs a series of LP optimizations on an ecModel,
 %   by first maximizing biomass, then fixing a suboptimal value and 
-%   proceeding to maximize a given production target reaction, last 
-%   a protein pool minimization is performed subject to the optimal 
-%   production level.
-%   
-%       model     (struct) ecModel with total protein pool constraint.
-%                 the model should come with growth pseudoreaction as 
-%                 an objective to maximize.
-%       target    (string) Rxn ID for the objective reaction.
-%       cSource   (string) Rxn ID for the main carbon source uptake 
-%                 reaction
-%       objCoeff  (integer) Coefficient for the target reaction in the
-%                 objective function
-%       alpha     (dobule) scalling factor for desired suboptimal growth
-%       tol       (double) numerical tolerance for fixing bounds
+%   proceeding to protein pool minimization, last a minimization and
+%   maximization of a given production target reaction is performed
+%   subject to the optimal production level.
 %
-% Usage: flux = simulateGrowth(model,target,C_source,alpha,tol)
+% Input:
+%   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
+%	target          rxn ID for the production target reaction, a exchange
+%                   reaction is recommended.
+%	cSource	        rxn ID for the main carbon source uptake reaction.
+%   alpha           scalling factor for desired suboptimal growth.
+%                   (Optional, defaul 1)
+%   tolerance       numerical tolerance for fixing bounds
+%                   (Optional, defaul 1E-8)
+%   modelAdapter    a loaded model adapter. (Optional, will otherwise use
+%                   the default model adapter)
 %
+% Output:
+%   minFlux         vector of minimum flux rates at minimum target production,
+%                   corresponding to model.rxns
+%   maxFlux         vector of maximum flux rates at maximum target production,
+%                   corresponding to model.rxns
+%
+% Usage:
+%   [minFlux, maxFlux] = simulateGrowth(model,target,C_source,alpha,tol)
 
-if nargin<6
-    tol = 1E-6;
-    if nargin<5
-        alpha = 1;
-        if nargin<4
-            objCoeff = 1;
-        end
-    end
-end
-%Fix a unit main carbon source uptake
-cSource_indx           = find(strcmpi(model.rxns,C_source));
-model.lb(cSource_indx) = (1-tol);
-model.ub(cSource_indx) = (1+tol);
-%Position of target rxn:
-posP      = strcmpi(model.rxns,target);
-growthPos = find(model.c);
-%Max growth:
-sol = solveLP(model);
-%Fix growth suboptimal and then max product:
-model.lb(growthPos) = sol.x(growthPos)*(1-tol)*alpha;
-flux                = optModel(model,posP,objCoeff,tol,sol.x);
-end
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-function sol = optModel(model,pos,c,tol,base_sol)
-%Find reversible reaction for the production target, if available 
-if any(pos)
-    rxn_code = model.rxns{pos};
-    if endsWith(rxn_code,'_REV')
-        rev_pos = strcmp(model.rxns,rxn_code(1:end-4));
-    else
-        rev_pos = strcmp(model.rxns,[rxn_code '_REV']);
+if nargin < 6 || isempty(modelAdapter)
+    modelAdapter = ModelAdapterManager.getDefaultAdapter();
+    if isempty(modelAdapter)
+        error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
     end
-    %Block reversible production target reaction
-    model.lb(rev_pos) = 0;
-    model.ub(rev_pos) = 0;
 end
-%Optimize for a given objective
-model.c      = zeros(size(model.rxns));
-model.c(pos) = c;
-sol          = solveLP(model);
-%If not successful then relax reversible rxn bounds with values
-%from the basal distribution
-if isempty(sol.x) && ~isempty(base_sol)
-    model.lb(rev_pos) = base_sol(rev_pos);
-    model.ub(rev_pos) = base_sol(rev_pos);
-    sol               = solveLP(model);
+params = modelAdapter.getParameters();
+
+if nargin < 5 || isempty(tolerance)
+    tolerance = 1E-8;
 end
-%Now fix optimal value for objective and minimize unmeasured proteins
-%usage
-if ~isempty(sol.x)
-    model.lb(pos)    = (1-tol)*sol.x(pos);
-    protPos          = contains(model.rxnNames,'prot_pool_');
-    model.c(:)       = 0;
-    model.c(protPos) = -1;
-    sol              = solveLP(model,1);
-    sol              = sol.x;
-else
-    sol = zeros(length(model.rxns),1);
+
+if nargin < 4 || isempty(alpha)
+    alpha = 1;
 end
+
+% Fix carbon source uptake
+cSourceIdx = strcmpi(ecModel.rxns,cSource);
+uptake = ecModel.lb(cSourceIdx);
+ecModel = setParam(ecModel, 'var', cSource, uptake, tolerance);
+
+% Max growth and fix growth suboptimal
+ecModel = setParam(ecModel, 'obj', params.bioRxn, 1);
+sol = solveLP(ecModel);
+bioRxnIdx = strcmpi(ecModel.rxns, params.bioRxn);
+ecModel = setParam(ecModel, 'lb', params.bioRxn, sol.x(bioRxnIdx) * (1-tolerance) * alpha);
+
+% Minimize prot_pool_exchange and fix
+poolIdx = strcmpi(ecModel.rxns, 'prot_pool_exchange');
+ecModel = setParam(ecModel, 'obj', 'prot_pool_exchange', 1);
+sol = solveLP(ecModel, 1);
+ecModel = setParam(ecModel, 'lb', params.bioRxn, sol.x(poolIdx));
+
+% Minimize target
+ecModel = setParam(ecModel, 'obj', target, -1);
+sol = solveLP(ecModel);
+minFlux = sol.x;
+
+% Maximize target
+ecModel = setParam(ecModel, 'obj', target, 1);
+sol = solveLP(ecModel);
+maxFlux = sol.x;
+
 end
\ No newline at end of file

From 855e213874de028b1103535413b300d58a1b39e0 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Tue, 4 Apr 2023 17:35:50 +0200
Subject: [PATCH 04/13] refactor: compatibility to GECKO3 ecFlux_scanning

---
 .../utilities/ecFSEOF/ecFlux_scanning.m       | 147 +++++++++---------
 1 file changed, 76 insertions(+), 71 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
index 030f53459..8ff2a81e2 100644
--- a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
+++ b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
@@ -1,130 +1,135 @@
-function FC = ecFlux_scanning(model,target,C_source,alpha,tol,filterG)
-%ecFlux_scanning
+function FC = ecFlux_scanning(ecModel,target,cSource,alpha,tolerance,filterG)
+% ecFlux_scanning
 %
-%       model     (struct) ecModel with total protein pool constraint.
-%                 the model should come with growth pseudoreaction as 
-%                 an objective to maximize.
-%       target    (string) Rxn ID for the production target reaction, 
-%                 a exchange reaction is recommended.
-%       cSource   (string) Rxn ID for the main carbon source uptake 
-%                 reaction
-%       alpha     (dobule) scalling factor for production yield 
-%                 for enforced objective limits
-%       tol       (double, optional) numerical tolerance for fixing 
-% 				  bounds
-%		filterG	  (logical) TRUE if genes K_scores results should
-%				  be filtered according to the alpha vector distribution
-%
-% Usage:  FC = compare_substrate(model,target,C_source,alpha,tol,filterG)
+% Input:
+%   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
+%	rxnTarget       rxn ID for the production target reaction, a exchange
+%                   reaction is recommended.
+%	cSource	        rxn ID for the main carbon source uptake reaction.
+%   alpha           scalling factor for production yield for enforced objective
+%                   limits
+%   tolerance       numerical tolerance for fixing bounds.
+%                   (Optional, defaul 1E-4)
+%   filterG	        logical value. TRUE if genes K_scores results should be
+%		            filtered according to the alpha vector distribution
+%                   (Optional, defaul false)
 %
+% Usage:
+%   FC = ecFlux_scanning(model,target,cSource,alpha,tolerance,filterG)
 
-if nargin<6
+if nargin < 6 || isempty(filerG)
 	filterG = false;
 end
 
-%Simulate WT (100% growth):
-cd ..
-FC.flux_WT = simulateGrowth(model,target,C_source,1,1);
-%Simulate forced (X% growth and the rest towards product) based on yield:
+if nargin < 5 || isempty(tolerance)
+    tolerance = 1E-4;
+end
+
+% Simulate WT (100% growth):
+[~, FC.flux_WT] = getFluxTarget(ecModel,target,cSource);
+
+% Simulate forced (X% growth and the rest towards product) based on yield:
 FC.alpha = alpha;
-%initialize fluxes and K_scores matrices
-v_matrix = zeros(length(model.rxns),length(alpha));
-k_matrix = zeros(length(model.rxns),length(alpha));
+
+% Initialize fluxes and K_scores matrices
+v_matrix = zeros(length(ecModel.rxns),length(alpha));
+k_matrix = zeros(length(ecModel.rxns),length(alpha));
 for i = 1:length(alpha)
     %disp(['Iteration #' num2str(i)])
-    FC.flux_MAX   = simulateGrowth(model,target,C_source,1,alpha(i));
+    [~, FC.flux_MAX] = getFluxTarget(ecModel,target,cSource,alpha(i));
     v_matrix(:,i) = FC.flux_MAX;
     k_matrix(:,i) = FC.flux_MAX./FC.flux_WT;
 end
-cd ecFSEOF
-%Take out rxns with no grRule:
-withGR   = ~cellfun(@isempty,model.grRules);
-%Generate rxn equations:
-rxnEqs   = constructEquations(model,model.rxns(withGR),true);
+
+% Take out rxns with no grRule:
+withGR   = ~cellfun(@isempty,ecModel.grRules);
+
+% Generate rxn equations:
+rxnEqs   = constructEquations(ecModel,ecModel.rxns(withGR),true);
 v_matrix = v_matrix(withGR,:);
 k_matrix = k_matrix(withGR,:);
-rxnGeneM = model.rxnGeneMat(withGR,:);
-FC.rxns  = [model.rxns(withGR),model.rxnNames(withGR),model.grRules(withGR),rxnEqs];
+rxnGeneM = ecModel.rxnGeneMat(withGR,:);
+FC.rxns  = [ecModel.rxns(withGR),ecModel.rxnNames(withGR),ecModel.grRules(withGR),rxnEqs];
 
-%Filter out rxns that are always zero -> k=0/0=NaN:
+% Filter out rxns that are always zero -> k=0/0=NaN:
 non_nan  = sum(~isnan(k_matrix),2) > 0;
 v_matrix = v_matrix(non_nan,:);
 k_matrix = k_matrix(non_nan,:);
 rxnGeneM = rxnGeneM(non_nan,:);
 FC.rxns  = FC.rxns(non_nan,:);
 
-%Replace remaining NaNs with 1s:
+% Replace remaining NaNs with 1s:
 k_matrix(isnan(k_matrix)) = 1;
-%Replace any Inf value with 1000 (maximum value is ~700):
+
+% Replace any Inf value with 1000 (maximum value is ~700):
 k_matrix(k_matrix>1000) = 1000;
 
-%Filter out values that are inconsistent at different alphas:
-always_down  = sum(k_matrix <= (1-tol),2) == length(alpha);
-always_up    = sum(k_matrix >= (1+tol),2) == length(alpha);
-%Identify those reactions with mixed patterns 
+% Filter out values that are inconsistent at different alphas:
+always_down  = sum(k_matrix <= (1-tolerance),2) == length(alpha);
+always_up    = sum(k_matrix >= (1+tolerance),2) == length(alpha);
+
+% Identify those reactions with mixed patterns
 incons_rxns  = always_down + always_up == 0;
-%Identify genes that are linked to "inconsistent rxns"
+
+% Identify genes that are linked to "inconsistent rxns"
 incons_genes = sum(rxnGeneM(incons_rxns,:),1) > 0;
-%Finally, inconsistent reactions are those that are not conected
-%to "inconsistent genes" from the original "inconsistent rxns" set
+
+% Finally, inconsistent reactions are those that are not conected
+% to "inconsistent genes" from the original "inconsistent rxns" set
 incons_rxns  = sum(rxnGeneM(:,incons_genes),2) > 0;
-%Keep results for the consistent rxns exclusively
+
+% Keep results for the consistent rxns exclusively
 v_matrix     = v_matrix(~incons_rxns,:);
 k_matrix     = k_matrix(~incons_rxns,:);
 rxnGeneM     = rxnGeneM(~incons_rxns,:);
 FC.rxns      = FC.rxns(~incons_rxns,:);
-%Get median k score across steps
+
+% Get median k score across steps
 FC.k_rxns   = mean(k_matrix,2);
-%Remove arm_rxns from the list of rxns with K_score>1
-armUP       = startsWith(FC.rxns(:,1),'arm_') & FC.k_rxns>1;
-FC.k_rxns   = FC.k_rxns(~armUP,:);
-FC.v_matrix = v_matrix(~armUP,:);
-FC.k_matrix = k_matrix(~armUP,:);
-rxnGeneM    = rxnGeneM(~armUP,:);
-FC.rxns     = FC.rxns(~armUP,:);
-
-%Order from highest to lowest median k_score (across alphas)
+
+% Order from highest to lowest median k_score (across alphas)
 [~,order]   = sort(FC.k_rxns,'descend');
 FC.k_rxns   = FC.k_rxns(order,:);
 FC.v_matrix = v_matrix(order,:);
 FC.k_matrix = k_matrix(order,:);
 rxnGeneM    = rxnGeneM(order,:);
 FC.rxns     = FC.rxns(order,:);
-%Create list of remaining genes and filter out any inconsistent score:
-%Just those genes that are connected to the remaining rxns are
-FC.genes     = model.genes(sum(rxnGeneM,1) > 0);
-FC.geneNames = model.geneShortNames(sum(rxnGeneM,1) > 0);
+
+% Create list of remaining genes and filter out any inconsistent score:
+% Just those genes that are connected to the remaining rxns are
+FC.genes     = ecModel.genes(sum(rxnGeneM,1) > 0);
+FC.geneNames = ecModel.geneShortNames(sum(rxnGeneM,1) > 0);
 FC.k_genes   = zeros(size(FC.genes));
 cons_genes   = false(size(FC.genes));
 rxnGeneM     = rxnGeneM(:,sum(rxnGeneM,1) > 0);
 
 for i = 1:length(FC.genes)
-	%Extract all the K_scores (from rxns across alphas) conected to
-	%each remaining gene
+	% Extract all the K_scores (from rxns across alphas) conected to
+	% each remaining gene
     k_set         = FC.k_rxns(rxnGeneM(:,i) > 0);
-    %Check the kind of control that gene i-th exerts over its reactions
-    always_down   = sum(k_set <= (1-tol)) == length(k_set);
-    always_up     = sum(k_set >= (1+tol)) == length(k_set);
-    %Evaluate if gene is always exerting either a positive or negative
-    %control
+    % Check the kind of control that gene i-th exerts over its reactions
+    always_down   = sum(k_set <= (1-tolerance)) == length(k_set);
+    always_up     = sum(k_set >= (1+tolerance)) == length(k_set);
+    % Evaluate if gene is always exerting either a positive or negative
+    % control
     cons_genes(i) = always_down + always_up == 1;
     FC.k_genes(i) = mean(k_set);
 end
-%Keep "consistent genes"
+% Keep "consistent genes"
 FC.genes     = FC.genes(cons_genes);
 FC.geneNames = FC.geneNames(cons_genes);
 FC.k_genes   = FC.k_genes(cons_genes);
 rxnGeneM     = rxnGeneM(:,cons_genes);
 
 if filterG
-	%Filter any value between mean(alpha) and 1:
-	unchanged    = (FC.k_genes >= mean(alpha) - tol) + (FC.k_genes <= 1 + tol) == 2;
+	% Filter any value between mean(alpha) and 1:
+	unchanged    = (FC.k_genes >= mean(alpha) - tolerance) + (FC.k_genes <= 1 + tolerance) == 2;
 	FC.genes     = FC.genes(~unchanged);
 	FC.geneNames = FC.geneNames(~unchanged);
 	FC.k_genes   = FC.k_genes(~unchanged);
     rxnGeneM     = rxnGeneM(:,~unchanged);
-    %Update results for rxns-related fields (remove remaining reactions
-    %without any associated gene in rxnGeneM)
+    % Update results for rxns-related fields (remove remaining reactions
+    % without any associated gene in rxnGeneM)
     toKeep      = (sum(rxnGeneM,2) > 0);
     FC.k_rxns   = FC.k_rxns(toKeep,:);
     FC.v_matrix = v_matrix(toKeep,:);
@@ -132,7 +137,7 @@
     FC.rxns     = FC.rxns(toKeep,:);
 end
 
-%Order from highest to lowest k:
+% Order from highest to lowest k:
 [~,order]    = sort(FC.k_genes,'descend');
 FC.genes     = FC.genes(order,:);
 FC.geneNames = FC.geneNames(order,:);

From d3e233e16d1bc45526e4f69927afce5e9596b112 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Tue, 4 Apr 2023 17:45:49 +0200
Subject: [PATCH 05/13] refactor: add compatibility to GECKO3

---
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m | 122 +++++++++++--------
 1 file changed, 71 insertions(+), 51 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index e54464131..fb7b33c9f 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -1,53 +1,72 @@
-function results = run_ecFSEOF(model,rxnTarget,cSource,alphaLims,Nsteps,file1,file2)
-% 
+function results = run_ecFSEOF(ecModel,rxnTarget,cSource,alphaLims,Nsteps,file_genes,file_rxns,modelAdapter)
 % run_ecFSEOF
 %
 % 	Function that runs Flux-scanning with Enforced Objective Function
 % 	for a specified production target.
-%   
-%		model     (struct) ecModel with total protein pool constraint.
-%		rxnTarget (string) Rxn ID for the production target reaction, 
-% 				  a exchange reaction is recommended.
-%		cSource	  (string) Rxn ID for the main carbon source uptake 
-%			      reaction (make sure that the correct directionality is
-%			      indicated).
-% 		alphaLims (vector) Minimum and maximum biomass yield [gDw/mmol Csource] 
-%				  for enforced objective limits
-%       Nsteps    (integer) Number of steps for suboptimal objective 
-%                 in FSEOF
-%		file1 	  (string) opt, file name for an output .txt tab-separated 
-% 				  file for the results at the genes level
-%		file2 	  (string) opt, file name for an output .txt tab-separated 
-% 				  file for the results at the reactions level		
 %
-% Usage: run_ecFSEOF(ecModel,rxnTarget,cSource,alphaLims,Nsteps)
+% Input:
+%   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
+%	rxnTarget       rxn ID for the production target reaction, a exchange
+%                   reaction is recommended.
+%	cSource	        rxn ID for the main carbon source uptake reaction.
+% 	alphaLims       vector of Minimum and maximum biomass yield [gDW/mmol Csource]
+%			        for enforced objective limits (e.g. [0.5 2])
+%   Nsteps          number of steps for suboptimal objective in FSEOF.
+%                   (Optional, default 16)
+%	file_genes 	    file name for results output at the genes level.
+%                   (Optional, default output in the model-specific 'output'
+%                   sub-folder taken from modelAdapter,
+%                   e.g. GECKO/userData/ecYeastGEM/output/ecFSEOF_genes.tsv)
+%	file_rxns 	    file name for results output at the reactions level.
+%                   (Optional, default output in the model-specific 'output'
+%                   sub-folder taken from modelAdapter,
+%                   e.g. GECKO/userData/ecYeastGEM/output/ecFSEOF_rxns.tsv)
+%   modelAdapter    a loaded model adapter. (Optional, will otherwise use
+%                   the default model adapter)
 %
+% Usage:
+%   results = run_ecFSEOF(ecModel,rxnTarget,cSource,alphaLims)
+%
+
+if nargin < 8 || isempty(modelAdapter)
+    modelAdapter = ModelAdapterManager.getDefaultAdapter();
+    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 < 7 || isempty(file_rxns)
+	file_rxns = fullfile(params.path,'output','ecFSEOF_rxns.tsv');;
+end
 
-if nargin < 7
-	file2 = [];
-	if nargin < 6
-		file1 = [];
-	end
+if nargin < 6 || isempty(file_genes)
+    file_genes = fullfile(params.path,'output','ecFSEOF_genes.tsv');;
 end
-%Check model format, if COBRA model then convert to RAVEN format
-if isfield('model','rules') && ~isfield('model','metComps')
-    model = ravenCobraWrapper(model);
+
+if nargin < 5 || isempty(Nsteps)
+    Nsteps = 16;
 end
-%Define alpha vector for suboptimal enforced objective values
+
+% Define alpha vector for suboptimal enforced objective values
 alphaV  = alphaLims(1):((alphaLims(2)-alphaLims(1))/(Nsteps-1)):alphaLims(2);
-%Standardize grRules and rxnGeneMat in model
-[grRules,rxnGeneMat] = standardizeGrRules(model,true);
-model.grRules        = grRules;
-model.rxnGeneMat     = rxnGeneMat;
-%run FSEOF analysis
-results = ecFlux_scanning(model,rxnTarget,cSource,alphaV,1E-4,true);
-%Create gene table:
+
+% Standardize grRules and rxnGeneMat in model
+[grRules,rxnGeneMat] = standardizeGrRules(ecModel,true);
+ecModel.grRules      = grRules;
+ecModel.rxnGeneMat   = rxnGeneMat;
+
+% run FSEOF analysis
+results = ecFlux_scanning(ecModel,rxnTarget,cSource,alphaV,1E-4,true);
+
+% Create gene table:
 results.geneTable      = cell(length(results.genes),3);
 results.geneTable(:,1) = results.genes;
 results.geneTable(:,2) = results.geneNames;
 results.geneTable(:,3) = num2cell(results.k_genes);
-%Create rxns table (exclude enzyme usage reactions):
-toKeep                 = find(~startsWith(results.rxns(:,1),'draw_prot_'));
+
+% Create rxns table (exclude enzyme usage reactions):
+toKeep                 = find(~startsWith(results.rxns(:,1),'usage_prot_'));
 results.k_rxns         = results.k_rxns(toKeep);
 results.k_matrix       = results.k_matrix(toKeep,:);
 results.v_matrix       = results.v_matrix(toKeep,:);
@@ -57,21 +76,22 @@
 results.rxnsTable(:,3) = num2cell(results.k_rxns);
 results.rxnsTable(:,4) = results.rxns(toKeep,3);
 results.rxnsTable(:,5) = results.rxns(toKeep,4);
-try
-    if ~isempty(file1)
-        varNames = {'gene_IDs' 'gene_names' 'K_score'}; 
-    	T        = cell2table(results.geneTable,'VariableNames',varNames);
-    	writetable(T,file1,'Delimiter','\t','QuoteStrings',false) 
-    end
-    if ~isempty(file2)
-        varNames = {'rxn_IDs' 'rxnNames' 'K_score' 'grRules' 'rxn_formula'};
-    	T        = cell2table(results.rxnsTable,'VariableNames',varNames);
-        writetable(T,file2,'Delimiter','\t','QuoteStrings',false) 
-    end
-catch
-    warning('Output files directory not found');
-end
-%Remove redundant output fields
+
+writetable(cell2table(results.geneTable, ...
+    'VariableNames', {'gene_IDs' 'gene_names' 'K_score'}), ...
+    file_genes, ...
+    'Delimiter', '\t', ...
+    'QuoteStrings', false);
+
+writetable(cell2table(results.rxnsTable, ...
+    'VariableNames', {'rxn_IDs' 'rxnNames' 'K_score' 'grRules' 'rxn_formula'}), ...
+    file_rxns, ...
+    'Delimiter', '\t', ...
+    'QuoteStrings', false);
+
+disp('ecFSEOF results stored at %s\n',fileparts(file_genes));
+
+% Remove redundant output fields
 results = rmfield(results,'k_rxns');
 results = rmfield(results,'rxns');
 results = rmfield(results,'genes');

From 04e672435a4b22395397e3c3d9c0cda5fa74d0ad Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Wed, 5 Apr 2023 09:24:08 +0200
Subject: [PATCH 06/13] fix: typo

---
 src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m | 2 +-
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
index 8ff2a81e2..ebfdffabb 100644
--- a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
+++ b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
@@ -17,7 +17,7 @@
 % Usage:
 %   FC = ecFlux_scanning(model,target,cSource,alpha,tolerance,filterG)
 
-if nargin < 6 || isempty(filerG)
+if nargin < 6 || isempty(filterG)
 	filterG = false;
 end
 
diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index fb7b33c9f..314459cc4 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -10,7 +10,7 @@
 %                   reaction is recommended.
 %	cSource	        rxn ID for the main carbon source uptake reaction.
 % 	alphaLims       vector of Minimum and maximum biomass yield [gDW/mmol Csource]
-%			        for enforced objective limits (e.g. [0.5 2])
+%			        for enforced objective limits (e.g. [0.5 1]). Max value: 1.
 %   Nsteps          number of steps for suboptimal objective in FSEOF.
 %                   (Optional, default 16)
 %	file_genes 	    file name for results output at the genes level.

From bc86f7122e77fbf2d3ee18376b76872cf953fb00 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Wed, 5 Apr 2023 10:39:52 +0200
Subject: [PATCH 07/13] fix: save file type file

---
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index 314459cc4..f1217dd8a 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -37,11 +37,11 @@
 params = modelAdapter.getParameters();
 
 if nargin < 7 || isempty(file_rxns)
-	file_rxns = fullfile(params.path,'output','ecFSEOF_rxns.tsv');;
+	file_rxns = fullfile(params.path,'output','ecFSEOF_rxns.tsv');
 end
 
 if nargin < 6 || isempty(file_genes)
-    file_genes = fullfile(params.path,'output','ecFSEOF_genes.tsv');;
+    file_genes = fullfile(params.path,'output','ecFSEOF_genes.tsv');
 end
 
 if nargin < 5 || isempty(Nsteps)
@@ -80,16 +80,18 @@
 writetable(cell2table(results.geneTable, ...
     'VariableNames', {'gene_IDs' 'gene_names' 'K_score'}), ...
     file_genes, ...
+    'FileType', 'text', ...
     'Delimiter', '\t', ...
     'QuoteStrings', false);
 
 writetable(cell2table(results.rxnsTable, ...
     'VariableNames', {'rxn_IDs' 'rxnNames' 'K_score' 'grRules' 'rxn_formula'}), ...
     file_rxns, ...
+    'FileType', 'text', ...
     'Delimiter', '\t', ...
     'QuoteStrings', false);
 
-disp('ecFSEOF results stored at %s\n',fileparts(file_genes));
+disp(['ecFSEOF results stored at: ' newline fileparts(file_genes)]);
 
 % Remove redundant output fields
 results = rmfield(results,'k_rxns');

From 88512798dd3b5dcf9d7b5125d3e973993061f7b9 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Wed, 5 Apr 2023 11:24:32 +0200
Subject: [PATCH 08/13] fix: filter tolerance

---
 src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
index ebfdffabb..c97aa9477 100644
--- a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
+++ b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
@@ -65,8 +65,8 @@
 k_matrix(k_matrix>1000) = 1000;
 
 % Filter out values that are inconsistent at different alphas:
-always_down  = sum(k_matrix <= (1-tolerance),2) == length(alpha);
-always_up    = sum(k_matrix >= (1+tolerance),2) == length(alpha);
+always_down  = sum(k_matrix <= 1,2) == length(alpha);
+always_up    = sum(k_matrix >= 1,2) == length(alpha);
 
 % Identify those reactions with mixed patterns
 incons_rxns  = always_down + always_up == 0;

From d2f2fcaf08de7c3eac34ab9b40e8f801caaa95e5 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Wed, 5 Apr 2023 11:25:25 +0200
Subject: [PATCH 09/13] fix: if not carbon uptake have been defined

---
 src/geckomat/utilities/getFluxTarget.m | 25 ++++++++++++++++++++-----
 1 file changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/geckomat/utilities/getFluxTarget.m b/src/geckomat/utilities/getFluxTarget.m
index 0ab939975..64d9e3bf1 100644
--- a/src/geckomat/utilities/getFluxTarget.m
+++ b/src/geckomat/utilities/getFluxTarget.m
@@ -47,11 +47,26 @@
 % Fix carbon source uptake
 cSourceIdx = strcmpi(ecModel.rxns,cSource);
 uptake = ecModel.lb(cSourceIdx);
-ecModel = setParam(ecModel, 'var', cSource, uptake, tolerance);
 
-% Max growth and fix growth suboptimal
-ecModel = setParam(ecModel, 'obj', params.bioRxn, 1);
-sol = solveLP(ecModel);
+% Check if a carbon source uptake rate have been defined
+if uptake ~= -1000
+    ecModel = setParam(ecModel, 'var', cSource, uptake, tolerance);
+    
+    % Maximize growth
+    ecModel = setParam(ecModel, 'obj', params.bioRxn, 1);
+    sol = solveLP(ecModel);
+else
+    % Maximize growth an use the uptake rate predicted to constraint carbon source
+    ecModel = setParam(ecModel, 'obj', params.bioRxn, 1);
+    sol = solveLP(ecModel, 1);
+
+    ecModel = setParam(ecModel, 'var', cSource, sol.x(cSourceIdx), tolerance);
+    
+    warning(['The maximum uptake rate for the carbon source is -1000. ' ...
+        num2str(abs(sol.x(cSourceIdx))) ' [mmol/gDW/h] was used instead'])
+end
+
+% Fix growth suboptimal
 bioRxnIdx = strcmpi(ecModel.rxns, params.bioRxn);
 ecModel = setParam(ecModel, 'lb', params.bioRxn, sol.x(bioRxnIdx) * (1-tolerance) * alpha);
 
@@ -59,7 +74,7 @@
 poolIdx = strcmpi(ecModel.rxns, 'prot_pool_exchange');
 ecModel = setParam(ecModel, 'obj', 'prot_pool_exchange', 1);
 sol = solveLP(ecModel, 1);
-ecModel = setParam(ecModel, 'lb', params.bioRxn, sol.x(poolIdx));
+ecModel = setParam(ecModel, 'lb', 'prot_pool_exchange', sol.x(poolIdx));
 
 % Minimize target
 ecModel = setParam(ecModel, 'obj', target, -1);

From dd5d295d988605018070b40f63cc8815fbf8ac01 Mon Sep 17 00:00:00 2001
From: Albert Tafur Rangel <39610893+ae-tafur@users.noreply.github.com>
Date: Wed, 5 Apr 2023 14:14:50 +0200
Subject: [PATCH 10/13] fix: add toleracne to protein pool

---
 src/geckomat/utilities/getFluxTarget.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geckomat/utilities/getFluxTarget.m b/src/geckomat/utilities/getFluxTarget.m
index 64d9e3bf1..bd5b20609 100644
--- a/src/geckomat/utilities/getFluxTarget.m
+++ b/src/geckomat/utilities/getFluxTarget.m
@@ -73,8 +73,8 @@
 % Minimize prot_pool_exchange and fix
 poolIdx = strcmpi(ecModel.rxns, 'prot_pool_exchange');
 ecModel = setParam(ecModel, 'obj', 'prot_pool_exchange', 1);
-sol = solveLP(ecModel, 1);
-ecModel = setParam(ecModel, 'lb', 'prot_pool_exchange', sol.x(poolIdx));
+sol = solveLP(ecModel);
+ecModel = setParam(ecModel, 'lb', 'prot_pool_exchange', sol.x(poolIdx) * (1+tolerance));
 
 % Minimize target
 ecModel = setParam(ecModel, 'obj', target, -1);

From 7c605e0d0b012b580b5e941a17ecf6d11879c217 Mon Sep 17 00:00:00 2001
From: edkerk <eduardk@chalmers.se>
Date: Tue, 18 Apr 2023 10:31:45 +0200
Subject: [PATCH 11/13] refactor: spaces instead of tabs

---
 .../utilities/ecFSEOF/ecFlux_scanning.m       | 28 +++++++++----------
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m  | 22 +++++++--------
 src/geckomat/utilities/getFluxTarget.m        |  8 +++---
 3 files changed, 29 insertions(+), 29 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
index c97aa9477..ea508d161 100644
--- a/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
+++ b/src/geckomat/utilities/ecFSEOF/ecFlux_scanning.m
@@ -1,24 +1,24 @@
 function FC = ecFlux_scanning(ecModel,target,cSource,alpha,tolerance,filterG)
-% ecFlux_scanning
+%ecFlux_scanning
 %
 % Input:
 %   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
-%	rxnTarget       rxn ID for the production target reaction, a exchange
+%   rxnTarget       rxn ID for the production target reaction, a exchange
 %                   reaction is recommended.
-%	cSource	        rxn ID for the main carbon source uptake reaction.
+%  cSource          rxn ID for the main carbon source uptake reaction.
 %   alpha           scalling factor for production yield for enforced objective
 %                   limits
 %   tolerance       numerical tolerance for fixing bounds.
 %                   (Optional, defaul 1E-4)
-%   filterG	        logical value. TRUE if genes K_scores results should be
-%		            filtered according to the alpha vector distribution
+%   filterG         logical value. TRUE if genes K_scores results should be
+%                   filtered according to the alpha vector distribution
 %                   (Optional, defaul false)
 %
 % Usage:
 %   FC = ecFlux_scanning(model,target,cSource,alpha,tolerance,filterG)
 
 if nargin < 6 || isempty(filterG)
-	filterG = false;
+    filterG = false;
 end
 
 if nargin < 5 || isempty(tolerance)
@@ -104,8 +104,8 @@
 rxnGeneM     = rxnGeneM(:,sum(rxnGeneM,1) > 0);
 
 for i = 1:length(FC.genes)
-	% Extract all the K_scores (from rxns across alphas) conected to
-	% each remaining gene
+    % Extract all the K_scores (from rxns across alphas) conected to
+    % each remaining gene
     k_set         = FC.k_rxns(rxnGeneM(:,i) > 0);
     % Check the kind of control that gene i-th exerts over its reactions
     always_down   = sum(k_set <= (1-tolerance)) == length(k_set);
@@ -122,11 +122,11 @@
 rxnGeneM     = rxnGeneM(:,cons_genes);
 
 if filterG
-	% Filter any value between mean(alpha) and 1:
-	unchanged    = (FC.k_genes >= mean(alpha) - tolerance) + (FC.k_genes <= 1 + tolerance) == 2;
-	FC.genes     = FC.genes(~unchanged);
-	FC.geneNames = FC.geneNames(~unchanged);
-	FC.k_genes   = FC.k_genes(~unchanged);
+    % Filter any value between mean(alpha) and 1:
+    unchanged    = (FC.k_genes >= mean(alpha) - tolerance) + (FC.k_genes <= 1 + tolerance) == 2;
+    FC.genes     = FC.genes(~unchanged);
+    FC.geneNames = FC.geneNames(~unchanged);
+    FC.k_genes   = FC.k_genes(~unchanged);
     rxnGeneM     = rxnGeneM(:,~unchanged);
     % Update results for rxns-related fields (remove remaining reactions
     % without any associated gene in rxnGeneM)
@@ -143,4 +143,4 @@
 FC.geneNames = FC.geneNames(order,:);
 FC.k_genes   = FC.k_genes(order,:);
 
-end
\ No newline at end of file
+end
diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index f1217dd8a..5cc679806 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -1,23 +1,23 @@
 function results = run_ecFSEOF(ecModel,rxnTarget,cSource,alphaLims,Nsteps,file_genes,file_rxns,modelAdapter)
-% run_ecFSEOF
+%run_ecFSEOF
 %
-% 	Function that runs Flux-scanning with Enforced Objective Function
-% 	for a specified production target.
+% Function that runs Flux-scanning with Enforced Objective Function
+% for a specified production target.
 %
 % Input:
 %   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
-%	rxnTarget       rxn ID for the production target reaction, a exchange
+%   rxnTarget       rxn ID for the production target reaction, a exchange
 %                   reaction is recommended.
-%	cSource	        rxn ID for the main carbon source uptake reaction.
-% 	alphaLims       vector of Minimum and maximum biomass yield [gDW/mmol Csource]
-%			        for enforced objective limits (e.g. [0.5 1]). Max value: 1.
+%   cSource         rxn ID for the main carbon source uptake reaction.
+%   alphaLims       vector of Minimum and maximum biomass yield [gDW/mmol Csource]
+%                   for enforced objective limits (e.g. [0.5 1]). Max value: 1.
 %   Nsteps          number of steps for suboptimal objective in FSEOF.
 %                   (Optional, default 16)
-%	file_genes 	    file name for results output at the genes level.
+%   file_genes      file name for results output at the genes level.
 %                   (Optional, default output in the model-specific 'output'
 %                   sub-folder taken from modelAdapter,
 %                   e.g. GECKO/userData/ecYeastGEM/output/ecFSEOF_genes.tsv)
-%	file_rxns 	    file name for results output at the reactions level.
+%   file_rxns       file name for results output at the reactions level.
 %                   (Optional, default output in the model-specific 'output'
 %                   sub-folder taken from modelAdapter,
 %                   e.g. GECKO/userData/ecYeastGEM/output/ecFSEOF_rxns.tsv)
@@ -37,7 +37,7 @@
 params = modelAdapter.getParameters();
 
 if nargin < 7 || isempty(file_rxns)
-	file_rxns = fullfile(params.path,'output','ecFSEOF_rxns.tsv');
+    file_rxns = fullfile(params.path,'output','ecFSEOF_rxns.tsv');
 end
 
 if nargin < 6 || isempty(file_genes)
@@ -99,4 +99,4 @@
 results = rmfield(results,'genes');
 results = rmfield(results,'geneNames');
 results = rmfield(results,'k_genes');
-end
\ No newline at end of file
+end
diff --git a/src/geckomat/utilities/getFluxTarget.m b/src/geckomat/utilities/getFluxTarget.m
index bd5b20609..67a58fd52 100644
--- a/src/geckomat/utilities/getFluxTarget.m
+++ b/src/geckomat/utilities/getFluxTarget.m
@@ -1,5 +1,5 @@
 function [minFlux, maxFlux] = getFluxTarget(ecModel,target,cSource,alpha,tolerance,modelAdapter)
-% getFluxTarget
+%getFluxTarget
 %
 %   Function that performs a series of LP optimizations on an ecModel,
 %   by first maximizing biomass, then fixing a suboptimal value and 
@@ -9,9 +9,9 @@
 %
 % Input:
 %   ecModel         an ecModel in GECKO 3 format (with ecModel.ec structure).
-%	target          rxn ID for the production target reaction, a exchange
+%   target          rxn ID for the production target reaction, a exchange
 %                   reaction is recommended.
-%	cSource	        rxn ID for the main carbon source uptake reaction.
+%   cSource         rxn ID for the main carbon source uptake reaction.
 %   alpha           scalling factor for desired suboptimal growth.
 %                   (Optional, defaul 1)
 %   tolerance       numerical tolerance for fixing bounds
@@ -86,4 +86,4 @@
 sol = solveLP(ecModel);
 maxFlux = sol.x;
 
-end
\ No newline at end of file
+end

From d382f1d1d40efbb3af1f917c26dbcd65fa0683c8 Mon Sep 17 00:00:00 2001
From: ae-tafur <ae.tafur@hotmail.com>
Date: Mon, 24 Apr 2023 09:25:41 +0200
Subject: [PATCH 12/13] fix: function input description

---
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index 5cc679806..4672995ed 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -9,8 +9,8 @@
 %   rxnTarget       rxn ID for the production target reaction, a exchange
 %                   reaction is recommended.
 %   cSource         rxn ID for the main carbon source uptake reaction.
-%   alphaLims       vector of Minimum and maximum biomass yield [gDW/mmol Csource]
-%                   for enforced objective limits (e.g. [0.5 1]). Max value: 1.
+% 	alphaLims       vector of Minimum and maximum biomass scalling factors for
+%			        enforced objective limits (e.g. [0.5 1]). Max value: 1.
 %   Nsteps          number of steps for suboptimal objective in FSEOF.
 %                   (Optional, default 16)
 %   file_genes      file name for results output at the genes level.

From 35c50ba4c45c0075538f2f27d8c3d914c4926e2f Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Fri, 12 May 2023 01:31:35 +0200
Subject: [PATCH 13/13] Update src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m

---
 src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
index 4672995ed..d64a1557a 100644
--- a/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
+++ b/src/geckomat/utilities/ecFSEOF/run_ecFSEOF.m
@@ -9,8 +9,8 @@
 %   rxnTarget       rxn ID for the production target reaction, a exchange
 %                   reaction is recommended.
 %   cSource         rxn ID for the main carbon source uptake reaction.
-% 	alphaLims       vector of Minimum and maximum biomass scalling factors for
-%			        enforced objective limits (e.g. [0.5 1]). Max value: 1.
+%   alphaLims       vector of Minimum and maximum biomass scalling factors for
+%                   enforced objective limits (e.g. [0.5 1]). Max value: 1.
 %   Nsteps          number of steps for suboptimal objective in FSEOF.
 %                   (Optional, default 16)
 %   file_genes      file name for results output at the genes level.