#215 GrowthCouplingPotential

                r.id if isinstance(r, cobra.Reaction) else r for r in exclude_reactions

    """

    Optimize the Growth Coupling Potential through a combination of knockouts, knockins and medium additions

    Use the 'run' method to optimize the MILP problem and find the best solution.

    In order to find multiple sub-optimal solutions through the use of solution pools, please refer to the

    relevant solvers documentation. The underlying solver object can be accessed through

        GrowthCouplingPotential(...).model.solver.problem

    The method is described in the paper "OptCouple: Joint simulation of gene knockouts,

    insertions and medium modifications for prediction of growth-coupled strain designs"

    https://www.sciencedirect.com/science/article/pii/S2214030118300373

    Attributes

    ----------

    model : cobra.Model

        A cobra model object containing the target reaction as well as all native and heterologous reactions.

    target : cobra.Reaction or string

        The reaction (or reaction ID) that is to be growth-coupled

    knockout_reactions : list

        A list of ID's for the reactions that the algorithm should attempt to knock out

    knockin_reactions : list

        A List of ID's for the reactions that are heterologous and that the algorithm should attempt to add

    medium_additions : list

        A list of metabolite ID's for the the metabolites that should potentially be added to the medium

    n_knockouts : int

        The maximum allowed number of knockouts

    n_knockins : int

        The maximum allowed number of knockins

    n_medium : int

        The maximum allowed number of medium additions

    use_solution_pool : Bool

        Set to True when using the Gurobi solver to identify multiple alternative solutions at once

    """

            self, model, target, biomass_id, knockout_reactions, knockin_reactions, medium_additions,

            n_knockouts, n_knockin, n_medium, remove_blocked=True, use_solution_pool=False, **kwargs

    ):

        # Check model for forced fluxes, which should not be present. The zero solution must be allowed.

                    "%s has a forced flux. Remove the reaction or change its bounds for the function to work." % r.id)

        # If a cobra.Reaction is given as target, convert to its ID string.

        # Add boundary reactions for these metabolites named "ADD_<metabolitename>" to the model, with bounds (-10,0).

            # r.type = "medium"

        # Set the solver to Gurobi for the fastest result. Set to CPLEX if Gurobi is not available.

            # The tolerances are set to the minimum value. This gives maximum precision.

                "Changing solver to CPLEX, as Gurobi is not available."

                "This may cause a big slowdown and limit options afterwards.")

            # The tolerances are set to the minimum value. This gives maximum precision.

        else:

                           model.solver.interface.__name__.split(".")[-1])

        # Remove reactions that are blocked: no flux through these reactions possible.

        # This will reduce the search space for the solver, if not done already.

        # Make dual

        # Combine primal and dual

        # Dictionaries to hold the binary control variables:

        # Now the fun stuff

        # Fill the dictionaries with binary variables

        # For the knockouts: only for the native reactions which are not exchanges

            # Add constraint variables

            # Constrain the primal: flux through reactions within (-1000, 1000)

                                                  name="primal_y_const_" + reaction.id + "_ub"))

                                                  name="primal_y_const_" + reaction.id + "_lb"))

            # Constrain the dual

            # Add to dictionary with binary variables for native reactions:

        # For the knockins and medium additions:

            # Add constraint variables

            # Constrain the primal: flux through reaction within (-1000, 1000)

                                                  name="primal_y_const_" + reaction.id + "_ub"))

                                                  name="primal_y_const_" + reaction.id + "_lb"))

            # Constrain the dual

            # Add y variable to the corresponding modifications dictionary

            else:

        # Set the objective

            dual_problem.objective, model=model.solver

        )

            *((c * v) for v, c in dual_objective.expression.as_coefficients_dict().items()

              if v not in constrained_dual_vars)

        )

        # Add number of knockouts constraint. ub=K+1 as the target will be forced to be 1 as well

            optlang.symbolics.Add(*native_y_vars), lb=0, ub=n_knockouts + 1, name="number_of_knockouts_constraint"

        )

        # Add number of knockins constraint

            optlang.symbolics.Add(*heterologous_y_vars), lb=0, ub=n_knockin, name="number_of_knockins_constraint"

        )

        # Add number of medium additions constraint

            optlang.symbolics.Add(*medium_y_vars), lb=0, ub=n_medium, name="number_of_medium_additions_constraint"

        )

        # Force target control variable to be 1, but allow flux in primal

                raise NotImplementedError("Solution pools are currently only available with the Gurobi solver")

            # Find the N best solutions (N is specified as pool_size in run())

            # Don't store solutions with objective = 0

        """

        Run the GrowthCouplingPotential method

        pool_size: None or int

            The number of alternative solutions to find (0 - 2,000,000,000)

        pool_time_limit: None or int

            The maximum time (in seconds) to spend finding solutions. Will terminate with a time_limit status.

        """

        else:

                reac.id for var, reac in self.native_y_vars.items() if var.primal > 0.9 and reac.id != self.target

            ]

                "knockouts": knockouts,

                "knockins": knockins,

                "medium": medium_additions,

                "obj_val": objective_value

            }

        """

        Extracts the individual solutions from an MILP solution pool created by the Gurobi solver.

        Returns the solution pool and a header with information to accompany the list of dicts.

        """

                "The solver for the MILP is not Gurobi. This function only works for Gurobi created solution pools.")

            # Note that the objective is the growth coupling potential, not the growth or production rate

            # Prepare dict for looping over variables

                    else:

                # Extract target flux and growth rate. Subtract the forward and reverse fluxes to get the net flux.

    def write_solutions_to_file(solutions, filename):