cycle free FBA implemented

xbanalysis/problems/fba_problem.py

 import numpy as np
+import math
 from scipy import sparse
 
 from xbanalysis.solvers.glpk_linear_problem import GlpkLinearProblem
@@ -192,6 +193,45 @@ class FbaProblem:
         """
         self.lp.remove_obj_coefs(self.get_coefficients(objective))
 
+    def update_lp_flux_bounds(self, flux_bounds):
+        """update flux bounds for selected reactions in the linear problem.
+
+        lp problem might need to be instantiated
+        overwrite bounds of protected functions will be cleaned from already removed reactions
+
+        :param flux_bounds: lower and upper flux bounds for selected reactions
+        :type flux_bounds: dict: key: reaction id, value: list-like of two float values
+        """
+        if self.lp is None:
+            self.lp = self.create_core_fba_lp()
+        for rid, bound_upd in flux_bounds.items():
+            if rid in self.rids:
+                temp_bounds = self.flux_bounds[:, self.rid2idx[rid]].copy()
+                for i in [0, 1]:
+                    if math.isfinite(bound_upd[i]):
+                        temp_bounds[i] = bound_upd[i]
+                self.lp.set_col_bnd(self.rid2idx[rid], np.fmax(temp_bounds, 0))
+                if rid in self.rids_rev:
+                    self.lp.set_col_bnd(self.rid2idx_rev[rid], np.fmax(np.flip(-temp_bounds), 0))
+
+    def reset_lp_flux_bounds(self, rids):
+        """restore original flux bounds for selected reactions in linear problem.
+
+        overwrite bounds of protected functions will be cleaned from already removed reactions
+
+        :param rids: reaction ids for which to restore the model orignal flux bounds
+        :type rids: list-like of str
+        """
+        assert self.lp is not None
+
+        for rid in rids:
+            if rid in self.rids:
+                assert rid in self.rids
+                bounds = self.flux_bounds[:, self.rid2idx[rid]]
+                self.lp.set_col_bnd(self.rid2idx[rid], np.fmax(bounds, 0))
+                if rid in self.rids_rev:
+                    self.lp.set_col_bnd(self.rid2idx_rev[rid], np.fmax(np.flip(-bounds), 0))
+
     def fba_optimize(self, short_result=False):
         """FBA optimization of current problem
 
@@ -252,6 +292,62 @@ class FbaProblem:
         res['reduced_costs'] = self.combine_fwd_rev(res['reduced_costs'])
         return res
 
+    def cyclefree_fba_optimize(self):
+        """cycle free FBA (remove loops from FBA solution)
+
+        Reference:
+        Desouki et al., 2015, CycleFreeFlux: efficient removal of thermodynamically
+        infeasible loops from flux distributions
+
+        Method:
+        - run FBA and determine a flux distribution
+        - implement model objective function as equality constraint with FBA optimum
+        - configure a new objective function ('minimize' absolut fluxes)
+            - exchange reactions are not considered and fluxes are fixed at FBA level
+            - reactions with positive FBA fluxes get a coefficient of 1.0 and fluxes bound between zero and FBA level
+            - reactions with negative FBA fluxes get a coefficient of -1.0 and fluxes bound between FBA level and zero
+                - i.e. internal reactions are only allowed to shrink, but not change direction
+            - run FBA on new objective and retrieve new flux values for reactions.
+                - use new flux values with model objective value as optimum
+        :return: result of cycleFree FBA solution
+        :rtype: dict
+        """
+        # get optimum value for model objective function
+        res = self.fba_optimize(short_result=False)
+        if res['success'] is False:
+            return res
+        optimum = res['fun']
+
+        cff_objective = {}
+        cff_flux_bounds = {}
+        for rid, idx in self.rid2idx.items():
+            flux = res['x'][idx]
+            if len(self.xba_model.reactions[rid].products) == 0:
+                cff_flux_bounds[rid] = [flux, flux]
+            elif flux >= 0:
+                cff_flux_bounds[rid] = [0.0, flux]
+                cff_objective[rid] = 1.0
+            else:
+                cff_flux_bounds[rid] = [flux, 0.0]
+                cff_objective[rid] = -1.0
+
+        # implement model main objective as an equality constraint with optimal value
+        row_bnds = np.array([optimum, optimum])
+        row_idx0 = self.lp.add_constraint(self.get_coefficients(self.objective), row_bnds)
+
+        # configure cyclefree FBA objective
+        self.set_lp_objective(cff_objective, 'minimize')
+        self.update_lp_flux_bounds(cff_flux_bounds)
+        res = self.lp.solve(short_result=False)
+        self.reset_lp_flux_bounds(set(cff_flux_bounds))
+        self.remove_lp_objective(cff_objective)
+        self.lp.del_constraint(row_idx0)
+
+        res['fun'] = optimum
+        res['x'] = self.combine_fwd_rev(res['x'])
+        res['reduced_costs'] = self.combine_fwd_rev(res['reduced_costs'])
+        return res
+
     def fva_optimize(self, rids=None, frac_of_opt=1.0):
         """Flux variability analyzis on reactions within fraction of optimum