From 3d59a4a920bd22a5183e29cee8bdafd7660e3744 Mon Sep 17 00:00:00 2001 From: Jon Lund Steffensen Date: Tue, 28 Feb 2017 11:26:20 -0500 Subject: [PATCH] gapfind: Alternative GapFind algorithm which is more robust Does not depend on v_max, instead sets a maximum flux when solving this problem to 1. This reduces the distance between v_max and epsilon making the problem more robust. This is able to solve larger models that fail with the previous implementation. This respects reaction directionality but ignores any reaction flux limits. This also lowers the maximum flux which scales the solution. This means that some reactions may not be able to produce epsilon flux. --- psamm/commands/gapcheck.py | 3 +-- psamm/gapfill.py | 15 ++++++++++++--- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/psamm/commands/gapcheck.py b/psamm/commands/gapcheck.py index fa4f8cbf..e10373c6 100644 --- a/psamm/commands/gapcheck.py +++ b/psamm/commands/gapcheck.py @@ -46,11 +46,10 @@ def run(self): solver = self._get_solver(integer=True) extracellular_comp = self._model.extracellular_compartment epsilon = self._args.epsilon - v_max = float(self._model.default_flux_limit) # Run GapFind on model logger.info('Searching for blocked compounds...') - result = gapfind(self._mm, solver=solver, epsilon=epsilon, v_max=v_max) + result = gapfind(self._mm, solver=solver, epsilon=epsilon) try: blocked = set(compound for compound in result diff --git a/psamm/gapfill.py b/psamm/gapfill.py index 8575353d..b59903f8 100644 --- a/psamm/gapfill.py +++ b/psamm/gapfill.py @@ -49,7 +49,7 @@ def _find_integer_tolerance(epsilon, v_max, min_tol): return int_tol -def gapfind(model, solver, epsilon=0.001, v_max=1000, implicit_sinks=True): +def gapfind(model, solver, epsilon=0.001, implicit_sinks=True): """Identify compounds in the model that cannot be produced. Yields all compounds that cannot be produced. This method @@ -58,7 +58,9 @@ def gapfind(model, solver, epsilon=0.001, v_max=1000, implicit_sinks=True): produced is the presence of the compounds needed to produce it. Epsilon indicates the threshold amount of reaction flux for the products - to be considered non-blocked. V_max indicates the maximum flux. + to be considered non-blocked. The method respects the reversibility of + reactions but does not consider flux bounds of reactions i.e. every + reaction is assumed to be able to run at maximum flux. This method is implemented as a MILP-program. Therefore it may not be efficient for larger models. @@ -75,6 +77,8 @@ def gapfind(model, solver, epsilon=0.001, v_max=1000, implicit_sinks=True): """ prob = solver.create_problem() + v_max = 1 + # Set integrality tolerance such that w constraints are correct min_tol = prob.integrality_tolerance.min int_tol = _find_integer_tolerance(epsilon, v_max, min_tol) @@ -96,7 +100,12 @@ def gapfind(model, solver, epsilon=0.001, v_max=1000, implicit_sinks=True): w.define([spec]) w_var = w(spec) - lower, upper = (float(x) for x in model.limits[reaction_id]) + lower, upper = 0, 0 + if model.limits[reaction_id].lower < 0: + lower = -v_max + if model.limits[reaction_id].upper > 0: + upper = v_max + if value > 0: dv = v(reaction_id) else: