def print_product_reactions_by(self,ref_type):
        """
        Prints in the output the list of all reactions this metabolite participates in 
        as a product in the model with the format specified by ref_type 
 
        INPUTS:
        -------
        ref_type: A string indicating the in what format the metabolites should be 
                  printed in the output. Current eligible choices are 'id', 'name'
                  and formula. If name or id was not provided, id used instead.
        """
        if ref_type.lower() not in ['id','name','equation']:
            raise userError("**Error! Invalid reference type (eligible choices are 'id' and 'name')")
        elif self.reactions == []:
            raise userError('**Error! List of the reactions is not defined for this metabolite ...')

        for reaction in self.product_reactions:
            if ref_type.lower() == 'id':
                rxn = reaction.id
            elif ref_type.lower() == 'name':
                if reaction.name is None:
                    rxn = reaction.id
                    userWARNING("Waarning! No name was provided for '" + reaction.id + "'. id is used instead.")
                else:
                    rxn = reaction.name
            elif ref_type.lower() == 'equation':
                if reaction.equation is None:
                    rxn = reaction.id
                    userWARNING("Waarning! No equation was provided for '" + reaction.id + "'. id is used instead.")
                else:
                    rxn = reaction.equation

            print rxn

    def print_reactions_by(self,ref_type, metab_ref = None):
        """
        Prints in the output the list of all reactions this metabolite participates in 
        as a reactant or product in the model with the format specified by ref_type 
 
        INPUTS:
        -------
         ref_type: A string indicating  in what format the reactions should be 
                   printed in the output. Current eligible choices are 'id', 'name'
                   and formula. If name or id was not provided, id used instead.
        metab_ref: A string indicating in what format the metabolites should appear
                   in a reaciton equation if ref_type = 'equation'. If None, metabolite
                   ids are used
        """
        if ref_type.lower() not in ['id','name','equation']:
            raise userError("Invalid reference type (eligible choices are 'id' and 'name')")
        elif self.reactions == None:
            raise userError('List of the reactions is not defined for this metabolite ...')

        for reaction in self.reactions:
            if ref_type.lower() == 'id':
                rxn = reaction.id
            elif ref_type.lower() == 'name':
                if reaction.name is None:
                    rxn = reaction.id
                    userWARNING("Waarning! No name was not provided for '" + reaction.id + "'. id is used instead.")
                else:
                    rxn = reaction.name
            elif ref_type.lower() == 'equation':
                if metab_ref != None:
                    rxn = reaction.get_eqn_by(metab_ref)
                else:
                    rxn = reaction.get_eqn_by('id')

            print rxn 

    def __init__(self,model, optSolverName = None, createModel = None, screenOutput = None): 

        """
        INPUTS (required):
        ------
                      model: An instance of class model containing the information
                             about the metabolic model

        INPUTS (optional):
        ------
              optSolverName: Name of the LP solver to be used to solve the LP. Current 
                             allowable choices are cplex and gurobi
               screenOutput: By default (on) writes  a summary including the solve 
                             status, optimality status (if not optimal), objective 
                             function value and the elapsed time on the screen.
                             if set to a value of 'off' no resuults are written on 
                             the screen, in which case The user can instead specifiy 
                             an output fiile using the option outputFile, or store 
                             them in the variable runOutput (see the 'run' method for
                             details)
                createModel: A parameter indicating whether a pyomo model should be 
                             created (1) or not. Default is 1. The options is useful
                             for the cases a model is already created and one just 
                             wants to change some model attributes (e.g., flux bounds)
                             and rerun FBA. Setting this parameter to zero will save 
                             soem runtime as the model need not to be created again.
        """
       
        # Metabolic model
        self.model = model

        # Solver name
        if optSolverName == None:
            self.optSolverName = 'cplex'
        else:
            if optSolverName.lower() in ['cplex','gurobi']:
                self.optSolverName = optSolverName
            else:
                raise userError('**Error! Invalid solver name (eligible choices are cplex and gurobi)\n')          

        # Whether to create a pyomo model
        if createModel == None: 
            self.createModel = 1   
        elif createModel not in [0,1]:    
            raise userError('**Error! Inalid value for createModel! The allowable values are 0 and 1 ')
        else:
            self.createModel = createModel
               
        # Output to the screen 
        if screenOutput == None:
            self.screenOutput = 'on'
        elif type(screenOutput) is not str:
            raise userError("**Error! screenOutput should be a string ('on' or 'off')")
        elif screenOutput.lower() not in ['on','off']:
            raise userError("**Error! The only eligible values for screenOutput are 'on' and 'off'")
        else:
             self.screenOutput = screenOutput

    def biomass_yield_calc(self,fba_model = None):
        """
        Calculates the biomass yield of this metabolite (if applicable).
        For the yield to be computed this metabolite should particiapte in
        in an exchange reaction (EX_m(e): m[e] <==>)

        INPUTS:
        -------
        fba_model: An instance of the class fba containing the fba model for the 
                   wild-type organism under the desired growth condition. If this is 
                   not provided, then the current fba_model of model is used 

        OUTPUS:
        -------
        The output is assigned to the global variable self.biomass_yield 
        """
        # Find the exchange reaction this metabolite participates in
        exch_rxn = [r for r in self.reactant_reactions if r.type.lower() == 'exchange']
        if len(exch_rxn) == 0:
            raise userError('**ERROR! Unable to compute the biomass yield for metabolite ' + self.id +'. The metabolite must have an excange reaction in reactant_reactions ...')
        elif len(exch_rxn) > 1:
             raise userError('**ERROR! metabolite ' + self.id + ' participates in more than one exchange reaciton ...')
        else:
            exch_rxn = exch_rxn[0]

            # Original LB on reaction flux 
            exch_rxn_LB = exch_rxn.flux_bounds[0]
  
            # Perform FBA for 10 mole uptake of this metabolite 
            exch_rxn.flux_bounds[0] = -10      

            if fba_model == None:
                # WARNING! If no fba_model is provided, the current model.fba_model is 
                # used to compute the biomass yield
                self.model.fba(create_model = False,store_opt_fluxes = True, screen_output = 'off')
                fba_solution = self.model.fba_model.solution
            else:
                fba_model.store_opt_fluxes = True
                fba_model.screen_output = 'off' 
                fba_model.run()
                fba_solution = fba_model.solution
 
            if fba_solution['exit_flag'] == 'globallyOptimal':
                if exch_rxn.flux < 0:
                    self.biomass_yield = fba_solution['objective_value']/(-exch_rxn.flux)
                else:
                    # if zero (no uptake) or positive (production)
                    self.biomass_yield = None 
                    print '\n**WARNING! The fba problem to compute the biomass yield for metabolite ' + self.id + ' resulted in a non-negative flux value for the corresponding exchange reaction. No value is assigned to biomass_yield'
            else:
                self.biomass_yield = None
                print '\n**WARNING! The fba problem to compute the biomass yield for metabolite ' + self.id + ' was not solved to optimality. No value is assigned to biomass_yield'

            # Set the lower bound on exchange reaction back to what it was before
            exch_rxn.flux_bounds[0] = exch_rxn_LB

    def get_compounds(self,ref_type = 'id'):
        """
        Reports the list of compounds participating in this reaction in the output 
        with the format specified by ref_type 
 
        INPUTS:
        -------
        ref_type: A string indicating the in what format the compounds should be 
                  printed in the output. Current eligible choices are 'id', 'name'
                  and formula. If name or id was not provided, id used instead.
        """
        if ref_type.lower() not in ['id','name','formula']:
            raise userError("**Error! Invalid reference type (eligible choices are 'id' and 'name')")

        for metab in self.compounds:
            if ref_type.lower() == 'id':
                mm = metab.id
            elif ref_type.lower() == 'name':
                if metab.name is None:
                    mm = metab.id
                else:
                    mm = metab.name
            elif ref_type.lower() == 'formula':
                if metab.formula is None:
                    mm = metab.id
                else:
                    mm = metab.formula

            return mm

    def get_products_by(self,ref_type):
        """
        Reports the list of products of this reaciton in the output 
        with the format specified by ref_type 
 
        INPUTS:
        -------
        ref_type: A string indicating the in what format the metabolites should be 
                  printed in the output. Current eligible choices are 'id', 'name'
                  and formula. If name or id was not provided, id used instead.
        """
        if ref_type.lower() not in ['id','name']:
            raise userError("**Error! Invalid reference type (eligible choices are 'id' and 'name')")

        for metab in self.products: 
            if ref_type.lower() == 'id':
                mm = metab.id
            elif ref_type.lower() == 'name':
                if metab.name is None:
                    mm = metab.id
                else:
                    mm = metab.name
            elif ref_type.lower() == 'formula':
                if metab.formula is None:
                    mm = metab.id
                else:
                    mm = metab.formula

            print mm

    def assign_props(self):
        """
        Assigns the properties of reaction
        """
        # A list of the form [dGmin,dGmax] containing the min and max values of deltaG
        # (Gibbs free energy change) for this reaction
        if self.deltaG_range == [] and self.deltaG != None and self.deltaG_uncertainty != None:
            self.deltaG_range = [min(self.deltaG - self.deltaG_uncertainty,self.deltaG + self.deltaG_uncertainty),max(self.deltaG - self.deltaG_uncertainty,self.deltaG + self.deltaG_uncertainty)]

        # List of compound objects containing all compounds participating in this reaction
        self.compounds = sorted(list(set([c for c in self.stoichiometry.keys() if self.stoichiometry[c] < 0] + [c for c in self.stoichiometry.keys() if self.stoichiometry[c] > 0])),key=lambda x:x.id) 

        # List of compound objects containing all reactants of this reaction
        self.reactants = sorted(list(set([m for m in self.stoichiometry.keys() if self.stoichiometry[m] < 0])),key=lambda x:x.id)

        # List of compound objects containing all products of this reaction
        self.products = sorted(list(set([m for m in self.stoichiometry.keys() if self.stoichiometry[m] > 0])),key=lambda x:x.id)

        # This is particularly hepful when a reaction occurs in more than one compartment
        # Here, we have a list instead of a single element because different compounds
        # in the reaction may participate in different compartment
        if self.compartment == []:
            self.compartment = list(set([c.compartment for c in self.compounds]))

        # model
        models = list(set([c.model for c in self.compounds if hasattr(c,'model') and c.model != None]))
        if len(models) == 1:
           self.model = models[0]
        elif len(models) > 1:
            raise userError('More than one model assigned to "compounds" of reaction ' + self.id + ': ' + str([m.id for m in models]))
        if len(models) == 0:
           self.model = None 

    def get_compounds(self,ref_type = 'id'):
        """
        Reports the list of compounds participating in this reaction in the output 
        with the format specified by show_cpds_by 
 
        INPUTS:
        -------
        show_cpds_by: A string indicating the in what format the compounds should be 
                  printed in the output. Current eligible choices are 'id', 'name', 'formula',
                  'ModelSEED_id', 'KEGG_id', 'BiGG_id'. If ModelSEED_id, KEGG_id or BiGG_id ir not
                  available for any compound participating in the reaction, it is replaced with 
                  its id.
        """
        if show_cpds_by.lower() not in ['id','name','formula','ModelSEED_id', 'KEGG_id', 'BiGG_id']:
            raise userError("**Error! Invalid reference type (eligible choices are 'id' and 'name')")

        for metab in self.compounds:
            if show_cpds_by.lower() == 'id':
                mm = metab.id
            elif show_cpds_by.lower() == 'name':
                if metab.name is None:
                    mm = metab.id
                else:
                    mm = metab.name
            elif show_cpds_by.lower() == 'formula':
                if metab.formula is None:
                    mm = metab.id
                else:
                    mm = metab.formula

            return mm

    def assign_flux_bounds(self, assignLB = True, assignUB = True):
        """
        Assigns general bounds to fluxes based on the reaction type 
        assignLB and assignUB indicate whether to assign the computed lower and upper bounds
        """
        if len(self.flux_bounds) == 2:
            flux_LB = self.flux_bounds[0]
            flux_UB = self.flux_bounds[1]
        else:
            flux_LB = None
            flux_UB = None

        if (not self.is_exchange) and self.reversibility.lower() in ['irreversible','irreversible_forward']:
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 1000
        # An irreversible reaction written in the backward direction, e.g., B <-- A
        elif (not self.is_exchange) and self.reversibility.lower() == 'irreversible_backward':
            if assignLB:
                flux_LB = -1000
            if assignUB:
                flux_UB = 0
        elif (not self.is_exchange) and self.reversibility.lower() == 'reversible':
            if assignLB:
                flux_LB = -1000
            if assignUB:
                flux_UB = 1000
        elif (not self.is_exchange) and self.reversibility.lower() == 'reversible_forward':
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 1000
        # Reverse part of a reversible reaciton wirttein the forward direction. For example, if A <==> B, 
        # A --> B is reversible_forward and B --> A is reversible_backward
        elif (not self.is_exchange) and self.reversibility.lower() == 'reversible_backward':
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 1000
        elif self.is_exchange and (self.reversibility.lower() == 'reversible' or self.reversibility.lower() == 'exchange'):
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 1000
        elif self.is_exchange and self.reversibility.lower() == 'rversible_forward':
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 1000
        elif self.is_exchange and self.reversibility.lower() == 'rversible_backward':
            if assignLB:
                flux_LB = 0
            if assignUB:
                flux_UB = 0
        else:
            raise userError('Unknown reaction reversibility: {} for reaction {} with is_exchange = {}'.format(self.reversibility, self.id, self.is_exchange))

        self.flux_bounds = [flux_LB,flux_UB]

 def objectiveFunc_rule(self,fbaModel):
     # Reactions for which the objective coefficient has not bee assigned
     non_obj_rxns = [j.id for j in fbaModel.J if j.objective_coefficient == None]
     if len(non_obj_rxns) >= 1: 
         print("**Error! 'objective_coefficient' has not been defined for the following reacitons:")
         print non_obj_rxns
         raise userError()
     return sum(j.objective_coefficient*fbaModel.v[j] for j in fbaModel.J)

    def __init__(self, id, name = '', name_aliases = [], domain = '', genus = '', species = '', strain = '', model_id = None, ModelSEED_type = '', gDW_per_cell = None, gWW_per_cell = None, cells_per_ml = None, gDW_per_ml = None, mu = None, random_mortality_rate = None, notes = ''): 
    
        # Organism id 
        self.id = id

        # organism complete name (case insensitive string)
        self.name = name

        # Name aliases
        self.name_aliases = name_aliases

        # Domain (case insensitive string).Example of choices are:  
        # bacteria, archaea,Eukaryotes 
        self.domain = domain

        # Genus (case insensitive string). Example: Escherichia 
        self.genus = genus

        # Species ((case insensitive string)). Example: coli
        self.species = species

        # Strain ((case insensitive string)). Example: MG1655
        self.strain = strain

        # model id
        self.model_id = model_id

        # ModelSEED type (bacteria_GramPositive, bacteria_GramNegative, Human and Plant)
        self.ModelSEED_type = ModelSEED_type

        # Gram of dry weight per (one) cell
        self.gDW_per_cell = gDW_per_cell

        # Gram of wet weight per (one) cell
        self.gWW_per_cell = gWW_per_cell

        # Cell mass concnetration in number of cells per ml of the culture
        self.cells_per_ml = cells_per_ml

        # Cell mass concentration in gram of dry weight per ml of the culture
        self.gDW_per_ml = gDW_per_ml

        # Specific growth rate (in 1/h)
        self.mu = mu

        # Mortality rate
        if random_mortality_rate > 0:
            raise userError('Mortality rate for organism ' + self.id + ' must be non-positive')
        else:
            self.random_mortality_rate = random_mortality_rate

        # Notes and comments
        if isinstance(notes,str):
            self.notes = notes
        else:
            self.notes = ''

def find_coopr_exchrxns_fluxRanges(results_filename = 'coopr_exchrxns_fluxRanges.py', warnings = True, stdout_msgs = True):
    """
    Finds the min and max value of each compound that can be produced by the wild-type strain
    """
    from auxoMetabs import auxoMetabsMutants

    # Model path
    model_path = '/usr2/postdoc/alizom/work/models/Escherichia_coli/iJO1366/'

    # Exchange reactions needed for cooperaiton in all pairs
    cooperative_exchrxn_ids = list(set([r for m in auxoMetabsMutants.keys() for r_list in auxoMetabsMutants[m] for r in r_list]))
  
    WT = create_model(warnings = warnings, stdout_msgs = stdout_msgs)

    # Growth medium
    set_specific_bounds(model = WT, file_name = model_path + 'iJO1366_minimal_glucose_anaerobic.py',limiting_nutrients = {'EX_glc(e)':[-10,1000]})

    with open(results_filename,'w') as f:
        f.write('coopr_exchrxns_fluxRanges = {}\n')

    for exchrxn_id in cooperative_exchrxn_ids:

        if stdout_msgs:
            print '\n------- ',exchrxn_id,' ---------\n'

        exchrxn = WT.reactions_by_id[exchrxn_id]

        # Find maximum
        for rxn in WT.reactions:
            rxn.objective_coefficient = 0
        exchrxn.objective_coefficient = 1
 
        WT.fba(store_opt_fluxes = False, warnings = warnings, stdout_msgs = stdout_msgs)
        if WT.fba_model.solution['exit_flag'] == 'globallyOptimal':
             max_exch = WT.fba_model.solution['objective_value']
        else:
             max_exch = 0

        # Find minimum
        for rxn in WT.reactions:
            rxn.objective_coefficient = 0
        exchrxn.objective_coefficient = -1
 
        WT.fba(store_opt_fluxes = False, warnings = warnings, stdout_msgs = stdout_msgs)
        if WT.fba_model.solution['exit_flag'] == 'globallyOptimal':
             min_exch = -WT.fba_model.solution['objective_value']
        else:
             min_exch = 0

        if min_exch > max_exch:
            raise userError('min_exch = {} > max_exch = {}'.format(min_exch,max_exch))

        with open(results_filename,'a') as f:
            f.write("coopr_exchrxns_fluxRanges['" + exchrxn_id + "'] = (" + str(min_exch) + ',' + str(max_exch) + ')\n') 

    print '\nResults were written into coopr_exchrxns_fluxRanges.py ...\n'

def integrate_results_files(results_filenames, results_var_name, mutant_pair, output_file_name):
    """
    Integrates the results of all scripts (when splitting one job to multiple smaller jobs) inot one sinble file 
 
    INPUTS:
    -------
    results_file_names: A list of strings containing the names of the files containing the results
      results_var_name: Name of the variable storing the results
      output_file_name: Name of the output file name containing the integration of all results
           mutant_pair: A tuple containiing the names of the mutant pairs
    """
    from imp import load_source

    # sum of the entries in all result files
    entries_sum = 0

    results = {}
    results[mutant_pair] = {}
    # Import the data in the module stored in file_name
    for file_name in results_filenames:
        if type(file_name) == str:
            if not os.path.isfile(file_name):
                raise IOError("No such file was found :'" + file_name + "'")
            else:
                # First delete the model dataFile if it already exists. If it is not deleted
                # the new module data is merged with the previous ones
                try:
                    del sys.modules['dataFile']
                except:
                    pass
                load_source('dataFile',file_name)
                import dataFile
        exec('results_curr = dataFile.' + results_var_name)
        
        entries_sum += len(results_curr[mutant_pair].keys())
        print 'Total # of entries in ', file_name,' = ',len(results_curr[mutant_pair].keys())

        results_keys = results[mutant_pair].keys()
        for k in results_curr[mutant_pair].keys():
            if k in results_keys:
                raise userError(str(k) + ' in ' + file_name + ' already in results_keys\n')
            else:
                results[mutant_pair][k] = results_curr[mutant_pair][k]

    print '\nThe total # of entries in results = {},  entries_sum = {} '.format(len(results[mutant_pair].keys()),entries_sum)

    # Write the results inot the specified output file
    with open(output_file_name,'w') as f:
        f.write(results_var_name + ' = {}\n')
        f.write(results_var_name + '[' + str(mutant_pair) + '] = {}\n')
        for k in results[mutant_pair].keys():
            f.write(results_var_name + '[' + str(mutant_pair) + '][' + str(k) + '] = ' + str(results[mutant_pair][k]) + '\n')
    print '\nResults were integrated and written into ',output_file_name,'\n'

def create_job_files(total_comb_num,interval_size, job_filename_base, joboutput_filename_base, max_walltime, coopr_level_m1_str, coopr_level_m2_str, glc_conc, glc_excess, use_DMMM_coopr):
    """
    Creates the job files given:
             total_comb_num: Total number of combinations 
              interval_size: The desired of iteration intervals
          outfile_base_name: The base name of the output files
                       task: The task that should be performed ('analyze_games' or 'analyze_cost')
          job_filename_base: Base name of the output job file
    joboutput_filename_base: The name of the file storing the dynamic output of the job. It is essentially the same as job_filename_base
                             with the file path deleted (in order to create .out files in the same directory as the job file)
    """
    if not glc_excess:
        if glc_conc == 111.01:
            results_filename_base = 'results/coopr_level_expGlc'
        elif glc_conc == 1000:
            results_filename_base = 'results/coopr_level_fixedGlc'
        else:
            raise userError('Unknonw glc_conc:' + str(glc_conc))
    else:
        results_filename_base = 'results/coopr_level_excessGlc'

    coopr_level_m1 = range(0,101,5)

    slices = create_intervals(total_comb_num = total_comb_num,interval_size = interval_size)

    for slice in slices:
        job_filename = job_filename_base + '_' + str(slice[0]) + '_' + str(slice[1]) + '.sh'
        joboutput_filename = joboutput_filename_base + '_' + str(slice[0]) + '_' + str(slice[1]) + '.out'
        print 'creating file ',job_filename,'...'
        with open(job_filename,'w') as outfile:
            outfile.write('#!/bin/bash -l\n')
            outfile.write('#$-l h_rt=' + str(int(max_walltime))+ ':00:00\n\n')
            outfile.write('cd /usr2/postdoc/alizom/work/EcoliPairs\n\n')

            # This is to merge the .e[jobid] and .o[jobid] files
            outfile.write('#$ -j y\n')

            # Set the output file
            outfile.write('\n#$ -o ' + joboutput_filename + '\n')

            outfile.write('\nsource activate /projectnb/bioinfor/SEGRE/alizom\n\n')

            # python -c "import time;print '\n**Job started at ',time.strftime('%c'),'\n'" > job_dynamic_yeast_stsrt_pos_end_pos.out 2>&1
            outfile.write("\npython -c \"import time;print '\\n**Job started at ',time.strftime('%c'),'\\n'\"\n")

            start_pos = slice[0]
            end_pos = slice[1]

            outfile.write('\npython -c \"from coopr_level import run_master_func;run_master_func(start_pos_pair = ' + str(start_pos) + ', end_pos_pair = ' + str(end_pos) + ', coopr_level_m1 = ' + coopr_level_m1_str + ', coopr_level_m2 = ' + coopr_level_m2_str + ', glc_conc = ' + str(glc_conc) + ", glc_excess = " + str(glc_excess) + ", use_DMMM_coopr = " + str(use_DMMM_coopr) + ", run_main = True, run_test = False, results_filename_base = '" + results_filename_base + "')\"\n\n")

            # python -c "import time;print '\n**Job ended at ',time.strftime('%c'),'\n'" >> job_dynamic_yeast_start_pos_end_pos.out 2>&1
            outfile.write("\npython -c \"import time;print '\\n**Job ended at ',time.strftime('%c'),'\\n'\"\n")

    def set_products(self, products):
        """
        Makes modificaitons to attribute reactants
        """
        if not isinstance(products,list)  and not isinstance(products,list):
            raise TypeError("Invalid 'products' format for reaction {}! Compounds must be a list of products but a {} object was provided instead".format(self.id, type(products))) 
        if len([n for n in products if not isinstance(n,compound.compound)]) > 0:
            raise TypeError("Invalid 'products' format for reaction {}! Compounds must be a list of 'compound' object but objects of {} were observed in the list instead. ".format(self.id, list(set([type(n) for n in products if not isinstance(n,compound.compound)]))))

        problem_cpds = [cpd.id for cpd in products if cpd not in [c for c in self.stoichiometry.keys() if self.stoichiometry[c] > 0]]
        if len(problem_cpds) > 0:
            raise userError('The following {} compounds appear in reactants of reaction {} but they do not appear in the reaction stoichiometry: {}'.format(len(problem_cpds), self.id, problem_cpds))

        if isinstance(products,tuple):
            self.__dict__['products'] = products
        elif isinstance(products,list):
            self.__dict__['products'] = list(products)

    def set_product_reactions(self, product_reactions):
        """
        Makes changes to attribute product_reactions
        """
        if product_reactions is not None and not isinstance(product_reactions,tuple) and not isinstance(product_reactions,list):
            raise TypeError("Invalid 'product_reactions' for compound " + str(self.id) + "! 'reactant_reactions'  must be a tuple or list of objects of type reaction. A " + str(type(product_reactions)) + " type object was entered instead")
        if len([n for n in product_reactions if not isinstance(n,reaction.reaction)]) > 0:
            raise TypeError("Invalid 'product_reactions' for compound " + str(self.id) + "! 'product_reactions'  must be a list of objects of type reaction. Objects that are not of type reaction found in the list: " + str([n for n in product_reactions if not isinstance(n,reaction.reaction)]))

        problem_rxns = [r.id for r in product_reactions if self not in [c for c in r.stoichiometry.keys() if r.stoichiometry[c] > 0]]
        if len(problem_rxns) > 0: 
            raise userError('The following {} reactions appears in "product_reactions" of compound {} but it does not appear in the stoichioemtry of this reaction as a product {}'.format(len(problem_rxns), self,id, problem_rxns))

        product_reactions = sorted(product_reactions,key=lambda x:x.id)

        if isinstance(product_reactions,tuple):
            self.__dict__['product_reactions'] = product_reactions
        elif isinstance(product_reactions,list):
            self.__dict__['product_reactions'] = tuple(product_reactions)

 def _assignFluxBounds(self):
     """
     Assigns general bounds to fluxes based on the reaction type 
     """
     if self.type.lower() == 'irreversible':
         self.flux_bounds = [0,1000]
     elif self.type.lower() == 'reversible':
         self.flux_bounds = [-1000,1000]
     elif self.type.lower() == 'reversible_forward':
         self.flux_bounds = [0,1000]
     elif self.type.lower() == 'reversible_backward':
         self.flux_bounds = [0,1000]
     elif self.type.lower() == 'exchange':
         self.flux_bounds = [0,1000]
     elif self.type.lower() == 'exchange_forward':
         self.flux_bounds = [0,1000]
     elif self.type.lower() == 'exchange_backward':
         self.flux_bounds = [0,0]
     else:
         raise userError('Unknown reaciton type')

def master_func(start_pos = None, end_pos = None, max_biomass_percentages = range(0,101,5), max_pyrrolysine_percentages = range(0,101,5), aeration = 'aerobic', media_type = 'minimal',results_filename = '', stdout_msgs = True, warnings = True):
    """
    Creates the model for E. coli harboring the pathways producing the non-standard amino acids
    and assesses its impact on the cell growth 

    INPUTS:
    -------
                      start_pos: Start position of the array containing all possble cases to consider (see all_cases variable)
                        end_pos: End position of the array containing all possble cases to consider (see all_cases variable)
        max_biomass_percentages: A vector containing the percentages of max biomass flux to consider
    max_pyrrolysine_percentages: A vector containing the percentage of max EX_pyrlys flux to consider
                     media_type: Type of growth mediam: M9 (minimal(, LB (rich)
                       aeration: Type of aeration 'aerobic' or 'anaerobic'
    """
    print '\n------------- {} , {} ---------------\n'.format(media_type, aeration)

    if not isinstance(aeration,str):
        raise TypeError('aeration must a a string')
    elif aeration.lower() not in ['aerobic','anaerobic']:
        raise ValueError('Invalid aeration value: {}'.format(aeration))

    if not isinstance(media_type,str):
        raise TypeError('media_type must be a string')
    elif media_type.lower() not in ['minimal','rich']:
        raise ValueError('Invalid media_type value! Allowed choices are: [M9, LB]')

    # Model path
    model_path = '/usr2/postdoc/alizom/work/models/Escherichia_coli/iJO1366/'

    if media_type.lower() == 'minimal' and aeration.lower() == 'aerobic':
        fluxBounds_dict = {'EX_glc(e)':[-10,1000], 'EX_o2(e)':[-20,1000]}
        fluxBounds_filename = model_path + 'iJO1366_minimal_glucose_aerobic.py'
    elif media_type.lower() == 'minimal' and aeration.lower() == 'anaerobic':
        fluxBounds_dict = {'EX_glc(e)':[-10,1000]} 
        fluxBounds_filename = model_path + 'iJO1366_minimal_glucose_anaerobic.py'
    elif media_type.lower() == 'rich' and aeration.lower() == 'aerobic':
        fluxBounds_dict = {'EX_glc(e)':[-10,1000], 'EX_o2(e)':[-20,1000]} 
        fluxBounds_filename = model_path + 'iJO1366_rich_glucose_aerobic.py' 
    elif media_type.lower() == 'rich' and aeration.lower() == 'anaerobic':
        fluxBounds_dict = {'EX_glc(e)':[-10,1000]} 
        fluxBounds_filename = model_path + 'iJO1366_rich_glucose_anaerobic.py' 
    else:
        raise ValueError('Unknown media_type and/or aeraiton type: {}, {}'.format(media_type,aeration))

    model_organism = organism(id = 'Ecoli', name = 'Escherichia coli',domain = 'Bacteria', genus = 'Escherichia', species = 'coli', strain = 'MG1655')

    # Orignal iJo1266 model
    model = create_model(model_organism = model_organism, model_info = {'id':'iJO1366', 'sbml_filename':model_path + 'iJO1366_updated.xml', 'biomassrxn_id':'Ec_biomass_iJO1366_core_53p95M'}, growthMedium_fluxBounds = {'fluxBounds_filename':model_path + fluxBounds_filename, 'fluxBounds_dict': fluxBounds_dict}, stdout_msgs = True, warnings = True)

    # Add the nsAA biosynthesis pathways to the model
    Ecoli_nsAA_producing = add_nsAA_pathways(model = deepcopy(model), stdout_msgs = stdout_msgs, warnings = warnings) 

    # Find the max biomass flux for Ecoli_nsAA_producing 
    print '\n--- FBA after adding new pathways ---\n'
    set_specific_bounds(model = Ecoli_nsAA_producing, file_name = fluxBounds_filename, flux_bounds = fluxBounds_dict)

    Ecoli_nsAA_producing.fba(stdout_msgs = stdout_msgs)
    if Ecoli_nsAA_producing.fba_model.solution['exit_flag'] == 'globallyOptimal':
        max_biomass = Ecoli_nsAA_producing.fba_model.solution['objective_value']
        print 'max biomass = {} , Pyrrolysine exch flux = {}  , p-amino-phenylalanine exch. flux = {}'.format(max_biomass,Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_pyrlys_L(e)'],Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_paphe(e)'])
    else:
        raise userError('FBA to find max biomass flux for Ecoli_nsAA_producing was not solved to optimality')

    # Find the max pyrrolysine exchange reaction flux 
    print '\n--- FBA to find max pyrrolysine exchange reaction flux ---\n'
    for rxn in Ecoli_nsAA_producing.reactions:
        rxn.objective_coefficient = 0
    Ecoli_nsAA_producing.reactions_by_id['EX_pyrlys_L(e)'].objective_coefficient = 1
    set_specific_bounds(model = Ecoli_nsAA_producing, file_name = fluxBounds_filename, flux_bounds = fluxBounds_dict)

    Ecoli_nsAA_producing.fba(stdout_msgs = stdout_msgs)
    if Ecoli_nsAA_producing.fba_model.solution['exit_flag'] == 'globallyOptimal':
        max_EX_pyrlys =  Ecoli_nsAA_producing.fba_model.solution['objective_value']
        print 'biomass = {} , max Pyrrolysine exch flux = {}  , p-amino-phenylalanine exch. flux = {}'.format(Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes'][Ecoli_nsAA_producing.biomass_reaction.id],Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_pyrlys_L(e)'],Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_paphe(e)'])
    else:
        raise userError('FBA to find max exchange reaction flux for L-pyrrolysine was not solved to optimality')

    # Find the max p-amino-phenylalanine exchange reaction flux 
    print '\n--- FBA to find max p-amino-phenylalanine exchange reaction flux ---\n'
    for rxn in Ecoli_nsAA_producing.reactions:
        rxn.objective_coefficient = 0
    Ecoli_nsAA_producing.reactions_by_id['EX_paphe(e)'].objective_coefficient = 1
    set_specific_bounds(model = Ecoli_nsAA_producing, file_name = fluxBounds_filename, flux_bounds = fluxBounds_dict)

    Ecoli_nsAA_producing.fba(stdout_msgs = stdout_msgs)
    if Ecoli_nsAA_producing.fba_model.solution['exit_flag'] == 'globallyOptimal':
        max_EX_paphe =  Ecoli_nsAA_producing.fba_model.solution['objective_value']
        print 'biomass = {} , Pyrrolysine exch flux = {}  , max p-amino-phenylalanine exch. flux = {}'.format(Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes'][Ecoli_nsAA_producing.biomass_reaction.id],Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_pyrlys_L(e)'],Ecoli_nsAA_producing.fba_model.solution['opt_rxnFluxes']['EX_paphe(e)'])
    else:
        raise userError('FBA to find max exchange reaction flux for p-amino-phenylalanine was not solved to optimality')

    # Set all objective function coefficients to zero
    for rxn in Ecoli_nsAA_producing.reactions:
        rxn.objective_coefficient = 0

    if results_filename != '':
        with open(results_filename,'w') as f:
            f.write('results = {}\n')

    all_cases = [(max_biomass_percent, max_pyrlys_percent) for max_biomass_percent in max_biomass_percentages for max_pyrlys_percent in max_pyrrolysine_percentages] 
#.........这里部分代码省略.........