Find the variability of each reaction while fixing certain fluxes.

    Yields the reaction id, and a tuple of minimum and maximum value for each
    of the given reactions. The fixed reactions are given in a dictionary as
    a reaction id to value mapping.

    This is an implementation of flux variability analysis (FVA) as described
    in [Mahadevan03]_.

    Args:
        model: MetabolicModel to solve.
        reactions: Reactions on which to report variablity.
        fixed: dict of additional lower bounds on reaction fluxes.
        tfba: If True enable thermodynamic constraints.
        solver: LP solver instance to use.

    Returns:
        Iterator over pairs of reaction ID and bounds. Bounds are returned as
        pairs of lower and upper values.

def flux_variability(model, reactions, fixed, tfba, solver):
    """Find the variability of each reaction while fixing certain fluxes.

    Yields the reaction id, and a tuple of minimum and maximum value for each
    of the given reactions. The fixed reactions are given in a dictionary as
    a reaction id to value mapping.

    This is an implementation of flux variability analysis (FVA) as described
    in [Mahadevan03]_.

    Args:
        model: MetabolicModel to solve.
        reactions: Reactions on which to report variablity.
        fixed: dict of additional lower bounds on reaction fluxes.
        tfba: If True enable thermodynamic constraints.
        solver: LP solver instance to use.

    Returns:
        Iterator over pairs of reaction ID and bounds. Bounds are returned as
        pairs of lower and upper values.
    """

    fba = _get_fba_problem(model, tfba, solver)

    for reaction_id, value in iteritems(fixed):
        flux = fba.get_flux_var(reaction_id)
        fba.prob.add_linear_constraints(flux >= value)

    def min_max_solve(reaction_id):
        for direction in (-1, 1):
            yield fba.flux_bound(reaction_id, direction)

    # Solve for each reaction
    for reaction_id in reactions:
        yield reaction_id, tuple(min_max_solve(reaction_id))