def pfba_gapfill(model, reaction_bag, obj=None, obj_lb=10.0, obj_constraint=False, iters=1, tasks=None, task_lb=0.05, add_exchanges=True, extracellular='e', cores=4):
    '\n    Function that utilizes iterations of pFBA solution with a universal reaction bag \n    in order to gapfill a model.\n    \n    Parameters\n    ----------\n    model : cobra.Model\n        Model to be gapfilled\n    reaction_bag : cobra.Model\n        Reaction bag reference to use during gapfilling\n    obj : string\n        Reaction ID for objective function in model to be gapfilled.\n    obj_lb : float\n        Lower bound for objective function\n    obj_constraint : bool\n        Sets objective as contstraint which must be maximized\n    tasks : list or None\n        List of reactions IDs (strings) of metabolic tasks \n        to set a minimum lower bound for\n    task_lb : float\n        Lower bound for any metabolic tasks\n    iters : int\n        Number of gapfilling rounds. Unique reactions from each round are \n        saved and the union is added simulatneously to the model\n    add_exchanges : bool\n        Identifies extracellular metabolites added during gapfilling that\n        are not associated with exchange reactions and creates them\n    extracellular : string\n        Label for extracellular compartment of model\n    cores : int\n        Number of processors to utilize during flux sampling\n    '
    start_time = time.time()
    orig_rxn_ids = set([str(x.id) for x in model.reactions])
    orig_cpd_ids = set([str(y.id) for y in model.metabolites])
    univ_rxn_ids = set([str(z.id) for z in reaction_bag.reactions])
    overlap_rxn_ids = univ_rxn_ids.intersection(orig_rxn_ids)
    if (obj == None):
        obj = get_objective(model)
    else:
        obj = obj
    new_rxn_ids = set()
    print('Creating universal model...')
    with reaction_bag as universal:
        for rxn in overlap_rxn_ids:
            universal.reactions.get_by_id(rxn).remove_from_model()
        if obj_constraint:
            universal.add_reactions([model.reactions.get_by_id(obj)])
            universal.objective = obj
            orig_rxn_ids.remove(obj)
            orig_rxns = []
            for rxn in orig_rxn_ids:
                orig_rxns.append(copy.deepcopy(model.reactions.get_by_id(rxn)))
        else:
            orig_rxns = list(copy.deepcopy(model.reactions))
        add_pfba(universal)
        universal.add_reactions(orig_rxns)
        if (not obj_constraint):
            universal.reactions.get_by_id(obj).lower_bound = obj_lb
        if (tasks != None):
            for task in tasks:
                try:
                    universal.reactions.get_by_id(task).lower_bound = task_lb
                except:
                    print((task + 'not found in model. Ignoring.'))
                    continue
        print('Optimizing model with combined reactions...')
        solution = universal.optimize()
        if (iters > 1):
            print('Generating flux sampling object...')
            sutil.fix_objective_as_constraint(universal, fraction=0.99)
            optgp_object = OptGPSampler(universal, processes=cores)
            print((('Sampling ' + str(iters)) + ' flux distributions...'))
            flux_samples = optgp_object.sample(iters)
            rxns = list(flux_samples.columns)
            for distribution in flux_samples.iterrows():
                for flux in range(0, len(list(distribution[1]))):
                    if (abs(list(distribution[1])[flux]) > 1e-06):
                        new_rxn_ids |= set([rxns[flux]]).difference(orig_rxn_ids)
        else:
            rxns = list(solution.fluxes.index)
            fluxes = list(solution.fluxes)
            for flux in range(0, len(fluxes)):
                if (abs(fluxes[flux]) > 1e-06):
                    new_rxn_ids |= set([rxns[flux]])
    if (obj in new_rxn_ids):
        new_rxn_ids.remove(obj)
    for rxn in orig_rxn_ids:
        try:
            new_rxn_ids.remove(rxn)
        except:
            continue
    print('Gapfilling model...')
    new_rxns = copy.deepcopy([reaction_bag.reactions.get_by_id(rxn) for rxn in new_rxn_ids])
    new_cpd_ids = set()
    for rxn in new_rxns:
        new_cpd_ids |= set([str(x.id) for x in list(rxn.metabolites)])
    new_cpd_ids = new_cpd_ids.difference(orig_cpd_ids)
    new_cpds = copy.deepcopy([reaction_bag.metabolites.get_by_id(cpd) for cpd in new_cpd_ids])
    new_model = copy.deepcopy(model)
    new_model.add_metabolites(new_cpds)
    new_model.add_reactions(new_rxns)
    if (add_exchanges == True):
        new_exchanges = extend_exchanges(new_model, new_cpd_ids, extracellular)
        if (len(new_exchanges) > 0):
            new_rxn_ids |= new_exchanges
    duration = int(round((time.time() - start_time)))
    print((((((('Took ' + str(duration)) + ' seconds to gapfill ') + str(len(new_rxn_ids))) + ' reactions and ') + str(len(new_cpd_ids))) + ' metabolites.'))
    new_obj_val = new_model.slim_optimize()
    if (new_obj_val > 1e-06):
        print((('Gapfilled model objective now carries flux (' + str(new_obj_val)) + ').'))
    else:
        print('Gapfilled model objective still does not carry flux.')
    return new_model
Function that utilizes iterations of pFBA solution with a universal reaction bag 
in order to gapfill a model.

Parameters
----------
model : cobra.Model
    Model to be gapfilled
reaction_bag : cobra.Model
    Reaction bag reference to use during gapfilling
obj : string
    Reaction ID for objective function in model to be gapfilled.
obj_lb : float
    Lower bound for objective function
obj_constraint : bool
    Sets objective as contstraint which must be maximized
tasks : list or None
    List of reactions IDs (strings) of metabolic tasks 
    to set a minimum lower bound for
task_lb : float
    Lower bound for any metabolic tasks
iters : int
    Number of gapfilling rounds. Unique reactions from each round are 
    saved and the union is added simulatneously to the model
add_exchanges : bool
    Identifies extracellular metabolites added during gapfilling that
    are not associated with exchange reactions and creates them
extracellular : string
    Label for extracellular compartment of model
cores : int
    Number of processors to utilize during flux sampling
def pfba_gapfill(model, reaction_bag, obj=None, obj_lb=10.0, obj_constraint=False, iters=1, tasks=None, task_lb=0.05, add_exchanges=True, extracellular='e', cores=4):
    '\n    Function that utilizes iterations of pFBA solution with a universal reaction bag \n    in order to gapfill a model.\n    \n    Parameters\n    ----------\n    model : cobra.Model\n        Model to be gapfilled\n    reaction_bag : cobra.Model\n        Reaction bag reference to use during gapfilling\n    obj : string\n        Reaction ID for objective function in model to be gapfilled.\n    obj_lb : float\n        Lower bound for objective function\n    obj_constraint : bool\n        Sets objective as contstraint which must be maximized\n    tasks : list or None\n        List of reactions IDs (strings) of metabolic tasks \n        to set a minimum lower bound for\n    task_lb : float\n        Lower bound for any metabolic tasks\n    iters : int\n        Number of gapfilling rounds. Unique reactions from each round are \n        saved and the union is added simulatneously to the model\n    add_exchanges : bool\n        Identifies extracellular metabolites added during gapfilling that\n        are not associated with exchange reactions and creates them\n    extracellular : string\n        Label for extracellular compartment of model\n    cores : int\n        Number of processors to utilize during flux sampling\n    '
    start_time = time.time()
    orig_rxn_ids = set([str(x.id) for x in model.reactions])
    orig_cpd_ids = set([str(y.id) for y in model.metabolites])
    univ_rxn_ids = set([str(z.id) for z in reaction_bag.reactions])
    overlap_rxn_ids = univ_rxn_ids.intersection(orig_rxn_ids)
    if (obj == None):
        obj = get_objective(model)
    else:
        obj = obj
    new_rxn_ids = set()
    print('Creating universal model...')
    with reaction_bag as universal:
        for rxn in overlap_rxn_ids:
            universal.reactions.get_by_id(rxn).remove_from_model()
        if obj_constraint:
            universal.add_reactions([model.reactions.get_by_id(obj)])
            universal.objective = obj
            orig_rxn_ids.remove(obj)
            orig_rxns = []
            for rxn in orig_rxn_ids:
                orig_rxns.append(copy.deepcopy(model.reactions.get_by_id(rxn)))
        else:
            orig_rxns = list(copy.deepcopy(model.reactions))
        add_pfba(universal)
        universal.add_reactions(orig_rxns)
        if (not obj_constraint):
            universal.reactions.get_by_id(obj).lower_bound = obj_lb
        if (tasks != None):
            for task in tasks:
                try:
                    universal.reactions.get_by_id(task).lower_bound = task_lb
                except:
                    print((task + 'not found in model. Ignoring.'))
                    continue
        print('Optimizing model with combined reactions...')
        solution = universal.optimize()
        if (iters > 1):
            print('Generating flux sampling object...')
            sutil.fix_objective_as_constraint(universal, fraction=0.99)
            optgp_object = OptGPSampler(universal, processes=cores)
            print((('Sampling ' + str(iters)) + ' flux distributions...'))
            flux_samples = optgp_object.sample(iters)
            rxns = list(flux_samples.columns)
            for distribution in flux_samples.iterrows():
                for flux in range(0, len(list(distribution[1]))):
                    if (abs(list(distribution[1])[flux]) > 1e-06):
                        new_rxn_ids |= set([rxns[flux]]).difference(orig_rxn_ids)
        else:
            rxns = list(solution.fluxes.index)
            fluxes = list(solution.fluxes)
            for flux in range(0, len(fluxes)):
                if (abs(fluxes[flux]) > 1e-06):
                    new_rxn_ids |= set([rxns[flux]])
    if (obj in new_rxn_ids):
        new_rxn_ids.remove(obj)
    for rxn in orig_rxn_ids:
        try:
            new_rxn_ids.remove(rxn)
        except:
            continue
    print('Gapfilling model...')
    new_rxns = copy.deepcopy([reaction_bag.reactions.get_by_id(rxn) for rxn in new_rxn_ids])
    new_cpd_ids = set()
    for rxn in new_rxns:
        new_cpd_ids |= set([str(x.id) for x in list(rxn.metabolites)])
    new_cpd_ids = new_cpd_ids.difference(orig_cpd_ids)
    new_cpds = copy.deepcopy([reaction_bag.metabolites.get_by_id(cpd) for cpd in new_cpd_ids])
    new_model = copy.deepcopy(model)
    new_model.add_metabolites(new_cpds)
    new_model.add_reactions(new_rxns)
    if (add_exchanges == True):
        new_exchanges = extend_exchanges(new_model, new_cpd_ids, extracellular)
        if (len(new_exchanges) > 0):
            new_rxn_ids |= new_exchanges
    duration = int(round((time.time() - start_time)))
    print((((((('Took ' + str(duration)) + ' seconds to gapfill ') + str(len(new_rxn_ids))) + ' reactions and ') + str(len(new_cpd_ids))) + ' metabolites.'))
    new_obj_val = new_model.slim_optimize()
    if (new_obj_val > 1e-06):
        print((('Gapfilled model objective now carries flux (' + str(new_obj_val)) + ').'))
    else:
        print('Gapfilled model objective still does not carry flux.')
    return new_model