fix: production envelope with unspecified c-source

Fix bug in production_envelopee occurring for models where the C-source
was formulated as -> x instead of x <-

Also support guessing the C-source reaction.

Author: hredestig

cobra/flux_analysis/phenotype_phase_plane.py

@@ -1,10 +1,12 @@
 # -*- coding: utf-8 -*-
 
-from __future__ import absolute_import
+from __future__ import absolute_import, division
 
 from warnings import warn
 from collections import defaultdict
+from operator import itemgetter
 from itertools import product
+from six import iteritems
 from multiprocessing import Pool
 import pandas as pd
 from optlang.interface import OPTIMAL
@@ -16,6 +18,7 @@
 from cobra.solvers import get_solver_name, solver_dict
 import cobra.util.solver as sutil
 from cobra.flux_analysis import flux_variability_analysis as fva
+from cobra.exceptions import OptimizationError
 
 # attempt to import plotting libraries
 try:
@@ -351,7 +354,7 @@ def production_envelope(model, reactions, objective=None, c_source=None,
         - carbon_yield: if carbon source is defined and the product is a
           single metabolite (mol carbon product per mol carbon feeding source)
 
-        - _mass_yield: if carbon source is defined and the product is a
+        - mass_yield: if carbon source is defined and the product is a
           single metabolite (gram product per 1 g of feeding source)
 
         - direction: the direction of the optimization.
@@ -371,36 +374,28 @@ def production_envelope(model, reactions, objective=None, c_source=None,
                          'based solver interfaces.')
     reactions = model.reactions.get_by_any(reactions)
     objective = model.solver.objective if objective is None else objective
+    result = None
     with model:
         model.objective = objective
-        carbon_io = _c_input_output(model, c_source)
+        if c_source is None:
+            c_input = get_c_input(model)
+        else:
+            c_input = model.reactions.get_by_any(c_source)[0]
+        objective_reactions = list(sutil.linear_reaction_coefficients(model))
+        if len(objective_reactions) != 1:
+            raise ValueError('cannot calculate yields for objectives with '
+                             'multiple reactions')
+        carbon_io = c_input, objective_reactions[0]
         min_max = fva(model, reactions, fraction_of_optimum=0)
         grid = [linspace(min_max.minimum[rxn.id], min_max.maximum[rxn.id],
                          points, endpoint=True) for rxn in reactions]
         grid_list = list(product(*grid))
-        result = _envelope_for_points(model, reactions, grid_list, carbon_io)
+        result = envelope_for_points(model, reactions, grid_list, carbon_io)
 
     return pd.DataFrame(result)
 
 
-def _c_input_output(model, c_source=None):
-    if c_source is None:
-        return None, None
-    c_input = _carbon_reaction(model.reactions.get_by_any(c_source))
-    objective_coefficients = sutil.linear_reaction_coefficients(model)
-    c_output = _carbon_reaction(list(objective_coefficients))
-    return c_input, c_output
-
-
-def _carbon_reaction(reactions):
-    exchanges = [rxn for rxn in reactions if len(rxn.reactants) == 1 and
-                 rxn.reactants[0].elements.get('C', 0) > 1]
-    if len(exchanges) != 1:
-        raise ValueError('cannot calculate yields for %s' % reactions)
-    return exchanges[0]
-
-
-def _envelope_for_points(model, reactions, grid, carbon_io):
+def envelope_for_points(model, reactions, grid, carbon_io):
     results = defaultdict(list)
     for direction in ('minimum', 'maximum'):
         sense = "min" if direction == "minimum" else "max"
@@ -416,39 +411,40 @@ def _envelope_for_points(model, reactions, grid, carbon_io):
                     results['direction'].append(direction)
                     results['flux'].append(model.solver.objective.value)
                     if carbon_io[0] is not None:
-                        results['carbon_yield'].append(
-                            _carbon_yield(carbon_io))
-                        results['mass_yield'].append(_mass_yield(carbon_io))
+                        results['carbon_yield'].append(carbon_yield(carbon_io))
+                        results['mass_yield'].append(mass_yield(carbon_io))
     for key, value in results.items():
         results[key] = array(value)
     if carbon_io[0] is not None:
         results['carbon_source'] = carbon_io[0].id
     return results
 
 
-def _carbon_flux(reaction):
-    """Carbon flux for a reaction
-
-    Parameters
-    ----------
-    reaction : cobra.Reaction
-        the reaction to carbon return flux for
+def carbon_yield(c_input_output):
+    """ mol product per mol carbon input
 
     Returns
     -------
     float
-        reaction flux multiplied by number of carbon in reactants"""
-    carbon = sum(metabolite.elements.get('C', 0) for
-                 metabolite in reaction.reactants)
+        the mol carbon atoms in the product (as defined by the model
+        objective) divided by the mol carbon in the input reactions (as
+        defined by the model medium) or zero in case of division by zero
+        arises
+    """
 
+    c_input, c_output = c_input_output
+    if c_input is None:
+        return nan
+    carbon_input_flux = total_carbon_flux(c_input, consumption=True)
+    carbon_output_flux = total_carbon_flux(c_output, consumption=False)
     try:
-        return reaction.flux * carbon
-    except AssertionError:
+        return carbon_output_flux / carbon_input_flux
+    except ZeroDivisionError:
         return nan
 
 
-def _carbon_yield(c_input_output):
-    """Mol product per mol carbon input
+def mass_yield(c_input_output):
+    """Gram product divided by gram of carbon input source
 
     Parameters
     ----------
@@ -459,39 +455,88 @@ def _carbon_yield(c_input_output):
     Returns
     -------
     float
-        the mol carbon atoms in the product (as defined by the model
-        objective) divided by the mol carbon in the input reactions (as
-        defined by the model medium) or zero in case of division by zero
-        arises
+        gram product per 1 g of feeding source
     """
     c_input, c_output = c_input_output
-    carbon_input_flux = _carbon_flux(c_input)
-    carbon_output_flux = _carbon_flux(c_output)
+    if input is None:
+        return nan
     try:
-        return carbon_output_flux / (carbon_input_flux * -1)
-    except ZeroDivisionError:
+        c_source, source_flux = single_flux(c_input, consumption=True)
+        c_product, product_flux = single_flux(c_output, consumption=False)
+    except ValueError:
         return nan
+    mol_prod_mol_src = product_flux / source_flux
+    x = mol_prod_mol_src * c_product.formula_weight
+    return x / c_source.formula_weight
 
 
-def _mass_yield(c_input_output):
-    """Gram product divided by gram of carbon input source
+def total_carbon_flux(reaction, consumption=True):
+    """summed product carbon flux for a reaction
 
     Parameters
     ----------
-    c_input_output : tuple
-        Two reactions, the one that feeds carbon to the system and the one
-        that produces carbon containing compound.
+    reaction : Reaction
+        the reaction to carbon return flux for
+    consumption : bool
+        flux for consumed metabolite, else produced
 
     Returns
     -------
     float
-        gram product per 1 g of feeding source
+        reaction flux multiplied by number of carbon for the products of the
+        reaction
+    """
+    direction = 1 if consumption else -1
+    c_flux = [reaction.flux * coeff * met.elements.get('C', 0) * direction
+              for met, coeff in reaction.metabolites.items()]
+    return sum([flux for flux in c_flux if flux > 0])
+
+
+def single_flux(reaction, consumption=True):
+    """flux into single product for a reaction
+
+    only defined for reactions with single products
+
+    Parameters
+    ----------
+    reaction : Reaction
+        the reaction to product flux for
+    consumption : bool
+        flux for consumed metabolite, else produced
+
+    Returns
+    -------
+    tuple
+        metabolite, flux for the metabolite
+    """
+    if len(list(reaction.metabolites)) != 1:
+        raise ValueError('product flux only defined for single metabolite '
+                         'reactions')
+    met, coeff = next(iteritems(reaction.metabolites))
+    direction = 1 if consumption else -1
+    return met, reaction.flux * coeff * direction
+
+
+def get_c_input(model):
+    """ carbon source reactions
+
+    Returns
+    -------
+    Reaction
+       The medium reaction with highest input carbon flux
     """
-    c_input, c_output = c_input_output
-    source_mass = sum(met.formula_weight for met in c_input.reactants)
-    product_mass = sum(met.formula_weight for met in c_output.reactants)
     try:
-        mol_prod_mol_src = c_output.flux / (c_input.flux * -1)
-    except AssertionError:
-        return nan
-    return (mol_prod_mol_src * product_mass) / source_mass
+        model.solver.optimize()
+        sutil.assert_optimal(model)
+    except OptimizationError:
+        return None
+
+    reactions = model.reactions.get_by_any(list(model.medium))
+    reactions_fluxes = [(rxn, total_carbon_flux(rxn, consumption=True))
+                        for rxn in reactions]
+    source_reactions = [(rxn, c_flux) for rxn, c_flux
+                        in reactions_fluxes if c_flux > 0]
+    try:
+        return max(source_reactions, key=itemgetter(1))[0]
+    except ValueError:
+        return None

cobra/test/test_flux_analysis.py

@@ -699,10 +699,15 @@ def test_envelope_one(self, model):
         df = production_envelope(model, ["EX_o2_e"])
         assert abs(sum(df.flux) - 9.342) < 0.001
 
+    def test_envelope_multi_reaction_objective(self, model):
+        obj = {model.reactions.EX_ac_e: 1,
+               model.reactions.EX_co2_e: 1}
+        with pytest.raises(ValueError):
+            production_envelope(model, "EX_o2_e", obj)
+
     def test_envelope_two(self, model):
         df = production_envelope(model, ["EX_glc__D_e", "EX_o2_e"],
-                                 objective="EX_ac_e",
-                                 c_source="EX_glc__D_e")
+                                 objective="EX_ac_e")
         assert abs(numpy.sum(df.carbon_yield) - 83.579) < 0.001
         assert abs(numpy.sum(df.flux) - 1737.466) < 0.001
         assert abs(numpy.sum(df.mass_yield) - 82.176) < 0.001

release-notes/next-release.md

@@ -6,6 +6,10 @@
   [was broken](https://github.com/opencobra/cobrapy/issues/537)
   following the release of 0.6.0 and has now been fixed (and now with
   unit tests).
+- `production_envelope` failed when model C-source was formulated as
+  -> x instead of x <-. Fixed added option to guess the C-source by
+  taking the medium reaction with the highest input C flux.
+- `model_to_pymatbridge` needs scipy and that's correctly handled now.
 
 ## New features
 - `flux_variability_analysis` now has the `pfba_factor` parameter